Numerical Modelling of Airflow Within and Above Forests and Forest Clearings Using Computational Fluid Dynamics Timothy James Phaneuf B.Sc. (Hons.), Nipissing University, 2000 Thesis Submitted in Partial Fulfilment of The Requirements for the Degree of Master of Science in Natural Resources and Environmental Studies (Environmental Science) The University of Northern British Columbia April 2009 © Timothy James Phaneuf, 2009 1*1 Library and Archives Canada Bibliotheque et Archives Canada Published Heritage Branch Direction du Patrimoine de I'edition 395 Wellington Street Ottawa ON K1A0N4 Canada 395, rue Wellington Ottawa ON K1A0N4 Canada Your file Votre reference ISBN: 978-0-494-48784-6 Our file Notre reference ISBN: 978-0-494-48784-6 NOTICE: The author has granted a nonexclusive license allowing Library and Archives Canada to reproduce, publish, archive, preserve, conserve, communicate to the public by telecommunication or on the Internet, loan, distribute and sell theses worldwide, for commercial or noncommercial purposes, in microform, paper, electronic and/or any other formats. AVIS: L'auteur a accorde une licence non exclusive permettant a la Bibliotheque et Archives Canada de reproduire, publier, archiver, sauvegarder, conserver, transmettre au public par telecommunication ou par Plntemet, prefer, distribuer et vendre des theses partout dans le monde, a des fins commerciales ou autres, sur support microforme, papier, electronique et/ou autres formats. The author retains copyright ownership and moral rights in this thesis. Neither the thesis nor substantial extracts from it may be printed or otherwise reproduced without the author's permission. L'auteur conserve la propriete du droit d'auteur et des droits moraux qui protege cette these. Ni la these ni des extraits substantiels de celle-ci ne doivent etre imprimes ou autrement reproduits sans son autorisation. In compliance with the Canadian Privacy Act some supporting forms may have been removed from this thesis. Conformement a la loi canadienne sur la protection de la vie privee, quelques formulaires secondaires ont ete enleves de cette these. While these forms may be included in the document page count, their removal does not represent any loss of content from the thesis. Bien que ces formulaires aient inclus dans la pagination, il n'y aura aucun contenu manquant. Canada Abstract The computational fluid dynamics program, FLUENT, was first tested to validate windtunnel measurements of a scaled 10 ha forest clearing in a two dimensional domain. A variety of domain and canopy configurations were examined along with processor settings. Validation of the CFD program produced excellent results for horizontal wind velocity. Conifer shaped tree elements for the forest stands performed well and similar to the more traditional way of representing forest canopies. Turbulent kinetic energy (TKE) values output by the program seem to over predict the values calculated by using wind tunnel statistics. Various sizes of forest clearings were simulated to determine the stress that would be experienced by a forest edge immediately downwind of a clearing. Shorter gaps (< 15 tree heights) seem to experience higher values of TKE over the downwind forest, compared to the stand upwind of the clearing; and lower stress values along the downwind forest edge. Larger gaps (>60 tree heights) saw higher stress values but TKE values no larger than those reported upwind of the clearing. From the stress values calculated from various input velocities and gap sizes, a new tool was produced which takes into account a sites endemic wind speed and canopy density to predict stress on forest edges downwind of clearings. 11 Table of Contents List of Tables ix List of Figures xi Acknowledgements xxi 1. Introduction 1 1.1. Thesis objectives 4 1.2. 5 Thesis layout 2. Methods 8 2.1. Validation of data sets 2.2. 8 2.1.1. Sicamous Creek research forest 9 2.1.2. UBC wind tunnel 10 2.1.2.1. UBC model forest 11 2.1.3. Instrumentation/Data collection 12 2.1.4. Validation of wind tunnel 13 Computational Fluid Dynamics 13 2.2.1. FLUENT 14 2.2.1.1. Continuity and momentum equations iii 15 2.2.1.2. Reynolds-averaged Navier-Stokes equations 15 2.2.1.3. Turbulence models 17 2.2.1.3.1. The standard K-e turbulence model 18 2.2.1.3.2. The RNG K-e turbulence model 20 2.2.1.3.3. The realizable K-S turbulence model 21 2.2.1.4. Method of solution 24 2.2.1.5. Post processing 25 2.2.2. Two dimensional domains 26 2.2.2.1. Forest stands 27 2.2.2.1.1. Trees 28 2.2.2.1.2. Blockstands 29 2.2.2.1.3. Clearing size 30 2.2.2.2. Mesh 30 2.2.2.2.1. Grid Refinement 31 2.2.2.2.2. Grid independence 33 2.2.2.3. Boundary and Continuum Types 33 2.2.2.3.1. Canopy properties 35 2.2.2.3.2. Grid Interface 38 2.2.3. Three dimensional domains 38 2.2.3.1. Forest stands 40 2.2.3.2. Mesh 40 2.2.3.3. Boundary and continuum types 41 2.3. Model evaluation methods 43 iv 2.3.1. Output collection for model validation 43 2.3.1.1. Filled isolines 44 2.3.1.2. Scatterplots 45 2.3.1.3. Vertical profiles 45 2.3.1.4. Quantitative analysis 45 2.3.1.5. Taylor diagrams 50 2.3.2. Evaluation of gap size on forest edges 3. Comparison and testing of turbulence models 50 65 3.1. Introduction 65 3.2. Results and Discussion 66 3.3. 3.2.1. Standard K-e turbulence model 67 3.2.2. Renormalized Group K-E turbulence model 70 3.2.3. Realizable K-S turbulence model 71 3.2.4. Combined summary of model performance 72 Conclusion 73 4. Comparison and testing of different canopies and domains 93 4.1. Introduction 93 4.2. Results 95 4.2.1. Porous tree stand 95 4.2.1.1. Qualitative description of profile 96 4.2.1.2. Scatterplot description 97 v 4.2.1.3. Velocity profile and description 98 4.2.1.4. Turbulent kinetic energy profile and description 99 4.2.1.5. Quantitative description 101 4.2.1.6. Discussion of porous tree stand 102 4.2.2. Two dimensional layered tree stand (Blockstand) 103 4.2.2.1. Qualitative description of profile 104 4.2.2.1.1. Standard K-e turbulence model 104 4.2.2.1.2. Renormalized Group K-e turbulence model 105 4.2.2.1.3. Realizable K-e turbulence model 106 4.2.2.2. Scatterplot description 106 4.2.2.2.1. Standard K-e turbulence model 107 4.2.2.2.2. Renormalized Group K-e turbulence model 107 4.2.2.2.3. Realizable K-e turbulence model 108 4.2.2.3. Vertical profile descriptions of horizontal velocity 108 4.2.2.4. Vertical profile descriptions of turbulent kinetic energy 110 4.2.2.5. Quantitative description Ill 4.2.2.6. Discussion of layered tree stand 112 4.2.3. Three dimensional layered blockstands 113 4.2.3.1. Qualitative description of profile 113 4.2.3.2. Scatterplot description 114 4.2.3.3. Velocity Profile and description 114 4.2.3.4. Turbulent kinetic energy profile and description 114 4.2.3.5. Quantitative description 115 vi 4.2.3.6. Discussion of the three dimensional layered stand 4.2.4. Layered tree stand (Blockstand) in 3 m tall domain 115 116 4.2.4.1. Vertical profiles of horizontal velocity between two domain heights 116 4.2.4.2. Vertical profiles of TKE between the standard and 3 m domain 118 4.2.4.3. Quantitative comparison 119 4.2.4.4. Discussion of 3m tall domain 120 4.2.5. Miscellaneous domain/canopy simulations 120 4.2.5.1. Stands only, no roughness elements with and without modified grid 121 4.2.5.2. Discussion of miscellaneous domain/canopy simulations 4.3. Discussion 4.4. 122 122 4.3.1. Porous tree stand 124 4.3.2. Layered tree stand (Blockstand) 124 4.3.3. Comparison of flow traits with other studies 127 4.3.4. Exit flows 128 Conclusion 130 5. Effects of gap size on forest edge stress 163 5.1. Introduction 163 5.2. Setup 165 5.3. Results and discussion 166 vii 5.3.1. Horizontal wind velocity 166 5.3.2. Turbulent kinetic energy 168 5.3.3. Edge stress on downwind stand 169 5.3.3.1. Logarithmic fits 170 5.3.4. Horizontal velocity profile 171 5.4. Conclusion 172 6. Summary of conclusions 211 6.1. Comparison and testing of K-e turbulence models 211 6.2. Comparison and testing of different canopies and domains 212 6.3. Effects of gap size on forest edge stress 214 6.4. Reflections and recommendations for future work 215 7. Nomenclature 217 8. Bibliography 220 vm List of Tables Table 2.1: Height ranges and accompanying leaf area densities for the different layers within the porous blockstand canopy, as reported by Liu et al. (1996) in their wind tunnel stand which was the same to the UBC validation data 58 Table 2.2: Clearing sizes used to determine stress at the downwind edge of a clearing between two porous blockstand forest stands 59 Table 3.1: Statistics for horizontal velocity for both short and long iterated steady state simulations for each of the K-e turbulence models, along with unsteady 4.5 and 6.0 second simulations using the standard (std) K-e turbulence model. Jn the table below, O is the averaged wind speed from all wind tunnel observations. P is the average wind speed output from sampled points from the simulation. The standard deviation for the observed data is noted by s0 and sp for the model output. N is the number of sampled points. The y-intercept is represented by b and m is the slope of the best fit line. MAE stands for Mean Absolute Error. RMSE is the Root Mean-Square Error with RMSES and RSMEU being the systematic and unsystematic components of the RSME, respectively. The Willmott d score, ranging from 0 to 1, is indicated by d, and finally, the correlation coefficient is r2. 81 Table 4.1: Statistics of horizontal velocity for the porous tree stand simulation using three versions of the K-e turbulence model._In the table below, O is the averaged wind speed from all wind tunnel observations. P is the average wind speed output from sampled points from the simulation. The standard deviation for the observed data is noted by s0 and sp for the model output. TV is the number of sampled points. The y-intercept is represented by b and m is the slope of the best fit line. MAE stands for Mean Absolute Error. RMSE is the Root Mean-Square Error with RMSES and RSMEU being the systematic and unsystematic components of the RSME, respectively. The Willmott d score, ranging from 0 to 1, is indicated by d, and finally, the correlation coefficient is r2 138 Table 4.2: Statistics of turbulent kinetic energy for the porous tree stand simulation using three versions of the K-e turbulence model (variables defined in table 4.1) 138 Table 4.3: Statistics of horizontal velocity for three K-e turbulence model schemes (standard, realizable and RNG) in the porous blockstand domain (variables defined in table 4.1). 149 Table 4.4: Statistics of turbulent kinetic energy for three K-e turbulence model schemes (standard, realizable and RNG) porous blockstand domain (variables defined in table 4.1). 149 IX Table 4.5: Statistics of horizontal velocity for two domain setups using the standard K-S turbulence model (two dimensional blockstand and three dimensional blockstand) (variables defined in table 4.1) 155 Table 4.6: Statistics of turbulent kinetic energy for two domain setups using the standard K-S turbulence model (two dimensional blockstand and three dimensional blockstand) (variables defined in table 4.1) 155 Table 4.7: Statistics of horizontal velocity for two domain setups using the standard K-S turbulence model (standard height and 3 m tall domain) (variables defined in table 4.1). 158 Table 4.8: Statistics of turbulent kinetic energy for two domain setups using the standard K-S turbulence model (standard height and 3 m tall domain) (variables defined in table 4.1). 158 Table 4.9: Statistics of horizontal velocity for two layouts, one with no roughness elements or bluff bodies using a blockstand style canopy (all canopy layout) and the other similar but utilising a boundary layer (BL) grid, using the standard K-e turbulence model (standard layout and all blockstand layout) (variables defined in table 4.1) 159 Table 4.10: Statistics of turbulent kinetic energy for two layouts, explained above, using the standard K-£ turbulence model (standard layout and all blockstand layout) (variables defined in table 4.1) 159 Table 4.11: Quantitative statistics of horizontal velocity from the full and bottom half of the analysed profile for the porous tree stand simulation using three versions of the K-e turbulence model (variables defined in table 4.1) 161 Table 4.12: Quantitative statistics of turbulent kinetic energy from the full and bottom half of the analysed profile for the porous tree stand simulation using three versions of the K-e turbulence model (variables defined in table 4.1) 161 Table 4.13: Quantitative statistics of horizontal velocity from the full and bottom half of the analysed profile for the porous blockstand simulation using three versions of the K-S turbulence model (variables defined in table 4.1) 162 Table 4.14: Quantitative statistics of turbulent kinetic energy from the full and bottom half of the analysed profile for the porous blockstand simulation using three versions of the K-e turbulence model (variables defined in table 4.1) 162 x List of Figures Figure 2.1: Sicamous Creek research forest was established by the British Columbia Ministry of Forests and was the site of many studies into forest health. Observations collected at the 10 ha plot, B-5, were used to validate wind tunnel measurements collected by the UBC team (Image from Puttonen and Murphy, 1997. Used with permission, British Columbia Ministry of Forests and Range, Research Branch) 52 Figure 2.2: Aerial view of B-5 plot. The Engelmann spruce and subalpine fir trees surrounding the 10 ha clearcut are approximately 30 m in height (Photo from Novak et al., 2001). 52 Figure 2.3: A detailed schematic of the upwind section of the wind tunnel used by Dr. Novak and his research team at UBC. Spires, bluff bodies, large roughness elements and small roughness elements are shown 53 Figure 2.4: Array of five Counihan spires used for generating large scale turbulence in the wind tunnel used by Dr. Novak and his team at UBC (Photo from Novak et al., 2001). . . . 54 Figure 2.5: Dimensions of the large and small roughness elements used in the wind tunnel. 54 Figure 2.6: Approximate dimensions of each tree in the wind tunnel 55 Figure 2.7: Artificial Christmas tree branches were used to make up the model forest in the wind tunnel (Photo from Novak et al., 2001) 55 Figure 2.8: RM Young Anemometer in plot B-5 at the Sicamous Creek research forest (Photo from Novak et al., 2001) 55 Figure 2.9: Measurements collected from both Sicamous field site and UBC wind tunnel (Image from Novak et al., 2001) 55 Figure 2.10: The two main types of numerical domains used for the simulations. A) Multiple tree shaped entities are used to make up the forested areas. B) Layered rectangular blocks represent the forest stands (blockstands) 56 Figure 2.11: Three metre tall domain used to help explain higher wind speeds seen in the 1.5 m tall domain. This domain also used the porous blockstand canopy 57 Figure 2.12: The "stand only" domain with a porous blockstand canopy and no upwind roughness elements was used to examine a domain with no upwind bluff bodies or XI roughness elements 57 Figure 2.13: The two dimensional numerical trees were based upon the silhouetted shape of the wind tunnel trees and had similar dimensions 58 Figure 2.14: Portion of the domain showing detail of the mesh along the upwind stand. Grey lines extending vertically are locations where sampling took place for validation with wind tunnel measurements. A) Mesh shown for the standard sized mesh used (adapted). Each grid segment of the triangular grid <2 z/ht is no larger than 1/10 ht or 0.015 m in length. B) To help determine grid independence a much finer mesh was produced, also a triangular grid, but with no grid segment larger than 1/30 ht or 0.005 m in length, producing a mesh nine times finer 60 Figure 2.15: Hanging node grid adaption used by FLUENT leaves a grid node unused on one side (after ANSYS, 2006) 61 Figure 2.16: Before (left) and after (right) illustrations of hanging node grid adaption applied to a triangular grid cell (after ANSYS, 2006) 61 Figure 2.17: This figure shows an example of conformal grid adaption used in FLUENT. Cell A is marked for refinement along it's longest edge. Because Cell B also shares the edge, it too will be split (after ANSYS, 2006) 61 Figure 2.18: Vertical profiles of horizontal velocity from both the standard grid ( ) and fine (GI) grid ( ) for determining grid independence. These profiles are shown with measurements from the wind tunnel (•) for comparison and will be described in more detail in subsequent chapters. The vertical profiles shown are from locations within the domain, with location "a" being the most upwind location sampled and location " j " the most downwind location sampled (see figure 2.21 to see where the specific positions of the sampling locations are within the wind tunnel/domain 62 Figure 2.19: A) Two continua cannot occupy the same space, so a "subtraction" and "retain" function must be used so that the two separate entities can exist. Note: white area only shows separation between two boundary conditions and is not actually present. B) By using a two sided boundary condition, identification and organisation of special boundary conditions or continua becomes easier and more manageable 63 Figure 2.20: An example of an offset grid interface between two zones. Interface zone 1 is made up of two faces, A-B and B-C, with interface zone 2 made up of faces D-E and E-F. The intersection of these faces produce new faces: a-d, d-b, b-e and e-c . Face a-d will form a wall zone, while faces d-b, b-e and e-c will form two-sided interior faces. To calculate the flow into cell IV, both faces d-b and b-e will be used to bring in the information from cell I and III, respectively (after ANSYS, 2006) 63 Figure 2.21: Positions of the ten sampling locations within the stands and clearing of both the xn wind tunnel and numerical domains. Moving left to right, Location A is the most upwind location, Locations B through I denoting the next eight locations with Location J being the most downwind location. Numbers indicate the Location's position (in metres) relative to the upwind edge of the clearing 64 Figure 3.1: Isotachs of horizontal velocity (flowing from left to right) from the shorter iterated standard K-S turbulence model are presented. The grey vertical lines indicate the ten sampling locations (A to J - refer back to figure 2.21 for explanation) where simulation output can be more directly compared to wind tunnel data (vertical profiles or quantitative analysis). White areas, exterior to the canopy elements, indicate reverse flow (negative horizontal velocity) 74 Figure 3.2: Filled isotachs of the converged, longer iterated K-e turbulence model simulation have a slightly smoother profile over the forest unlike the isotach intervals from the shorter one (figure 3.1). This profile compares much better with the isotachs from the wind tunnel (figure 3.5 a). As in figure 3.1, white areas exterior to the canopy elements indicate a reverse horizontal velocity 75 Figure 3.3: Turbulent kinetic energy from the lower iterated K-S turbulence model simulation. White areas within the canopy indicate that no domain is present and therefore no flow is permitted through the canopy. A large area of TKE is seen over the stands and clearing. This pattern seems unrealistic and is strong evidence that the simulation needs to be iterated longer 76 Figure 3.4: Similar to figure 3.3, this figure shows TKE from the longer iterated simulation using the standard K-e turbulence model. The TKE from this simulation is much lower above the clearing and stands, but notably higher above the upwind edge of the second stand 77 Figure 3.5: Observations from the wind tunnel were interpolated, plotted and isolines drawn from the data. In the plots above, the most upwind extent (-12 x/ht) would represent location "A" in figure 2.21; the most downwind extent (+14.25 x/ht) would represent location "J". Dashed lines show the extent of the canopy. A) Isotach of horizontal velocity show that the flow over the upwind stand is relatively streamlined with a tighter gradient close to the top of the forest. Above the upwind portion of the clearing wind speed increases, this increase becomes more pronounced with increasing height. Ahead of the second stand, there is a sudden decrease in horizontal velocity. The intensity of the decrease above the height of the canopy is mitigated with increasing height. Above the second stand, the speed increases as the flow adjusts to the new surface conditions. Also in the second stand, stem flow can be seen close to the ground. B) Isolines of turbulent kinetic energy calculated from wind tunnel measurements yield an interesting pattern. The area of greatest turbulence values is seen above and downwind of the first stand. There is also tight gradient that bounds the forest perimeters 78 Figure 3.6: Scatterplots of the shorter iterated (a) and long iterated (b) simulations show xni inconsistent trends at lower speeds, which were likely from those points extracted within or just above the canopy. The scatterplot with the shorter iterated simulation shows a fair degree of scatter at the higher speeds, while the converged simulation has a tighter profile at the faster speeds 79 Figure 3.7: Vertical profile plots of horizontal velocity for short and long simulations of all three K-S turbulence models from the domain with the solid tree shaped elements, as well as wind tunnel observations. Positions of each location (a-j) are in relation to the upwind edge of the clearing and can be visually referenced in figure 2.21 80 Figure 3.8: Isotachs of horizontal velocity. Similar setup to the simulation seen in figure 3.1 with the exception that this output was from the shorter iterated simulation using the RNG K-e turbulence model 82 Figure 3.9: Isotachs of horizontal velocity, this image is similar to figure 3.8 (RNG K-S turbulence model used) except that it was simulated for more iterations compared to the previous figure 83 Figure 3.10: Turbulent kinetic energy isolines from the lower iterated RNG K-£ turbulence model simulation. Aside from the turbulence model used, the setup for this simulation was very similar to the simulation used to generate figure 3.3 84 Figure 3.11: The main difference between the simulation from which this figure was produced from and the simulation from figure 3.10 was the number of iterations. This simulation was also used the RNG K-e turbulence model. Isolines of turbulent kinetic energy are shown 85 Figure 3.12: Scatterplots of the lesser iterated (a) and long iterated (b) RNG K-S turbulence model simulations appear fairly similar to one another. A little less scattering can be seen in the longer iterated plot 86 Figure 3.13: Similar to figure 3.1, isotachs of horizontal velocity for the lower iterated simulation are shown except that this image was from the simulation using the realizable K-S turbulence model 87 Figure 3.14: Side profile of horizontal velocity isotachs from the longer iterated realizable K-S turbulence model simulation, refer to figure 3.2 for additional details regarding setup. White areas exterior to the canopy indicate areas of reverse flow 88 Figure 3.15: Similar to figure 3.3 except that the turbulent kinetic energy isolines from in this image are from the lower iterated realizable K-S turbulence model simulation 89 Figure 3.16: Similar to figure 3.15 except that turbulent kinetic energy isolines are from the longer iterated realizable K-E turbulence model simulation 90 xiv Figure 3.17: Scatterplots of the short iterated (a) and long iterated (b) simulations using the realizable K-S turbulence model. These two plots exhibit similar trends seen in the other two turbulence models. The two plots are fairly similar, except for a tighter clustering of points at the higher windspeeds 91 Figure 3.18: Taylor diagram of all six simulations, low and high iterated simulations from each of the three K-e turbulence models (standard, RNG and realizable). Quantitative analysis from table 3.1 provided the data to graphically show how all six simulations relate to one another and to a perfect simulation (Ref.) 92 Figure 4.1: Side view of the porous tree shaped stands and clearing showing isotachs (ms 1 ) from the simulation using standard K-8 turbulence model with an input velocity of 8.7 ms"1. Vertical lines indicate sampling locations (left to right) A to J (refer to figure 2.21 for more details) 133 Figure 4.2: Side view of porous tree stands and clearing, similar to figure 4.1 except showing isolines of turbulent kinetic energy (mV 2 ). The standard K-£ turbulence model was also used for this simulaiton 134 Figure 4.3: Scatterplot distribution of horizontal velocity from all sampling points in both the wind tunnel and porous tree domain simulated using the standard K-8 turbulence model with an input velocity of 8.7 ms"1 135 Figure 4.4: Vertical profiles of horizontal wind speed from both wind tunnel observations and output from porous tree elements simulation using the three variants of the K-8 turbulence models. Locations a and b are within the first stand, upwind of the clearing; c, d and e are from the first half of the clearing; f, g, and h are locations from the second half of the clearing and i and j are locations found in the stand downwind of the clearing. Locations of each plot in the legend are from the upwind edge of the clearing (refer to figure 2.21). 136 Figure 4.5: Vertical profiles of turbulence kinetic energy calculated from wind tunnel observations and output from the porous tree elements simulation using the standard K-E turbulence model and two variants. Letters in top left of each plot reference the location of that plot as demonstrated in figure 2.21 137 Figure 4.6: Side profile view of layered blockstands showing isotachs from the standard K-8 turbulence model with an input velocity of 8.7 ms"1. Vertical grey lines represent the sampling locations (A-J) as defined in figure 2.21. Stem flow evidence can be seen in the second stand downwind of the forest edge near the ground 139 Figure 4.7: Simulation output from the standard K-8 turbulence model. Side profile view of the layered canopy domain showing the stands and clearing. Isolines of turbulent kinetic energy show areas of increased turbulence are found just above the tops of the canopy and close to the leeward end of the stands. TKE above the second stand is greater than what xv 140 was simulated above the first stand. Figure 4.8: Side profile showing isotachs of horizontal velocity from the porous blockstand simulation much like the one in figure 4.6 but instead using the RNG K-S turbulence model 141 Figure 4.9: The simulation from figure 4.8 except showing isolines of turbulent kinetic energy. TKE values over both stands seem to be fairly equal, with a small patch of TKE >6 m 2 s 2 over the second stand 142 Figure 4.10: Similar to figure 4.6, horizontal velocity isotachs from the simulation using the realizable K-e model in the porous blockstand domain 143 Figure 4.11: Turbulent kinetic energy isolines from the layered domain simulated using the realizable K-e turbulence model and an input velocity of 8.7 ms"1. Greater TKE was predicted over the second stand 144 Figure 4.12: Scatterplot distribution of horizontal velocity from observed and modelled sampling points for two unsteady simulations run for 4 seconds (a) and 6 seconds (b). Both simulations were identical in every other way using the porous blockstand canopies, the standard K-S turbulence model and an input velocity of 8.7 ms"1. These plots show virtually no difference in output for the longer iterated simulation 145 Figure 4.13: Scatterplot distribution for observed wind tunnel values versus modelled output. Aside from the use of the RNG K-e turbulence model and the time of simulation (4.6 seconds) all other settings were identical to those using in figure 4.12 145 Figure 4.14: Scatterplot distribution of horizontal velocity from observed wind tunnel data and output from the porous blockstand canopy domain using a realizable K-e turbulence model and an input velocity of 8.7 ms"1. No difference is seen between the four (a) and 4.6 (b) second simulations 146 Figure 4.15: Vertical profiles of horizontal wind speed from the wind tunnel observations and output from three versions of the K-S turbulence model. The three versions of the K-e turbulence models are: the standard (std) model, the realizable (real) model and the renormalized group (RNG) model. The ten plots represent the ten sampling locations (a-j) as shown in figure 2.21. Positions mentioned in the legend are from the upwind edge of the clearing 147 Figure 4.16: Vertical profiles of turbulent kinetic energy. Wind tunnel values were calculated using wind statistics from the wind tunnel observations. Output from the porous blockstand domain includes: standard (std), realizable (real) and renormalized group (RNG) K-e turbulence models. Each plot is identified by a letter which is also the location where the data and output were sampled from in the wind tunnel or domain, respectively (refer to figure 2.21) 148 xvi Figure 4.17: Filled horizontal isotachs from the simulation of the three dimensional porous blockstand domain using the standard K-S turbulence model. Input velocity was 8.7 ms"1. 150 Figure 4.18: From the same simulation in figure 4.17, isolines of turbulent kinetic energy are displayed. A larger area of TKE was predicted in the downwind stand 151 Figure 4.19: Scatterplot of horizontal wind speed, observed results measured in the wind tunnel and predicted output from the three dimensional porous blockstand domain simulated using the K-e turbulence model 152 Figure 4.20: Vertical profiles of horizontal wind velocity from the wind tunnel along with output from the three dimensional porous blockstand domain using the standard K-S turbulence model. The two dimensional porous blockstand output (from figure 4.15) is also included for comparison. Letters at the top right of each plot indicate the locations as shown in figure 2.21 153 Figure 4.21: Vertical profiles of turbulent kinetic energy calculated from wind tunnel observations. TKE profiles from the three dimensional porous blockstand domain are plotted along with output from the two dimensional porous blockstand simulations, both using the standard K-S turbulence model. There is a large difference between simulation output and those results calculated from the wind tunnel 154 Figure 4.22: Vertical profiles of horizontal velocity from the ten sampling locations (as shown in figure 2.21) from the wind tunnel (•) and output from the standard K-e turbulence model from the standard height domain ( ) and 3 m tall domain ( ) both using a porous blockstand canopy 156 Figure 4.23: Vertical profile for turbulent kinetic energy (TKE). Plots show wind tunnel results calculated from wind statistics (•) along with output from the standard height domain ( ) and the 3 m tall domain ( ). Both simulations were run using the standard K-S turbulence model and used a porous blockstand canopy to represent the forest stands. 157 Figure 4.24: Taylor diagram of all ten simulations done in this chapter and how they all relate to one another and to the "perfect simulation" or wind tunnel (Ref.) 160 Figure 5.1: Exit flow across a two-dimensional forest edge. The approach flow (A) upwind of the forest edge and at the forest edge (Al). Downwind of the stand, the flow field in the open is divided into quiet zone (D), mixing zone, (E) and re-equilibration zone (F) (after Lee, 2000) 173 Figure 5.2 a: Horizontal velocity isotachs of airflow over two porous blockstands separated by a 4 x/ht clearing and simulated with an input velocity of 8.7 ms"1 and using the standard K-S turbulence model 174 XVI1 Figure 5.2 b: Horizontal isotachs through a canopy with an 8 x/ht clearing 175 Figure 5.2 c: Clearing of 10 x/ht showing an increase in horizontal wind speed penetration towards the ground 176 Figure 5.2 d: Horizontal velocity isotachs showing a 10.9 x/ht clearing, which at full scale (l:200h t )would be just over 10 ha 177 Figure 5.2 e: Forest gap of 14 x/ht, 3 m/s isotach is well entrained along the ground . . . . 178 Figure 5.2 f: Clearing representing a 20 tree height gap 179 Figure 5.2 g and h: Increased clearing sizes clearly shows the increased infiltration between the 30 ht (g) and 50 h, (h) gap size. In g, the 6 ms"1 isotach is still above the height of the forest, in h, the isotach has moved down to just above the surface 180 Figure 5.2 i,j and k: Sixty, 80 and 100 ht. I,j and kail show the progression of the higher wind speeds infiltrating further down towards the surface 181 Figure 5.3 a: Horizontal isotachs flowing over the 4 x/ht gap size from the 4.35 ms"1 (input velocity) simulation also using the porous blockstand canopy and standard K-e turbulence model 182 Figure 5.3 b: Clearing size of 10 ht 183 Figure 5.3 c: Clearing size of 30 h, 184 Figure 5.3 d and e: Top image (d) shows a gap length of 60 ht, while the bottom (e) is 100 ht. It appears that the flow has fully adjusted to the new surface in figure e 185 Figure 5.4 a: Isotachs of a cross section with a 4 h, gap in a simulation with a 17.4 ms"1 input horizontal velocity using the standard K-e turbulence model and the porous blockstand canopy 186 Figure 5.4 b: Isotachs from the 10 ht gap domain 187 Figure 5.4 c: A 30 ht clearing with isotachs 188 Figure 5.4 d and e: top image (d) shows isotachs of the stands and 60 h, clearing; bottom image (e) shows a 100 ht clearing 189 Figure 5.5 a: Filled isolines showing turbulent kinetic energy from the domain with a 4 ht gap between two porous blockstands and an input velocity of 8.7 ms"1, in a simulation using the standard K-e turbulence model. Higher values of TKE are predicted over the second xvin forest stand, downwind of the clearing 190 Figure 5.5 b: Isolines of TKE in the domain with the 8 ht clearing 191 Figure 5.5 c: A 10 ht clearing and the resulting TKE isolines in and above the clearing and stands 192 Figure 5.5 d: Filled TKE isolines in the domain with the 10.9 h, clearing 193 Figure 5.5 e: A clearing of 14 h, and the resulting filled isolines of TKE 194 Figure 5.5 f: Clearing of 20 ht with the second stand just out of range of this figure, showing TKE 195 Figure 5.5 g and h: TKE shown in the 30 ht (g) and 50 ht (h) clearings 196 Figure 5.5 i, j and k: Domains, with gap sizes of 60 ht (i), 80 ht (j) and 100 ht (k), showing TKE. Figures i and j shows that the maximum TKE value over both stands is similar. Figure k shows that the peak TKE value above the second (downwind) stand is lower than what was predicted in the first (upwind) stand 197 Figure 5.6 a: Filled isolines of TKE from the simulations with the 4.35 ms"1 input velocity, this figure showing the domain with the 4 ht clearing using the porous blockstands and simulated using the standard K-e turbulence model. Unlike the other simulation with an 8.7 ms"' input velocity, these two stands have similar maximum values over upwind and downwind stands 198 Figure 5.6 b: TKE in and above two stands and a clearing 10 ht in length 199 Figure 5.6 c: Filled isolines showing the distribution of TKE in and adjacent to a 30 ht clearing. 200 Figure 5.6 d and e: TKE drops below 0.5 m2s"2 at around 40 ht in both the 60 ht (d) domain and 100 ht (e) domain. Even at these large gap sizes, the upwind and downwind stands have similar peak TKE values 201 Figure 5.7 a: Filled isolines of turbulent kinetic energy from the porous blockstand canopy domain with the 4 ht gap. This simulation had the inlet velocity set at 17.4 ms'1 and was simulated using the standard K-£ turbulence model. The peak TKE seen and the extent of it over the second stand is noticeably higher than the first stand. Quite high values of TKE are seen between stands above the clearing 202 Figure 5.7 b: TKE from the domain with the 10 ht gap xix 203 Figure 5.7 c: Gap from the 30 ht domain with isolines of TKE 204 Figure 5.7 d and e: Turbulent kinetic energy is still present far downwind of the first stand in both 60 ht (d) and 100 h, (e) clearings 205 Figure 5.8: Normalised stress values on the downwind edge of various gap sizes for three input velocities 206 Figure 5.9: Normalised stress from each input velocity versus the time taken to cover each gap distance 207 Figure 5.10: Normalised stress from each clearing size and each input velocity simulated as a function of time required to flow from start to end of clearing 208 Figure 5.11: Horizontal velocity profile of the 100 x/ht clearings at 1/3 stand height between stands 209 Figure 5.12: Similar to figure 5.11, but with horizontal velocity profile plotted against gap transit time 210 xx Acknowledgements I would first like to thank my Mother and Father for all their love and support. Many thanks to the members of my thesis advisory committee. My advisor, Peter L. Jackson, who with great patience, stuck with me to see a completed thesis. You Qin (Jean) Wang, who provided priceless support with her CFD and FLUENT knowledge. Michael D. Novak, who provided all of my validation data; my thesis would have been vastly different without your invaluable contribution. I also extend thanks to Ian D. Hartley, who made certain that I would (finally) have the opportunity to defend this semester. I would like to thank Stephen J. Mitchell, for agreeing to be the external examiner. Finally, I would like to thank all of the friends I've made here in Prince George and at UNBC, you all have enriched my time at UNBC. xxi 1. Introduction Forests make up a large area of Canada, covering 402.1 million ha and represent 29 billion m3 of total growing stock (CCFM, 2006). The abundance of this resource makes forestry an important industry to Canada and in particular, British Columbia. In 1992, the forest sector in BC was worth $11.5 billion annually to the provincial economy, representing $8000 for every household in British Columbia (BC Ministry of Forests, 1992). Natural factors that affect the health and productivity of forests include insect pests, forest fires and wind damage. While forest fires and more recently insect infestations garner most of the attention, wind damage losses are an important aspect to consider. Prior to the recent Mountain Pine Beetle epidemic, the average amount of timber damaged annually in British Columbia by wildfire and insects was over 5 million m3 and 6.5 million m3, respectively (BC Ministry of Forests, 1992). Wind damage was over 3 million m3 or equivalent to 4 % of the annual allowable cut on crown land in British Columbia (Mitchell, 1995). In forestry, wind is studied so that the impact on forested environments can be better understood and predicted. Wind - tree interactions can be catastrophic and lead to either a bole (stem) break or anchoring failure of trees. This failure of trees is referred to as windthrow. A great deal of literature on windthrow examines wind and tree interactions. Windthrow can be described as either catastrophic or endemic (Rollinson, 1987; Quine, 1995 and Stathers et al., 1994). The winds causing catastrophic windthrow are usually caused by exceptional storms that have very high winds (gusts > 40 ms"1) and are infrequently recorded 1 in an area (Quine, 1995). The winds that cause endemic windthrow are the annual peak winds observed in a region with gusts around 30 ms"1 (Quine, 1995; Quine and Bell, 1998). It is these endemic winds that cause more damage to timber supplies than catastrophic winds (20 % versus 6 %, respectively) (Rollinson, 1987 and Atterson, 1980 as cited in Rollinson, 1987). Harvesting methods can affect wind firmness in a forest. The two principle methods of harvesting timber are clearcuts, which result in cut edge boundaries, and thinning or selective harvesting of a stand. In British Columbia, 181 530 ha of forests were harvested in 1990, of this 166 565 ha were harvested by clearcutting (91.8%) while 14 965 ha were selectively harvested (8.2%) (BC Ministry of Forests, 1992). Over a 26 year period from 1970 to 1996, the proportion of harvesting by clearcutting and partial cutting systems averaged 87% and 13%, respectively (BC Ministry of Forests and Range, 2006). In the late 1990's, due to public pressure and new scientific findings, a new modified method of clearcutting was introduced called clearcutting with reserves (BC Ministry of Forests and Range, 2006). Clearcutting with reserves is a slight modification of clearcutting that retains a variable number of trees within the cutblock, either alone or in small groups for purposes other than regeneration (BC Ministry of Forests, 1995). In 2006, approximately 88 % of harvesting was with clearcutting: 58 % clearcutting with reserve and approximately 30% was standard clearcutting (BC Ministry of Forests and Range, 2006). Size of clearcut areas have also decreased since the late 1990's. While the changes to cutblock design may have had certain public and ecological benefits, they do pose a problem with wind firmness. To have the same yield from smaller and reserved cutblock designs, more cutblocks are required, which increases the amount of forest edges exposed to wind. As the amount of cutblock edges increase, so to do the concerns surrounding windthrow as recently harvested forest edges are more prone to windthrow (eg. Stacey et al., 1994; Ruel et al. 2001). Rowan et 2 al. (2003) and Gardiner et al. (1997) have studied ways of reducing damage along forest edges. Gardiner et al. (1997) investigated thinning in potential harvest areas which would allow trees to increase their taper and root mass. Rowan et al. (2003) suggest that proper pruning/topping of the crown can help alleviate wind damage along newly exposed forest edges. Other researchers have approached windthrow from many different perspectives: physiological, topographical, meteorological andpedological, using field measurements as well as wind tunnels and numerical simulations. Field studies are important sources of real world data, but data are often expensive to collect and it is hard to control conditions (Meroney, 1968). Wind tunnels can be less costly and data can be collected at many possible points and conditions can be better controlled than in a field. Wind tunnels have many different uses; however, they have limited availability. While wind tunnel studies may be less expensive than field studies, the cost of performing numerical studies has been rapidly decreasing due to progressively more powerful computers and improved software that can be used in many fields of study (Anderson, 1996). Numerical simulations can be run on a wide range of platforms, from personal computers to powerful compute servers. Three early and notable, numerical simulations performed using forested landscapes included Li et al. (1990), Green (1992) and Liu et al. (1996). All predicted airflow through canopies, with Li et al. (1990) and Liu et al. (1996) simulating airflow at a forest edge while Green (1992) examined airflow through a forest of widely spaced trees. The study by Li et al. (1990) reproduced Raynor's (1971) observations of flow into and out of a forest, while the study by Liu et al. (1996) attempted to reproduce wind tunnel measurements with airflow moving from a forest to a clearing. Green's (1992) simulation was also compared to wind tunnel measurements. 3 In the present study, horizontal airflow and turbulence were numerically simulated using a commercial computational fluid dynamics (CFD) package. This software package, FLUENT® (by ANSYS, Inc. of Canonsburg, Pennsylvania, USA) hereafter referred to as FLUENT, is a robust suite of tools that can perform all aspects involved with numerical modelling, from domain creation and meshing to simulation and post processing. Such programs have the ability to recreate real-world properties and simulate the characteristics of fluids. Much like water in a river, air flows over the landscape and interacts with elements in that landscape. Slight topographic and vegetative variations in these landscapes are able to affect airflow on a local scale. The majority of numerical studies use specialised code specifically developed for modelling airflow through and above forests, but the present study will be using a commercially available CFD package for all of the simulations. Turbulence will be handled by a specialised model built into the CFD package that attempts to resolve turbulence at scales smaller than the grid allows for. There are three different K-S turbulence models that were tested during the course of this thesis. One version is the original, standard version of the K-£ turbulence model (Jones and Lauder, 1972), while two others are refinements of the original. 1.1. Thesis objectives The main objective of this thesis is to determine whether a commercial CFD program like FLUENT provides an adequate representation of real world conditions in and above a forest stand and cutblock. The specific research objectives are: • To assess the use of FLUENT, in reproducing airflow conditions measured within and 4 above a model forest stand and cutblock. • To assess the K-e turbulence model and two of it's derivatives used by the FLUENT processor. • To assess the use of porous tree elements as a new way of representing forest stands in numerical simulations. • 1.2. To develop a means of predicting stress on a forest edge downwind of a clearing. Thesis layout Following this introductory first chapter, the second chapter covers the methodology, with the three chapters following detailing the results of the thesis objectives. The concluding chapter will recap the major findings from this thesis. While this thesis resembles a classic style thesis there are some noticeable differences. Instead of a literature review section, pertinent literature is reviewed in each chapter, in relation to the topic being discussed. Also, within each chapter, figures are inserted at the end for ease of reference, since several of them are referred to multiple times in different parts of the thesis. The second chapter of this thesis will detail the methods, beginning with the wind tunnel data sets used for validating the CFD simulations, followed by all the elements that go into producing and successfully simulating the domains created. Lastly, the model evaluation methods will be described. The third chapter compares the simulation output from the three turbulence models to determine which model provides the best representation of conditions measured in and above the 5 model forest stand and cutblock. This chapter will also help determine an appropriate simulation time for the domains. Results of this chapter will be used as an indication of the viability of FLUENT and the turbulence models for the simulations being undertaken in subsequent chapters. The fourth chapter compares and tests different canopy and domain setups. This chapter will introduce and compare the results from the domain that had the forest stand made up of individual porous tree elements along with the layered blockstand canopy used in all other numerical simulations of forest stands. Also included are results from simulations on a domain twice the height of the other domains used. As well, a narrow three dimensional blockstand canopy and a domain consisting entirely of a blockstand canopy, without using standard wind tunnel roughness elements that are used to generate proper conditions ahead of the area under investigation. The results of this chapter will address the three objectives: if the simulations performed by FLUENT are adequate in their reproduction of horizontal velocity and turbulence, if one version of an unmodified K-8 turbulence model is superior to another, and if the use of many porous tree elements together in a stand are as good or better than the porous layered blockstands used in other similar simulations. The fifth chapter will assess the effects of gap size on forest edge stress. Many domains with forest stands having various clearing sizes are simulated. These are done to help determine an equation that could be used by forest managers while considering the size of a cut based upon the endemic (peak annual) winds for an area. The results of this chapter will be used to produce a means of determining optimum cutblock size to minimise wind damage on the windward side of a clearing. Finally, the last chapter will be a summary of conclusions from the results presented in chapters three through five. This chapter will also include reflections of the research completed 6 and provide recommendations for future research. 7 2. Methods This chapter describes many of the programs and tools used to collect data and simulate winds throughout the entire thesis. Some items may be introduced or repeated in subsequent chapters for readability. This chapter is mainly for concepts and materials that will be applicable for all chapters or provide pertinent information that applies to the entire thesis. This chapter starts (section 2.1) by describing the validation data collected in the field and wind tunnel by Dr. Novak, then describes (section 2.2) the CFD model used and finally discusses the statistical methods used to compare model results with observations (section 2.3). 2.1. Validation of data sets This section describes the observed data sets used for validation. Measurements from wind tunnel experiments were provided by Dr. Michael Novak and members of his research group at the University of British Columbia (UBC). Dr. Novak modelled conditions at the Sicamous Creek Research Forest in the UBC wind tunnel (Novak et al., 2001). Although data collected in the research forest was not used in this thesis, it is still important as it was the basis for the wind tunnel work which was determined to be a good model of real world conditions. Wind tunnel data were collected by Dr. Novak and his research team from UBC. A description of that work is included here because the data they collected was used in the validation of the present 8 work. Field observations were collected by Dr. Novak and his group at the Sicamous Creek Silvicultural Systems Research Area in the southern interior of British Columbia, north of Kamloops. Wind tunnel modelling was performed at the UBC Department of Mechanical Engineering located in Vancouver. Dr. Novak and his research team found that wind tunnel measurements were in good agreement with field measurements (Novak et al., 2001). 2.1.1. Sicamous Creek research forest The research conducted at the Sicamous Creek Silvicultural Systems Research Area, provided the data necessary to validate Dr. Novak's wind tunnel model (Novak et al., 2001). Sicamous Creek was the site of many research projects as it was part of a much larger investigation by the BC Ministry of Forests. The main objective of the broader research initiative was to investigate alternative methods of clearcutting (ie. various partial cutting applications) for their biological, operational, social and economic implications (Puttonen and Murphy, 1997). The research area was located in an Englemann Spruce Subalpine Fir (ESSF) biogeoclimatic zone (Vyse, 1997). Each plot in the research area was approximately 30 ha in size; the plot used by Dr. Novak was B-5, which consisted of one clearcut, 10 ha in area (specifically 334 m x 326 m) (figure 2.1) (Vyse, 1997 and Novak et al., 2001). The plot itself was oriented about +10° off true north with a north facing slope of about 12% (Novak et al., 1997). In order to maintain plots that were independent of one another a 100 m buffer strip of mature trees bordered the 10 ha clearcut on all sides (figure 2.2) (Vyse, 1997). The trees surrounding the clearing were estimated to be about 30 m in height with more than 55% (both Engleman spruce and subalpine fir) being between 100 to 150 years old (Novak, 1997 and Parish, 1997). Parish (1997) reported a stand 9 density of 1526 stems per hectare with just under half of those (720) having a diameter at breast height (dbh) greater than 15 cm; the stand was made up of mainly subalpine fir (90%) with some Engelmann spruce (10%). In terms of standing volume of live timber, Vyse (1997) reported that Engelmann spruce represented 35% of the live timber volume compared to subalpine fir which was 65%. 2.1.2. UBC wind tunnel The research conducted using the wind tunnel at UBC provided the data necessary to validate the modelling done using the CFD program FLUENT. The wind tunnel used by Dr. Novak and his research team, was an open-return blow-through wind tunnel operated by the UBC Department of Mechanical Engineering; this wind tunnel is 25 m long (x axis), 2.4 m wide (y axis) and 1.5 m high (z axis) (Chen et al., 1995; Liu et al., 1996; Novak et al., 2000 and Novak et al., 2001). A series of upwind elements were used to recreate the proper flow conditions ahead of the forest section; these upwind elements consisted of Counihan spires, bluff bodies, large and small roughness elements (figure 2.3) Counihan spires (Counihan, 1969) were used to generate large-scale turbulence. Each Counihan spire is roughly equivalent to 1/4 of an ellipse, with the curved edge facing into the direction of airflow. An array of five Counihan spires were used in the wind tunnel (figure 2.4). The array was positioned on top of a 0.03 m thick piece of wood, which gave the spires a total height of 1.23 m. The length and width of the spire decreased as the height increased, at its base, each spire was 0.61 m long and 0.105 m wide; at the peak, the spire width was reduced to 0.013 m, which was the thickness along the outside edge of the spire (refer to figure 2.3 for measurements). A small cylindrical brace connected the spires together; the diameter of the brace was 0.005 m and was located 0.145 m down from the top and 0.095 m 10 ahead of the back edge. Bluff bodies were placed 0.5 m and 1.0 m down stream of the spires to produce intermediate turbulence; they were 0.05 m long, 0.25 m in height and extended across the entire width of the wind tunnel. Large roughness elements were used to generate smallerscale turbulence and were arranged in a diamond shaped pattern, located approximately 1.2 m to 4.8 m downwind of the spires with a density of approximately 21 m~2 (Chen et al., 1995). Each large rectangular element measured 0.15 m in height, 0.05 m wide and 0.025 m thick (figure 2.5 a). Six metres downwind from the spires, small roughness elements followed for 1.2 m with a density of 120 m"2; these were used to help reduced the total length of the forested section (Chen et al., 1995). The small roughness elements were the same height as the larger ones, but were half the width (0.025 m) and only 0.013m thick (figure 2.5 b). 2.1.2.1. UBC model forest The model forest constructed for the wind tunnel was based upon the Engelmann spruce. Model trees were scaled down 200 times to represent trees 30 m in height. They were arranged in a diamond shape pattern giving a density of 500 stems m"2. Each model tree was made from an artificial Christmas tree branch (Barcana, Granby, Quebec) with a height approximately 0.15 m (figure 2.6). Each model tree consisted of a pair of twisted steel wires (0.0009 m) as the main stem, with flexible plastic strips 0.03 m in length and 0.001 m wide, angled at 40°, to act as the foliage (figure 2.7) (Chen et al., 1995 and Novak, 2003a). The bottom 10% of the tree had the foliage removed and the top 30% was trimmed to a point. The diameter of these trees was approximately 0.045 m. The resulting model tree had a leaf area index of 6.3. The model forest began seven metres downwind of the spires, adjacent to the small roughness elements. A clearing measuring 1.63 m by 1.63 m was present 2.4 m into the forest and represented the 10 11 ha plot at Sicamous Creek. 2.1.3. Instrumentation/Data collection Six towers were erected in the Sicamous Creek research forest to observe wind speed and direction. The six towers (1-6) were laid out in an east-west transect through the centre of the B5 plot and placed 25 m, 75 m, 125 m, 200 m, 250 m and 300 m downwind of the edge. Each tower was 8 m in height and capped with R. M. Young anemometers (05103 and 05305, Traverse City, Michigan, USA) (figure 2.8). Instantaneous wind speed and direction were sampled every two seconds by Campbell Scientific 21X data loggers (Logan, Utah, USA); accumulated samples were then averaged and saved at 30 minute intervals. Data collected were only used for analysis if wind speed was within 15° of the transect direction, at all six stations. Three dimensional wind speeds were sampled in the wind tunnel using a Dantec tri-axial fibre-film, hot wire anemometer (55R91) (Dantec Dynamics, Denmark). Samples were taken at ten different locations in the wind tunnel; six locations were within the clearing and four locations within the forest, two each upwind and downwind of the clearing. Thirty points were sampled (vertically) at each location; samples began at a height of 0.02 m and following the second sample at 0.03 m, were recorded at 0.02 m intervals up to 0.59 m (ie. 0.02 m, 0.03 m, 0.05 m, ... 0.59 m). The first two locations within the clearing did not have data from all 30 points. The first site downwind of the clearing was missing data from heights 0.49 m, 0.51 m and 0.53 m. The second clearing location had two values recorded at height 0.31 m, but none at 0.33 m, so the second value was removed and 0.33 was left with no value. In total, over 300 points were sampled in the wind tunnel, by UBC researchers. Those points were sampled for 20 s at 500 Hz, which means each point was sampled 10000 times, this was done in order to provide 12 enough detail to identify turbulent fluctuations. 2.1.4. Validation of wind tunnel Measurements collected from the wind tunnel were deemed to be representative of field conditions. Novak et al. (2001) reports that wind tunnel measurements agreed with measurements taken in the field (figure 2.9). These wind tunnel measurements are therefore used to assess the viability of using FLUENT to predict wind and turbulence within and above forests and forest clearings. 2.2. Computational Fluid Dynamics Discussion of methodology for the computational fluid dynamics portion will be split into three main sections: FLUENT, two dimensional domains and three dimensional domains. The reason for treating two dimensional and three dimensional domains separately was due to the differences that existed between the design and building of these domains. The three dimensional wind tunnel had to be simplified into a two dimensional domain. Computational fluid dynamics (CFD) is a method of numerically modelling fluid flow and heat transfer in simple and complex geometries. CFD was employed to determine if it is capable of adequately predicting flow recorded by Dr. Novak and his research team at UBC in both the wind tunnel, and, by extension, the field. There are four parts to producing a successful simulation: domain production, pre-processing, simulation and post-processing. Domain production entails the creation of a geometry where fluid will flow or interact with objects under 13 investigation. Once nodes (points), edges (lines), faces (surfaces) and volumes (if three dimensional) are produced, a grid/mesh is applied, followed by defining boundary conditions and continua. Pre-processing is done to assign attributes to the domain (including boundary conditions) and to the fluid. Such attributes can include setting fluid velocity or selecting the type of turbulence parameterisation to be used. Simulation calculates all of the governing equations for each grid cell. This process can be time consuming depending on the type of problem you are attempting to solve. Post-processing involves exporting the simulated output for both qualitative and quantitative examination. All aspects of the CFD modelling, were performed on either a 28 processor, SGI Origin 3400 or 64 processor, SGI Altix compute server, in the high performance computing laboratory at the University of Northern British Columbia in Prince George, British Columbia. 2.2.1. FLUENT The computational fluid dynamics software package used for simulations in this thesis is FLUENT, which is released by ANSYS, Inc. (Canonsburg, Pennsylvania, USA). The two main components of this software package are GAMBIT® (hereafter referred to as GAMBIT) and FLUENT. GAMBIT is the unified domain building and meshing program. FLUENT performs the pre-processing, simulation and post-processing. Pre-processing in FLUENT assigns attributes to the domain (including boundary conditions) and to the fluid. Such attributes can include setting fluid velocity or selecting the type of turbulence parameterisation to be included. FLUENT also performs the actual simulation, including compartmentalising of the domain for multiprocessor simulations. Post-processing is the presentation or analysis of output from the simulation. 14 2.2.1.1. Continuity and momentum equations For all flow simulations FLUENT solves mass conservation (continuity) and momentum conservation equations. For flows similar to the ones encountered in this thesis, FLUENT (ANSYS, 2006) uses the following continuity equation: - £ + V • (pu) = 0 (2.1) at where p is the density, t is the time in seconds, u is the vector velocity. The gradient, V , is equivalent to ——. For the conservation of momentum, FLUENT (ANSYS, 2006) uses the follow equation: — (pu) + V-(puu) = -Vp + V • (rj + pg+F (2.2) where/? is static pressure, ng and J? are gravitational and external body forces, respectively. p can also contain additional model-dependant terms (eg. porous media). FLUENT (ANSYS, 2006) notes the stress tensor T and defines it as, T = jU (Vu + VuT)--V-uI (2.3) with JU being molecular viscosity, lis the unit tensor and the term in the parentheses is the effect of volume dilation. 2.2.1.2. Reynolds-averaged Navier-Stokes equations There are two general methods that have been devised to model turbulent flow, which FLUENT supports. These methods do not attempt to directly simulate turbulent flow; they either 15 average or filter turbulence. The Reynolds-averaged Navier-Stokes (RANS) equations employ a set of transport equations that attempt to determine the mean flow quantities for all scales. With this kind of averaging, the exact Navier-Stokes equations are decomposed into mean and fluctuating components for velocity and other scalar quantities like pressure and energy. Large eddy simulations (LES) do not average the Navier-Stokes equations but use a set of filtering equations that remove eddies smaller than a predetermined size, usually the grid size. The LES can be more computationally intensive than the RANS approach, and is only suggested for flows involving simple geometries (ANSYS, 2006; Ferziger and Peric, 1999). Also, LES is only available for use in three dimensional domains. Due to the complexity of the expected flows, a RANS approach was employed. To computationally predict airflow within the domain, FLUENT takes the instantaneous variables from the Navier-Stokes equations and decomposes them into mean (ensemble or time averaged) and fluctuating components. For horizontal velocity, this is: Ui - Ui + U. (2.4) where u , the instantaneous velocity, is equal to the sum of the mean velocity, S"., and the velocity fluctuation, u', from the mean. Other scalar quantities, 0 , like pressure, p, can also be represented: (p=^ + f. (2.5) Using these Reynolds averaged flow variables in the instantaneous continuity and momentum equations become, in Cartesian tensor form: ! + £(/»,) = <> 16 <2.«) and + (/*0 3—(/*<«•";) dt dx. dp d •+ C?X, • <^X, <9k A Kdxi •+ • N de a£er^eff M, dX, ,2 + Cl£-(GK K + (2.14) C,£Gh)-C2ep—-Re fC GK, Gh and YM have the same meaning as the standard K-e model. Inverse Prandtl number numbers, aK and ae, are equal to a, 20 Q (2.20) =Q ij - 2 f ij &> iJK K and Q «,- = " , y - W (2.21) with Q.. being the mean rotation rate as viewed in a rotating reference frame with the angular velocity (coK) and with A0 = 4.04 and As = V 6 COS 0 (2.22) ^ = -cos_1(V6W) (2.23) S:;S j^S^ '^ J * "' , (2.24) when and where W = 5= ft^: and 22 (2.25) —- +—'- u 2 dx: (2.26) d.'x)J The eddy viscosity coefficient still goes into the same eddy viscosity term as in the standard K-S model Mt=pCu K (2.27) — The model transport equations for the realizable K-S model are: (P*)+-£;(P*UJ) * dt (2.28) dK OJ dx- d Mt dX; + GK+Gh+pe-YM+SK for the K equation (the same as the standard K-S model, and for the e equation: ft(pe) + ^-(p£u:) = d JU +Ik. (2.29) de pClSe-pC2e-^-r=+CiekC^Gh + 9J K+^V£ K where C, = max 0.43, 7] r/+5 (2.30) and ?J=S- (2.31) £ Terms carried over from the standard K-S included GK, Gb and YM. Model constants were 1.44 23 for CIs and 1.9 for C2e. Turbulent Prandtl numbers for K and 8 (aK and at) were 1.0 and 1.2, respectively. 2.2.1.4. Method of solution FLUENT uses a control volume and employs an implicit segregated solution method to solve the governing equations. The implicit linearization method solves a given variable in each cell using one equation based on the governing equation for that variable. This is done using both known and unknown variables from surrounding cells. The segregated solution method means that these governing equations are solved in a specific order or segregated into a number of steps. The first step is to update the fluid properties based on the current solution, or if the simulation has just begun, the initialized solution. The second step involves solving the momentum equations which use the current values for pressure and mass fluxes in order to update the velocity field. The third step uses a "Poisson type" equation to perform a pressure correction which is used to update velocity fields and mass fluxes so that continuity in the calculations can be satisfied. The fourth step is where the turbulence models are solved using the updated values from the other variables. The final step, to determine convergence, is then made; if the solution is not yet converged, then the process is repeated again. Convergence is determined by examining and determining whether the residual sum of each conserved variable (continuity, u, v, K, S and w if 3D) is less than a predetermined value. The convergence criteria for residual values for these simulations is 1 x 10"4 (dimensionless property). Residual values for each variable is determined by: 24 Hcel,p\apM Where (J) is the variable of cell P, and b is the contribution of the constant part of the source term. The coefficient of the cell centre is aP and anh is the influence coefficient of neighbouring cells. A second order upwind scheme was used in the discretization of the domain. This is done by using a Taylor series expansion of the cell, where the gradients to adjacent upwind cells is determined so that appropriate values for the cell faces can be determined. Due to early convergence errors in some simulations, the solver was changed from solving the equations as a steady state solution to an unsteady solution. Unsteady solutions were operated at a frequency of 1000 Hz and 2500 Hz (1000 and 2500 timesteps per second, respectively) and run for either 4.0 or 4.5 seconds, which is more than enough timesteps to allow for disturbances caused by the increasing wind speed though the domain to find equilibrium. A successful simulation is determined when residuals (and output values) have levelled off, or are no longer decreasing. Unsteady simulations were run for 4 to 6 seconds of simulation time. Two dimensional unsteady simulations using one processor took, on average, 9.5 hours to complete. Successful three dimensional simulations, using one processor, took 48 to 72 hours. 2.2.1.5. Post processing For post processing, The CFD program used offers many tools to help users view simulation output. Graphically, various isoline plots can be produced along with vector plots. Output values at specific points, along lines and surfaces can be marked and plotted graphically. 25 Values can also be exported for further analysis. Horizontal velocity and turbulent kinetic energy output were marked and exported from simulations in chapters three and four at 286 points. This was done using a script written specifically for these simulations to mark the simulated output at nearly the same points as in the wind tunnel. Ten locations, with roughly 29 points at each location were marked. Also, in the stands made up of tree elements, slight alterations in the position of the four locations within the stands were done to ensure that the output was not from within the tree. The script written for exporting the output was similar to the two dimensional one previously discussed. This script also designated points from 0.03 m to 0.59 m to be collected at intervals of 0.02 m from each location. Output from the marked areas of the domain were exported to text files for analysis. As previously mentioned, two locations from the wind tunnel data set were missing a total of four points, these points were not included in the script. 2.2.2. Two dimensional domains A flow domain similar to the UBC wind tunnel used by Novak et al. (2001) was reproduced using GAMBIT software, a mesh authoring program packaged with FLUENT. The two dimensional domain most used was a total of 12.59 m in length (x) and 1.5 m in height (z), and referred to as the standard domain, based on length and height (figure 2.10). There was no cross-wind dimension. Due to the two dimensional properties, some changes were employed to express the wind tunnel in this manner. Counihan spires, used in the wind tunnel, could not be included in the two dimensional simulations as they would effectively function as a large bluff body and block the flow. This would greatly alter the dynamics of the flow causing compression and a build up of pressure 26 before the spire and greatly increasing the wind speed flowing over the spire, resulting in a simulation that would be no different if the fluid inlet was 30 cm from the top of the tunnel with a wind speed much greater than what was set in the UBC wind tunnel. Two bluff bodies, 0.55 m apart from one another were placed 1.5 m downwind of the domain inlet. Each bluff body was 0.25 m in height and 0.05 m wide. Two arrays of roughness elements followed, both of which were identical in height at 0.15 m, the same height as the model trees (ht). Large roughness elements were wider and spaced more sparsely than the small roughness elements. Large roughness elements began 2.39 m downwind, were 0.025 m thick and were positioned 0.38 m apart from one another. Small roughness elements began 5.8 m downwind of the inlet. They were about half the thickness of the large roughness elements (0.013 m) with subsequent rows 0.1846 m apart. In chapter four, some simulations were completed with a domain 3.0 m in height, twice as high as the standard domain (figure 2.11). This was done to understand the higher peak wind speeds towards the upper portions of the locations where output was collected. Also in chapter four, a domain was created with dimensions similar to the standard domain but lacking the bluff bodies and roughness elements. These structures were replaced by an extended blockstand canopy (figure 2.12). In chapter five, ten additional domains were created similar to the standard domain used. These other domains were used to test various clearing sizes between forest stands and will be described in more detail in that chapter and below. 2.2.2.1. Forest stands Using the standard domain layout, three types of forest stands were constructed to be used 27 in simulations. Of these, two were porous, like vegetation canopies (eg. Li et al., 1990), and the other was comprised of solid elements. One common trait, across all domains, is stand height. The height of the stand in these simulations is 0.15 m (ht), the same height used by the team at UBC (Chen et al., 1995; Liu et al., 1996; Novak et al., 2000 and Novak et al., 2001). The forest stands used in the domains are divided by a clearing; the length of the upwind forest stand is 2.38 m (15.87 ht) and the length of the downwind forest stand is 0.985 m (6.57 ht ). In chapter four, the length of the upwind forest stand for those simulations testing an extended canopy layer, instead of using bluff bodies and roughness elements, is 7.14 m (47.6 ht). All domains, except those in chapter five, have a forest clearing equal to 1.63 m (10.87 ht ). 2.2.2.1.1. Trees Two different homogeneous canopy representations were used to simulate the canopy (figure 2.10). One involves using tree like objects that have a similar dimension to those used in the wind tunnel (figure 2.10 a). The trees used in the two dimensional numerical simulations were slightly different from those in the wind tunnel. Current CFD techniques do not allow for such a fine degree of detail such that individual foliage elements could be replicated from the wind tunnel models. A shape based upon the silhouette of the wind tunnel trees was used as a template to author numerical trees (figure 2.14). Each numerical tree crown had a maximum diameter of 0.04 m and consisted, from the ground up, of a 0.001 m diameter stem that was 0.035 m tall, an inverted frustum 0.02 m in height with an angle of 45° was positioned next. On top of the frustum was a cylinder with a height of 0.05 m. The crown of the tree was a second frustum that was 0.045 m tall. Rows were separated by 0.045 m (eg. a 0.005 m gap between crowns). 28 In order to test the design concept of using tree shaped elements in a forest stand, the first use of these trees were as solid elements within the simulation. To create these, the trees were assembled from the shapes described above. The tree was then simply reproduced using a copy and paste function in Gambit. When the stands were completed, the "subtract" command was used to remove the trees from the domain. In Gambit, the subtract function removes one or more smaller volume(s) from a larger volume. When that occurs the exterior structure of the smaller volume(s) are retained. The structures left behind were merged with the bottom of the domain, effectively making the trees solid objects. To make the trees in the domain more representative of actual trees, which do allow air to move through the canopy; the same steps of building and replicating trees were repeated and the subtract tool in Gambit was used again, but the "retain" function was also used in order to not completely remove the objects from the domain. This allows for the trees to retain both their external and internal attributes without causing any logic conflicts that will be explained in section 2.2.2.3. Once completed, these porous trees will have a special boundary condition applied. 2.2.2.1.2. Blockstands The second type of forest stand used did not have the tree shaped representations, instead a layered rectangular block acted as the canopy. As there are two stands, there are two blocks making up the canopy upwind and downwind of the clearing. These layered blockstands occupy the same length (2.36 m for stand 1 and 0.98 m for stand 2) and height (0.15 m) as the tree stands (figure 2.10 b). The blockstands were divided into five layers, each layer represented a specific leaf area density (LAD) based upon the model trees used in the wind tunnel (further described 29 in 2.2.2.3.1.)- Table 2.1 identifies the height (z) and LAD of each layer. 2.2.2.1.3. Clearing size As previously mentioned, the size of the clearing was 10.87 ht, based on the size of the clearing in the wind tunnel which attempted to reproduce the 10 ha B-5 clearing at the Sicamous Creek research forest. Ten additional domains with different clearing sizes were created for the portion of the study described in chapter five. Each domain had large and small roughness elements, preceding the forest section. The size of all elements (large and small roughness elements and forest stands both upwind and downwind of the clearing) were identical to those used in previously mentioned domains (figure 2.10); the only difference being the size of the clearing. There were 11 clearing sizes used, the different sizes can be found in table 2.2. The stand type used in each of the domains was the porous layered blockstand. 2.2.2.2. Mesh A non-conformal irregular triangular grid scheme was utilised in meshing the domain. A boundary layer mesh was tested and used but only in the domain that featured an extended canopy layer instead of bluff bodies and roughness elements. It was not selected for the other simulations as it was deemed too difficult to implement for those domains which featured individual trees as the forest stand. Edges were first assigned nodes and from that, Gambit would apply a mesh to a face. In previous studies (eg. Green, 1992; Li et al., 1990 and Liu et al., 1996) grid size was smallest within the canopy and 1/10 the size of tree height. All domains in this thesis had distance between mesh nodes less than or equal to 1/10 ht (0.015 m) and used above roughness elements, 30 stands and the centre clearing up to a height of 2 ht. The mesh size at the ceiling of the domain was larger at 0.3 m, given the insignificance of the flow along that boundary. The boundary where air flowed into the domain, had tight node spacing close to the ground and gradually expanded towards the ceiling. This type of mesh grading scheme, called "Last First Ratio", had a ratio value ®) of 1.23 and was determined by: y/(„-i) f R= (2.33) u vv where /„ is the last interval on the edge (0.3 m), /, is the first interval (0.015 m), and n is the number of intervals/nodes along the edge (15). Along the outflow edge, fine node spacing (< 0.015 m) was used up to 2 z/ht (0.3 m) from the ground, then expanded to 0.3 m at the top. A last first ratio of 1.31 was also used along this edge, but only 12 nodes were used in total. After the edges had been meshed, the faces making up the forest stands were meshed along with the wind tunnel itself. By using the existing meshed edges to extrapolate, GAMBIT applies a fine mesh to the domain that is more representative of the 0.015 m mesh applied to the forest stands and ground rather than the coarse mesh found along the ceiling. Most of the mesh in the domain remains fine even as it approaches the ceiling (figure 2.14). 2.2.2.2.1. Grid Refinement The advantage of using an unstructured mesh, like the one used for this thesis, is the ability to save time setting up the grid compared to structured grids, and being able to incorporate "solution-adaptive" refinement to the grid. This tool aids with grid independence by expediting the refinement process. The CFD program used includes this refinement feature that analyses 31 the solution to determine where more cells can be added to improve the simulation. This allows the grid to be refined specifically where it is needed without the need to completely generate a new grid. Two types of grid adaption were offered by the CFD program. The first type, known as "hanging node" adaption, is a type of refinement that isotropically subdivides one marked cell into four, therefore cells adjacent to the refined cell will have their edge split in two to accommodate the two new cells (figure 2.15). A triangular cell, the cell types of this grid, would be subdivided into four triangles (figure 2.16). The other type of grid refinement offered is known as "conformal" adaption. With this type of adaption, edges of a marked cell are examined to determine the longest edge to split (figure 2.17). Neighbouring cells that shared that edge are also searched to see if any of that cells edges are longer and should be split. This continues until the there are no longer any edges in the adjacent cells that are longer than the last edge split. Hanging node adaption was ultimately used. Also, refinement was limited to one level of refinement between neighbouring cells. Gradient adaption is a method used to identify areas for grid adaption and was used for that purpose. Gradient adaption uses the gradient of a selected field variable (ie. x velocity) to determine which cells would benefit from refinement. In the CFD program used, gradient adaption displays the maximum difference of a variable between two adjacent cells. Generally, cells that have a difference with their neighbours that is > 10% of the maximum difference for that variable, are marked for refinement. This CFD program also offers a registry that allows the user to mark the cells of many different variables. This organises all of the cells to make sure that they are not repeatedly being marked and refined. Nine variables were examined for refinement after a converged solution with the K-e 32 turbulence model . These variables, static pressure, dynamic pressure, velocity magnitude, xvelocity, y-velocity, turbulent kinetic energy, turbulence intensity, turbulence dissipation and turbulence viscosity, were marked using the 10% threshold previously described. 2.2.2.2.2. Grid independence A second finer mesh was applied to the standard domain containing the forest stand with trees. Distance between nodes was made 1/3 smaller than the regular mesh, yielding a mesh nine times finer (904068 cells). The results of grid sensitivity between the fine and standard mesh are plotted as a vertical profile (described in more detail in section 2.3.1.3) in figure 2.18 along with the wind tunnel measurements at ten locations within the forest stands and clearing, from the ground up to 4 z/ht. 2.2.2.3. Boundary and Continuum Types In CFD, a domain is made up of boundaries and continua; boundaries interact with the flow and continua are the areas (2D) where flow properties are assigned and occur. To give a simplified example, in a numerical wind tunnel, there is a large rectangular object that represents the wind tunnel or domain. The rectangle has edges bounding the space or continuum. Specific attributes of that continuum are applied over the rectangular areas. If another object were to be added to the rectangle, there would be a conflict as the new object brings another set of attributes to the area/continuum it occupies. Two continua cannot occupy the same space. The continua have to be separated by boundaries. When objects are required to occupy the same space, there are a couple of different techniques that can be used to satisfy the logic requirements of the program. 33 The first technique involves removing one of the continua so that the other can remain. The "subtraction" function in Gambit is used to subtract an area where two or more objects overlap. Normally, under default settings an empty space remains that is bounded but contains no continuum, which is fine for objects that are solid. This is how the trees in the solid forest stand were made. For the other case involving the porous trees in the forest stand, the tree shaped object is subtracted from the wind tunnel, but the "retain" function is selected allowing the normally empty space to retain the tree shaped object (figure 2.19 a). In this case, two separate faces are made, the wind tunnel face and the tree face, each with their own continuum and flow properties, and separated by a two-sided "internal" boundary condition. This same method was used for the blockstand configuration, with each layer representing a distinctive continuum. In order to mimic the conditions of the wind tunnel, the simplest options were used when assigning boundary conditions to the domains. The inlet was designated with a very simple velocity inlet, which was set at 8.7 m s"1 to correspond to the velocity used in the Novak et al. (2001) wind tunnel. The outlet boundary was set as a pressure outlet boundary condition, which is used to improve results in instances where reverse flow occurs during iteration. The roof and floor (including bluff bodies and roughness elements), just like the wind tunnel, were solid and interacted with the flow accordingly. The domain continuum, where the flow occurs, was designated as air, with a density of 1.225 kgmf3 and a dynamic viscosity of 1.7894 x 10"5 kgnr's"1. Temperature, gravity and buoyancy were not factored into the simulations. 34 2.2.2.3.1. Canopy properties Aside from the forest stands made up of solid trees, special continua had to be assigned for the porous tree and block forest stands. There are conceivably four different properties that could be given to the canopy in the numerical wind tunnel based on what was previously outlined. The first one (and most useless) would be if the canopy of the trees was just given the same properties as the continuum which surrounds the trees, in essence it would act as if there was no canopy at all. The second, also a simplified interpretation, would be to have the canopy act as a solid and therefore no airflow through the canopy. The third and fourth possibilities deal with a canopy that would have porous properties. As previously stated, the authoring of tree limbs, branches and foliage is currently too complex to be produced in a computer simulation. The CFD package used, has two ways to deal with diffuse flow through an object; two modules called "Porous Jump" and "Porous Media", which allow for the movement of fluid through a barrier. The first module, porous jump, is a ID membrane that has known velocity/pressure drop characteristics, like a screen (ANSYS, 2006). This module is applied to an objects cell face or a designated interface between grid cells. These velocity/pressure drop characteristics are uniformly applied over the entire object face; therefore, the varying thickness of an object would not cause any additional decrease in velocity. It does not appear that this would fairly portray the movement of air through a model of a cone shaped conifer tree as air flowing through the upper sections of a real tree would not be impeded as much as the air flowing through lower sections due to differences in crown diameter (thickness of the barrier) at any given height. A greater thickness of foliage toward the bottom of a crown would be more capable of slowing wind, than the upper section of the crown where the thickness of the foliage is much less. The second module, porous 35 media, seems more appropriate, as the module applies a determined flow resistance to a region (cell zone/volume/thickness) that is defined as "porous" (ANSYS, 2006). This module seems like it would be better able to take into account the properties of a cone shaped structure such as a tree. However, Fluent (2001) recommends the use of the porous jump utility whenever possible because it yields better convergence of the solution. The forest stand with the tree shaped elements had a continuum assigned to the interior of the crown. The porous media continuum (fluid property) was assigned to tree stands (upwind and downwind stands) to handle the diffuse flow through the trees. The porous media module is essentially a momentum sink term made up of two components, a viscous loss term and an inertial loss term. The momentum sink equation used in The CFD program is, u S- = - ; + C0 KLAD 2 —du\u. ' ' ' (2.34) where /u is the viscosity of the air, p is the air density, u-t is the velocity and \u\ is the absolute value of the spatially and time averaged wind speed and C2=2CdLAD (2.35) where Cd is the drag coefficient. As momentum sink terms for canopies in the meteorological literature are not multiplied by one-half (eg. Amiro, 1990; Raupach and Thorn, 1981), the C2 parameter in the momentum equation needed to be multiplied by two. Leaf area density {LAD) which is the area of foliage contained within a given volume of canopy is expressed as m2m"3 or m"1. An error in the interpretation of the porous media term was only caught after all the simulations were analysed. In the viscous loss portion of the porous media module, the FLUENT manual had the denominator as alpha (a) which in the forestry literature commonly represents 36 leaf area density (m"1) but the manual had it representing permeability (m2). Confounding the problem further was that this mistake was not caught while balancing the equation, so no difference was found between the viscous and inertial loss terms. This error was not cause for concern as the permeability value needed to be many orders of magnitude smaller to influence the sum of the equation. This is because the viscous loss term, JU (dividend), was so small (1.7894 x 10"5) that it would also take an incredibly small permeability value (divisor) to affect the quotient. Another error that may have had a larger impact in this model was the way C2 was calculated in equation 2.35. The drag coefficient used in this term was 1 and based on what was reported in Novak et al. (2001); however, this value was calculated using the frontal area of the model tree used (Af = 0.0043 m2). This maybe incorrect for the current use which should be based on a volumetric drag coefficient in which case the surface area of the model tree (As = 0.011 m2), and would result in the drag coefficient being equal to 0.4. A smaller drag coefficient would ultimately result in less momentum sink by the canopy. The leaf area density for the porous tree forest strands was 55.5 m2 ~ 3 based upon values given in Liu et al. (1996), which used the same artificial tree foliage, for different levels of the canopy. Other authors (eg. Wilson, 1988; Green,1992; Liu et al., 1996 and Foudhil et al., 2005) have also included a term in canopy properties to account for vegetation wake turbulence. This turbulence is generated at fine length scales and are characteristic of canopy elements (Foudhil et al., 2005). This term could not be applied to the porous media module. 37 2.2.2.3.2. Grid Interface Meshing problems were experienced in some of the early domains. The inability to successfully mesh an entire domain in GAMBIT led to the use of the "grid interface" tool. This tool in FLUENT is used to unite the grids of two adjacent domains so that fluid can freely flow from one domain to another. The CFD program analyses the grid on both domains and will refine the grid along each edge so that grid point nodes match-up along the interface boundary (figure 2.20). 2.2.3. Three dimensional domains Many aspects of flow setup for three dimensions is similar to those in two dimensions. This section will describe those portions of the thesis where three dimensional domains were created and meshed in preparation for simulation. Three dimensional domains, like the two dimensional ones, were authored using Gambit. The size and complexity of the simulation that was being attempted made it impossible for the entire wind tunnel to be wholly created together. Therefore, some changes were made in order to successfully simulate airflow within the three dimensional domain. As with the two dimensional domain, locations in the domain were referenced from the centre of the most upwind tree (x axis), the centre width of the domain (y axis) and the floor (z axis). Due to the size and resolution, a full scale three dimensional wind tunnel could not be produced. After consideration, a 0.088 m wide domain was produced, with the same length and height measurements as the two dimensional domain (12.59m long and 1.5 m high). This width was decided upon, as it maintained the layout of many of the upwind elements seen in the wind tunnel; it also allowed a row with two complete trees and open space beside a tree and the 38 domain edge. After a lengthy construction and testing trial, it was determined that actual tree representations with porous canopies could not be successfully simulated in three dimensions. This was because the solution either diverged or never converged and accompanying output was unrealistic (ie. horizontal velocity was >100 ms"1). Instead a three dimensional layered blockstand canopy with similar dimensions to the two dimensional domain was created. In order to successfully mesh the volume, the domain was actually composed of four distinct volumes that were joined together, these volumes included: bluff body - large roughness elements (BB-LRE), small roughness elements (SRE), the first stand of trees and about two thirds of the clearing (Stand 1) and the second stand (Stand 2). It was later learnt that this division may not have been necessary, but it did help by identifying where meshing problems were occurring. Counihan spires could not be included as the base of one spire was greater than the width of the domain. Regardless, their presence in such a small width, would have caused large, undesirable disruptions in the flow both upwind and downwind of the spire, similar to the two dimensional scenario, though perhaps not quite as severe. The first section, which contained the bluff bodies and the large roughness elements was 5.625 m in length and located from -7 m to -1.375 m (x) (with x = 0 at the most upwind location of the forest canopy). Bluff bodies were 0.5 m apart from each other with the most upwind body at -5.5 m (x) and had identical dimensions as the UBC wind tunnel bluff bodies. The roughness elements were also identical to those in the UBC wind tunnel, they began at -4.61 m (x) and were arranged in a staggered diamond shaped pattern with a density equivalent to 21 elements m"2. The second section contained the small roughness elements. It was positioned from -1.375 m to -0.0497 m (x). The downwind boundary was positioned halfway between the last small 39 roughness elements and the first row of trees. The first elements began at -1.2 m (x) and retained the same dimensions as the UBC wind tunnel counterparts with a density equivalent to 120 elements m"2. The third section was the most complicated section. It contained the first stand of trees and 64% of the clearing. This section ranged from -0.0497 m to +3.4 m (x) with the stand 2.38 m in length (-0.02 m to 2.36 m (x)). The fourth section was 2.19 m in length and contained the second stand, bordered by two clearings. The first clearing was part of the main clearing under investigation while the second one, as previously mentioned, was to satisfy flow requirements along the outlet. This section was positioned from 3.4 m to 5.59 m (x) and contained the second stand which was nearly one metre in length (0.985 m) and 1.63 m downwind of the first stand. The second stand had the same characteristics as the first stand. 2.2.3.1. Forest stands The stands in the three dimensional domain were blockstands and identical to the two dimensional stands described in section 2.2.2.1.2. The only difference between the three and two dimensional blockstands was the extra dimension, which was the same width as the domain, 0.088 m. 2.2.3.2. Mesh A non-conformal irregular tetrahedral grid scheme was used in meshing the three dimensional domain. As with the two dimensional mesh, a boundary layer mesh was not used. Edges were first meshed by designating the distance between grid nodes, then faces were meshed 40 using an irregular triangular scheme by using the edges as a template. Once all surfaces were meshed, the surrounding volume was meshed using a four-node tetrahedral approach. Within the stand, distance between mesh nodes was set to no more than 1/10 ht (0.015 m). Boundaries between the four sections were meshed at this value up to 2 ht, after which distance between mesh nodes expanded to 2 ht along the ceiling of the domain. This was done by selecting both vertical edges on each face. The domain boundary where air flowed into the domain, had tight node spacing close to the ground and expanded as it went towards the ceiling. A last first ratio (refer to eq. 2.33 in 2.2.2.2) of 1.23 was used in Gambit along with 15 nodes along that edge so that node spacing was 0.015 m at the bottom and approximately 0.3 m at the top. At the boundaries between sections and at the outlet boundary, a last first ratio of 1.31 was also used along these edges, but only 12 nodes were used on each vertical edge. After all the edges had been meshed the faces making up the stands and the layers within the stands were meshed followed by all the domain boundaries. Once the faces were meshed, volumes were meshed. To help simplify the process, all of the layers in the first stand were meshed followed by the second stand. Then each of the four wind tunnel sections were meshed. Meshing volumes is the most computationally intensive portion of the authoring process. This is where a fair amount of difficulty arose, Gambit did have problems meshing the domain before it was split into four sections. 2.2.3.3. Boundary and continuum types In CFD, a domain is made up of two types of zones which help to define the characteristics of the domain. Boundary conditions are used to define how the simulation will interact with a physical or operational domain boundary. Continuum types define the physical or operational 41 characteristics within the domain. In three dimensions, boundaries are defined by faces and continua are the enclosed volumes where the flow occurs, unless they are designated solid. Boundary conditions were assigned to the different faces of the wind tunnel. These are properties/parameters that were given to the domain boundaries so that The CFD program will know how the fluid (air) will interact with these surfaces. The boundary conditions used in the wind tunnel setup were: velocity inlet, pressure outlet, wall, symmetry and interface. Velocity inlet was used to define the velocity of the fluid at the inlet. This was the equivalent velocity inlet of the UBC wind tunnel; for simulations, velocity was set at 8.7 m s"1. Pressure outlet boundary condition is known to offer better convergence when reverse flows along the outlet have been problematic. The wall boundary condition has the properties of a solid surface. The flow along this boundary will interact so that flow neither goes into or out of the wall boundary. The ground along with the bluff bodies and roughness elements were assigned the property of a wall given that these are all elements where the air is to flow over and/or around them. The roof of the domain was also assigned as a wall. The three dimensional setup consisted of a lateral boundary condition called "symmetry", this function "mirrors" the physical contents of the wind tunnel. Due to the complexity of the numerical setup of the three dimensional wind tunnel, a very narrow strip (0.088 m) of the wind tunnel was modelled. The sides of the chamber were assigned as symmetrical boundaries. Assigning the sides of the chamber as symmetrical boundaries is more consistent with the physical wind tunnel and field conditions, where the forest would continue on after the domain walls. The symmetrical boundary layer condition allows for a hypothetical continuation (mirror) of conditions after the defined parameters of the chamber/study area (ANSYS, 2006). Also the 42 flow does not interact with this boundary like a wall, it interacts with it as if there were no boundary. The symmetrical boundary conditions ensures that there is a zero flux between both sides of the boundary, so what would flow out of the domain would also flow back into the domain. The "interface" boundary condition was used when one single domain could not be successfully meshed. Two interface boundaries or zones are a way of joining volumes together that are perpendicular to the flow. The CFD program then aligns the mesh of the two surfaces together so that the transition from one volume to another is seamless (figure 2.20). 2.3. Model evaluation methods Most analysis completed for this thesis was to validate model output with wind tunnel observations. Therefore this section of the chapter will be split between the model evaluation used for chapters three and four and the evaluation for chapter five. 2.3.1. Output collection for model validation A variety of field variables were extracted from all simulations in chapters three and four. Output from simulations was extracted at nearly the same locations as in the wind tunnel (figure 2.21) and were noted as locations A to J, with location A the most upwind location in the first stand and location J the most downwind location in the second stand. Four locations, all within stands (locations A, B and I, J), were positioned slightly different from the wind tunnel so that sampling points were not within trees for domains that used trees (as opposed to blockstands). 43 This resulted in only a slight alteration from the actual positions from the UBC wind tunnel (<0.01 m). Each location, with the exception of C and D, had 29 sampling points ranging from 0.03 m to 0.59 m in height at 0.02 m intervals. Due to an absence of data from the UBC wind tunnel at sample heights of 0.49 m, 0.53 m and 0.57 m for location C and 0.33 m for location D, the total number of points sampled from all locations was 286. Creating all the sampling points for one domain was very time consuming but after the first one was completed, a script was produced which expedited this process for subsequent domains. Once all of the points were marked, the modelled values were exported as text files. Turbulent kinetic energy was not directly collected in the wind tunnel. TKE values for the wind tunnel were calculated using wind statistics that were available. Using the variance of each wind component collected, TKE can be determined by: 1 /—;2 TKE = -\u' —;2 + v' ;2\ +w' (2.36) Five methods of analysis were used to interpret the output from the simulations and compare them with one another and with the wind tunnel results. 2.3.1.1. Filled isolines Descriptions of output showing horizontal velocity and turbulent kinetic energy isolines were provided to examine how well each value compared to one another based upon the length of the simulation (number of iterations) and differences in how the models were simulated (ie. type of K-S model used and whether simulation was steady or unsteady). In many upcoming figures (also figure 2.21 in which they are illustrated as dotted lines), the grey vertical lines denote the ten locations used to extract points from the simulations to compare with very similar 44 points in the wind tunnel. 2.3.1.2. Scatterplots Horizontal velocity output was extracted from all the K-e turbulence models mentioned above. From ten locations, 286 points were compared with the wind tunnel and simulations. Scatterplots of model output versus observations were produced and showed how well the simulations matched the wind tunnel results. 2.3.1.3. Vertical profiles Vertical profiles of horizontal velocity for each of the ten locations were produced. The wind tunnel profiles were those recorded in the wind tunnel while the profiles from the simulations were extracted as line output; output extracted along the line was from the node centres through which the line intersected. Plotted lines will show how well the models predicted wind speed at a specific vertical location. Due to the canopy properties used in the simulations in chapter 3, the simulated profiles and observed data were expected to differ within the forest canopy. 2.3.1.4. Quantitative analysis Quantitative data analysis used the same 286 points from both the wind tunnel and simulations to produce a variety of statistical outputs based on horizontal wind velocity and turbulent kinetic energy. While many modelling studies have qualitative information from plots and graphs showing how well their models compare with observed data, only a handful show quantitative evidence. As a large component of my thesis entails model evaluation, it was 45 important to have a robust and quantitative way of analysing output from simulations. There are two main measures used in the evaluation of model performance: correlation and difference measures. Most of the work presented in this thesis will be on evaluating models using difference measures rather than correlation measures as these have been criticised by researchers in this field (eg. Willmott, 1981 and 1982). Pearson's product-moment correlation coefficient has been widely used in the past to provide a quantitative index of association as it describes the colinearity between predicted and observed values (Willmott, 1981 and 1982). Willmott (1982) explains that the main problem with r and r2 in published reports is that they are not consistently or properly related to the accuracy of the model prediction. While r and r2 do illustrate the proportional increases and decreases between model output and observations, they do not show the differences between the type and/or magnitude of potential covariation (Willmott, 1981). In a study by Willmott and Wicks (1980), statistically significant values of r and r2 were found to be misleading, as they were unrelated to the sizes of the differences between observed and predicted values. Willmott (1981) explained this by stating that r and r2 do not inform the modeller enough about the actual size of the error that could be produced by a model as a high standard deviation in either (or both) the observed or predicted values could mask a large error. Willmott (1982) finally states that r and r2 should not be used because the relationships between them and model performance are not usually well-defined. Despite this, r2 values will still be reported as they are still integral in the production of Taylor diagrams (section 2.3.1.5). Difference measures, also known as error indices, are favoured in the atmospheric sciences for model performance as they better summarise the mean difference between the observed and predicted values. Root mean square error (RMSE) and mean absolute error (MAE) have been 46 supported as the best method of reporting difference measures (Fox, 1981 and Willmott, 1982). Mean absolute error and RMSE are calculated from (Fox, 1981): N y MAE = N- Y\Pi-°>\ ( 2 - 37 > i=i -,0.5 RMSE w-Z(fl-o.) i=\ : (2.38) where N is the number points sampled and P, is the predicted value from the simulation and 0, is the observed value from the wind tunnel. Fox (1981) noted that MAE was simpler which made it more appealing but that it was not as sensitive to extreme values as RMSE. That fact led Willmott (1982) to conclude that while MAE or RMSE could either be reported by modellers, RMSE may have, in some in-depth statistical analysis, a slight advantage over MAE. While RMSE is a good overall measure of model performance, the modeller still needs to know what types/sources of error are present in order to help refine a model. Therefore, it can be beneficial for the modeller to know how much of the RMSE is systematic and how much is unsystematic. When describing the components of the RMSE, it is easier to remove the rooted component of RMSE to get mean square error (MSE). Systematic MSE (MSES) is composed of three measures: additive, proportional and interdependence, which together give an estimation of a models linear bias. This bias is, in part, a result of the constant over and under prediction from a model (Willmott, 1981). Systematic mean square error is calculated from: N MSES = N-^^-Oi) i=i when 47 (2.39) Pi = b + mOi (2.40) and where b is the intercept and m is the slope from a least-squares regression. The unsystematic error (MSEU) is equally important as it measures how much of the difference between the predictions and observations are due to random phenomena and/or other influences that occur outside of the legitimate range of the model (Doty et al., 2002). It is calculated from: MSE^N^tiP-^f. (2.41) The MSE system is conservative so, MSE = MSEs + MSEu. (2.42) Willmott (1982) does suggest, for reporting purposes, that RMSE be used as it can report in the same units as the predicted and observed values. In addition to RMSE, it was also recommended that the systematic and unsystematic components be reported: RMSES = ^MSES (2.43) RMSEU = ^MSEU. (2.44) and If the RMSE components are not reported, then it may be interpreted that the RMSE value is entirely unsystematic, which could be misleading. If RMSES is high (approaches RMSE) and is greater than RMSEU, then the model could be acceptable, but further refinement should be performed. If RMSEU is close to RMSE then the model may be as good as it can get and no further refinement will be possible. In other words, a successful model should have a ratio of RMSEU to RMSES greater than one. 48 Based upon their mathematical layout, MAE and RMSE provide the modeller and reader with an average error. Recognising the need for a simple value that would allow a reader to easily decipher how well the model predicted, Willmott (1981) and Willmott and Wicks (1980) developed an alternative measure that illustrates how well a measured value was predicted by a model. They described this measure as an "index of agreement" (d), hereafter referred to as Willmott d, that measures the degree to which a model's predictions are error free: d=l- Zte-4)2/Z(M+M) , 0d the intercept (b) and In addition it is recommended that "descriptive 49 measures" be reported (MAE, RMSE, RMSES, RMSEU, and Willmott d) and their interpretation be descriptive and based on scientific grounds and not on the basis of the measures' statistical significance (Willmott, 1982). Also, for reasons stated previously, r2 values will also be reported. 2.3.1.5. Taylor diagrams Finally, summarising model performance of multiple simulations in a single diagram is done using a Taylor diagram (Taylor, 2000). Taylor diagrams use normalised standard deviation and correlation to graphically show how well models compare to one another in relation to observations or reference data. All simulations are plotted on one Taylor diagram to show how each model compares with each other and to the wind tunnel observations. 2.3.2. Evaluation of gap size on forest edges An analysis is performed on each of the domains to determine the wind stress on the forest edge. This was to determine if there is an optimal clearing size in cut block design to minimize the wind stress. Output was extracted from each simulation at 0.23 x/ht ahead of the downwind edge, the same position used by Novak et al. (2001) in a study examining drag on a forest edge downwind of a clearing. The output extracted from the simulations was then fed into the following equation (Landsberg and Jarvis, 1973; Thorn, 1971): F = PZ(U? • Cdiui}-LAD-Azi) (2.48) z=0 Where F is the momentum absorbed (stress) (Pa or kg m"1 s"2), p is the air density (kg m"3), Cd is the drag coefficient (dimensionless) of the trees at speed ut, LAD is the leaf area density (m2 m~3), 50 u is the horizontal wind velocity squared, and Azt is the change in height. Constants from equation 2.48 were used in each of the domains. Air density is 1.225 kg m"\ and the leaf area density is variable depending on the height at which the output is extracted (Table 2.1). The drag coefficient used in this equation was 1, as determined by Novak et al. (2001) using a twostage strain-gauge balance to measure drag on a model tree used in wind tunnel tests; however, as explained in section 2.2.2.3.1., the drag coefficient used in Novak et al. (2001) was calculated differently. For these simulations the drag coefficient should have been 0.4; therefore, stress values determined using equation 2.48 would be about half the reported value. The drag coefficient on the CFD model forest stand based on speed was not determined, but it is documented that trees in nature can decrease their drag in response to increasing wind speeds (Rudnicki et al., 2004). Horizontal wind speed is extracted from the simulation output and AZJ is based upon how the grid points were arranged between 0.03 m to 0.15 m. This is done using the rake/line tool in The CFD program. Points below 0.03 m were not included as this would have represented the stem area and it was also a concern that the floor may also include some unwanted influence on the total value. Between 0.03 m to 0.15 m, a range of 15 to 19 points were typically extracted for each domain. 51 T"' Figure 2.1: Sicamous Creek research forest was established by the British Columbia Ministry of Forests and was the site of many studies into forest health. Observations collected at the 10 ha plot, B-5, were used to validate wind tunnel measurements collected by the UBC team (Image from Puttonen and Murphy, 1997. Used with permission, British Columbia Ministry of Forests and Range, Research Branch). Figure 2.2: Aerial view of B-5 plot. The Engelmann spruce and subalpine fir trees surrounding the 10 ha clearcut are approximately 30 m in height (Photo from Novak etal., 2001). 52 0.088 m 0.50 m 19ir0 D 0.24 m &r 0.15 m Bluff Body Section (1.2 m) 0.50 m ^ m Scale is 1:20 Large Roughness Elements Section (3.6 m) Upwind section of Novak et al. wind tunnel D D 0.088 mlt B D D D D Section (1.2 m) Small Roughness Elements 0.092 m 53 Figure 2.3: A detailed schematic of the upwind section of the wind tunnel used by Dr. Novak and his research team at UBC. Spires, bluff bodies, large roughness elements and small roughness elements are shown. Upwind Spire Section 'm 0.100 r O.'lSJiWJ 0.067 m M 1 0.105 m ; 0.4175 m 1 0.1525 m 0.1525 m 0.61m Trees Trees I •• t r "'. Figure 2.4: Array of five Counihan spires used for generating large scale turbulence in the wind tunnel used by Dr. Novak and his team at UBC (Photo from Novak etal., 2001). B Large Roughness Element Small Roughness Element A- _4=n E E 10 i ,^0.025 m 0.05 m M ^0.013 ^.om 0.025 m Figure 2.5: Dimensions of the large and small roughness elements used in the wind tunnel. 54 0.15 m Figure 2.7: Artificial Christmas tree branches were used to make up the model forest in the wind tunnel (Photo from Novak et al., 2001). 0.045 m| Figure 2.6: Approximate dimensions of each tree in the wind tunnel. Sicamous Bb 10 ru Cleattity •- Mean Wind Spfrec l.4r 1? to c o | 9' CD 0 8 o f 0,4 0.2 0, Figure 2.8: RM Young Anemometer in plot B-5 at the Sicamous Creek research forest (Photo from Novak etal., 2001). b *'• O c Sicamous West (199F. Sicamous t a s t (1996) W l (March 1 & 4, 1996) • WT (March 5 & 6. 1E«6} 10 12 Figure 2.9: Measurements collected from both Sicamous field site and UBC wind tunnel (Image from Novak etal., 2001). 55 0\ 0.55 m 0.55 m 0.38 m 0.38 m I Scale is 1:55 12.59 m - 0.15 m 0.184 m B ^ l I I M I If -12.59 m- 0.15 m 0.184 m n^i i • i • • i 2.38 m 2.38 m 1.63 m 1.63 m 0.985 m 0.985 m _a Figure 2.10: The two main types of numerical domains used for the simulations. A) Multiple tree shaped entities are used to make up the forested areas. B) Layered rectangular blocks represent the forest stands (blockstands). o m o Ji i -- 0.55 m 0.38 m 12.59 m 0.15 m 0.184 m n V i i i i~l i ii 2.38 m 1.63 m 0.985 m 2.38 m I2.59m 1.63 m 07985 m Figure 2.12: The "stand only" domain with a porous blockstand canopy and no upwind roughness elements was used to examine a domain with no upwind bluff bodies or roughness elements. o Scale is 1:55 Figure 2.11: Three metre tall domain used to help explain higher wind speeds seen in the 1.5 m tall domain. This domain also used the porous blockstand canopy. CO o o I 0.04 m 0.045 m E 0.05 m 0.02 m 0.035 m Figure 2.13: The two dimensional numerical trees were based upon the silhouetted shape of the wind tunnel trees and had similar dimensions. Table 2.1: Height ranges and accompanying leaf area densities for the different layers within the porous blockstand canopy, as reported by Liu et al. (1996) in their wind tunnel stand which was the same to the UBC validation data. Layer Height (m) Leaf Area Density (m2m"3) 1 0 - 0.02 6 2 0.021 - 0.05 57.5 3 0.051 -0.08 55.5 4 0.081-0.11 53.5 5 0.111-0.14 38 6 0.141 -0.15 28 58 Table 2.2: Clearing sizes used to determine stress at the downwind edge of a clearing between two porous blockstand forest stands. Clearing Size (m) Size (ht) 0.60 4 1.20 8 1.50 10 1.63 10.9 2.10 14 3.00 20 4.50 30 7.50 50 9.00 60 12.00 80 15.00 100 59 A B Figure 2.14: Portion of the domain showing detail of the mesh along the upwind stand. Grey lines extending vertically are locations where sampling took place for validation with wind tunnel measurements. A) Mesh shown for the standard sized mesh used (adapted). Each grid segment of the triangular grid <2 z/ht is no larger than 1/10 h, or 0.015 m in length. B)To help determine grid independence a much finer mesh was produced, also a triangular grid, but with no grid segment larger than 1/30 h, or 0.005 m in length, producing a mesh nine times finer. 60 > 11 i Hanging Node v i > II i Figure 2.15: Hanging node grid adaption used by (FLUENT i I leaves a grid node unused on one side (after ANSYS, 2006). Figure 2.16: Before (left) and after (right) illustrations of hanging node grid adaption applied to a triangular grid cell (after ANSYS, 2006). • s x B Figure 2.17: This figure shows an example of conformal grid adaption used in FLUENT. Cell A is marked for refinement along it's longest edge. Because Cell B also shares the edge, it too will be split (after ANSYS, 2006). 61 . | U (m/s) U (mis) U (mis) U (mis) Figure 2.18: Vertical profiles of horizontal velocity from both the standard grid ( ) and fine (Gl) grid ( ) for determining grid independence. These profiles are shown with measurements from the wind tunnel (•) for comparison and will be described in more detail in subsequent chapters. The vertical profiles shown are from locations within the domain, with location "a" being the most upwind location sampled and location "j" the most downwind location sampled (see figure 2.21 to see where the specific positions of the sampling locations are within the wind tunnel/domain. U (m/s) B Al Wind tunnel Continuum Wind tunnel Continuum A Wind tunnel Boundary A Tree continuum Two sided Boundary Tree continuum Tree boundary I Figure 2.19: A) Two continua cannot occupy the same space, so a "subtraction" and "retain" function must be used so that the two separate entities can exist. Note: white area only shows separation between two boundary conditions and is not actually present. B) By using a two sided boundary condition, identification and organisation of special boundary conditions or continua becomes easier and more manageable. cell zone 1 interface zone 1 interface zone 2 E \ D F> A IV cell zone 2 A VI Figure 2.20: An example of an offset grid interface between two zones. Interface zone 1 is made up of two faces, A-B and B-C, with interface zone 2 made up of faces D-E and E-F. The intersection of these faces produce new faces: a-d, d-b, b-e and e-c . Face a-d will form a wall zone, while faces d-b, b-e and e-c will form two-sided interior faces. To calculate the flow into cell IV, both faces d-b and b-e will be used to bring in the information from cell I and III, respectively (after ANSYS, 2006). 63 4^ ON B -0.50 D C 0.375 E 0.13 0.625 F 1.005 H G 1.505 I 1.255 . 1.76 J 2.13 Figure 2.21: Positions of the ten sampling locations within the stands and clearing of both the wind tunnel and numerical domains. Moving left to right, Location A is the most upwind location, Locations B through I denoting the next eight locations with Location J being the most downwind location. Numbers indicate the Location's position (in metres) relative to the upwind edge of the clearing. A -1.36 3. Comparison and testing of turbulence models 3.1. Introduction This chapter introduces and tests three variations of the K-S turbulence model on a domain that includes a forest of solid tree shaped elements. This is to ensure that the CFD program would be able to perform adequately prior to subsequent chapters in which the airflow will be interacting with the forest canopy. Airflow in forests, particularly when it is strong enough to cause windthrow, is an important consideration in the harvesting and management of forests. Harvesting of trees can greatly impact the ability of trees adjacent to those harvested to withstand strong wind for the first few years after the cut (Stathers et al., 1994), as the removal of the neighbouring trees alters their exposure to the wind. In order to better understand airflow along forest edges, many people have attempted to model airflow over these forest clearings. In the past, numerical research of airflow through forests, for example Li et al. (1990) and Liu et al. (1996), has been done with models using "in-house" codes, developed especially for simulating wind through tree canopies, usually employing a K-8 turbulence closure scheme. Most numerical airflow simulations utilise a turbulence model as knowledge of turbulent flow is useful in forestry and other environmental fields. Specifically in forestry, windthrow or blow down can occur as a result of dynamic wind 65 loading on trees, usually related to turbulence intensity and frequency (Liu et al., 1996). The turbulence models selected for this study are similar to those used in previous studies involving forested environments (Green, 1992 and Liu et al., 1996) or domains resembling such environments (eg. cylinders, roughness elements) (Belcher et al, 2003 and Lakehal, 1999). This study differs from most others as a commercially available CFD program was used. Another study to use a commercially available program, was Green (1992) who used PHOENICS to simulate turbulent flow in a stand of widely spaced trees. The CFD package used in the present study, FLUENT, did not have any modification made to its programming or to the turbulence models used. The purpose of this chapter is to determine if FLUENT, using K-S turbulence models, is capable of predicting airflow in a complex domain with a setup similar to the UBC wind tunnel. This chapter will also determine whether the K-S turbulence model can be employed and successfully run in FLUENT for study in subsequent chapters. 3.2. Results and Discussion Given the complexity of the wind tunnel domain used, the CFD program and turbulence models performed well. A total of eight different simulations of the solid tree canopy were analysed. These simulations all dealt with software settings and not the domain layout. As stated earlier, three different variations of the K-S turbulence model were used, each simulated for a minimal number of iterations until convergence, as well as a longer amount of time to a point where residuals had levelled off. Using the standard K-S turbulence model, unsteady simulations 66 of the domain were also undertaken: one equivalent to a 4.5 second flow and the other a 6 second flow. The unsteady flows were found to be equivalent to each other and equivalent to the longer iterated steady state standard K-S turbulence model. Analysis was mainly carried out on the horizontal velocity. Vertical velocity was not considered due to its much lower values in a primarily horizontal flow. Turbulent kinetic energy was not examined as it was expected that these values would have been greatly impacted by high shearing due to the solid tree design used. Turbulent kinetic energy will be examined more closely in subsequent chapters where porous canopy structures are employed. 3.2.1. Standard K-E turbulence model The two steady state standard K-e turbulence models were simulated and their output was analysed using several different methods and compared to wind tunnel results. The difference between the two models was the length of time for simulation, or number of iterations. The short one was simulated until all of the residuals became smaller than a specified threshold (R^ < 1 x 10"4) and the longer one was run so that the solution converged; residuals were no longer becoming any smaller. Between the short and long simulations, there was a slight difference in horizontal velocity (figures 3.1 and 3.2) and an even greater difference in turbulent kinetic energy (figure 3.3 and 3.4). The longer iterated simulation looked closer to the wind tunnel flow pattern (figure 3.5 a and 3.5 b). As one would expect with solid trees, the clearings downwind of forest edge have large areas of reversed horizontal flow that extends about five tree heights downwind of the edge. Both of these characteristics are seen in the short and long simulations. This was not observed in the wind tunnel because those trees were not solid objects. The longer 67 run simulation also seems to have a better aligning of wind profiles above the canopy and a smoother, more natural reaction to the clearing when compared to the wind tunnel. Visual comparisons of the turbulent kinetic energy profiles between the short and long simulations do show differences, both of which were different than what was seen in the wind tunnel (figure 3.5 b). Well above the stand and clearing, an area of higher TKE values are apparent in the short simulation but not the longer one, which is more similar to the wind tunnel. Also in the long simulation, an area of higher TKE values are seen at canopy height in the clearing and increase and peak just above the first two trees downwind of the clearing (figure 3.4). This pattern is seen at a much lesser scale in the wind tunnel where a peak in TKE was observed just above and downwind of the first edge (figure 3.5 b). Scatterplots show how well the predicted points compared with the observed wind tunnel values. Figure 3.6 show that both simulations had poorer results at the lower horizontal wind speeds which may be indicative of the points in or just above the canopy. With faster wind speeds, the longer iterated simulation had a tighter clustering of points. Vertical profiles of horizontal velocity are probably one of the most useful tools that show the differences between models and wind tunnel observations. They also indicate where the models over and under predicted conditions. For the longer simulated K-S model, above the canopy, the profiles at the two upwind locations (figure 3.7 a and b) are slightly different compared to the observed data. The modelled velocity increases rapidly over a short elevation then increases more slowly for the rest of the profile. In both upwind profiles, horizontal wind speed was over predicted above stand height. Even though wind speed was over predicted, the standard K-S turbulence model did appear to have a vertical windshear, indicated by the slope that was very similar to that of the observed data. 68 The six locations in the clearing (figure 3.7 c-h) continued to show a similar trend seen in the first two locations. At and above 1.5 z/ht the modelled velocity profile became parallel with the wind tunnel observations. Also, as distance from the edge increased, so did the minimum speeds predicted close to the ground. The last two locations, within the second stand, had a similar pattern seen within the first stand. At location i (figure 3.7 i) the model over predicts, and has a velocity profile with a slope parallel to the wind tunnel profile, but at a greater speed. The shorter run K-S turbulence model did seem to have some anomalies based on the quantitative results (table 3.1). Compared to the other simulations with solid tree canopies, the results for the "quicker" model show a much different set of statistics. Upon initial impression, the results seem quite fair, considering that the trees are solid objects. The root mean square error for the horizontal velocity was 1.20 ms"', with the systematic RMSE (RMSES) 0.702 ms"' and the unsystematic (RMSEU) 0.974 ms"1. The Willmott d score was 0.937, which was the highest for the group. These values do not seem to correspond well with what was seen in the scatterplot or the velocity profiles. The r2 value was only 0.868, a value that does seem more consistent with what was seen in the qualitative analysis. The longer run model seemed to have quantitative results that were more inline with the other models. This model had an RMSE of 1.35 ms"1 and with a higher systematic error (RMSES) of 1.21 ms"1 and a RMSEU of 0.611 ms"', indicating that more refinement of the model could have occurred. The Willmott d score was 0.931, which again is good considering that a solid canopy was used. The unsteady simulations using the standard K-S turbulence models show that the higher iterated steady state simulations had a more thorough solution than the lesser iterated simulations and that the unsteady simulations perform similarly to steady state ones (table 3.1). 69 3.2.2. Renormalized Group K-S turbulence model The Renormalized Group (RNG) K-e turbulence model had a short simulation run at just under 3000 iterations and a longer simulation at about 13000 iterations. Both short and long appear similar. The most noticeable differences is that the long simulation has a higher horizontal wind speed well above the forest stand compared to the short simulation (figure 3.8 and 3.9). Another difference between the two simulations has been seen in the other two turbulence models: a smoother profile over the first stand for the longer simulated runs. Slight changes in the isotachs can be seen over the first stand in the shorter run, this anomaly is not present in the longer simulated case. Horizontal velocity profiles above the clearing seem similar in both short and long simulations. Compared to other turbulence models, the area of reversed horizontal flow downwind of the first stand seems shorter at about 4 x/ht. The TKE appears to be quite similar in both short (figure 3.10) and long (figure 3.11) simulations in the RNG K-S model. The only difference, as seen with most of the TKE profiles so far, is the higher value in the longer simulation which appear above the beginning of the first stand. The scatterplots for the two steady state simulations using the renormalized group turbulence model, do have a slight difference (figure 3.12). The difference between the two is much less than the standard K-S turbulence models. Again, the difference is seen at the higher wind speeds. The velocity scattering of the two models does show that the longer simulation is slightly less scattered around the best fit line. The vertical profiles of horizontal velocity for the RNG turbulence model are quite similar to those previously described from the standard K-S turbulence model (refer to figure 3.7). Results of the quantitative analysis for the RNG K-S simulation fell between the standard 70 and realizable K-S models (table 3.1). The shorter simulation had a Willmott d score of 0.913 and a RMSE of 1.56 ms 1 . The RMSES was 1.39 ms 1 and the RMSEU was 0.723 ms 1 . For the longer simulation, the Willmott d score was better at 0.930. The RMSE was better too at 1.37 ms"', but despite the better scores, the systematic RMSE still comprised a greater portion of that error at 1.23 ms"1, while the unsystematic RMSE was 0.611 ms"1. These values are very similar to those from the standard K-e turbulence model, which can be seen in the similar vertical profiles of horizontal velocity. 3.2.3. Realizable K-S turbulence model Two simulation lengths are also discussed for the realizable K-e turbulence model. The short simulation, converged just before 2500 iterations; at this point the convergence criteria were modified (lowered) so that the simulation would continue iterating until the residuals levelled off. The long simulation ended before reaching 12500 iterations. With this turbulence model, the two horizontal velocity profiles appeared much more similar to one another (figure 3.13 and 3.14) compared to the velocity profiles from the two standard K-e simulations. The longer run model did appear to have a more organised velocity profile where the interval between the speeds seemed more predictable/regular. The shorter run simulation was more influenced by the clearing than the longer run simulation. As with the standard K-S simulation, reversed horizontal flows were seen for about 5 ht downwind of the edges. The TKE for the short (figure 3.15) and long (figure 3.16) simulations involving the realizable K-e turbulence model also seem fairly similar. A larger difference is seen between the long realizable K-e profile and the long standard K-e profile, as the standard K-e profile seems to predict higher values of TKE within the clearing and above the second stand. 71 Scatterplots for the realizable K-e turbulence models were fairly similar to one another (figure 3.17 a and b). The longer iterated simulation seemed to have less scatter when plotted against the observed values. The vertical profiles of horizontal velocity showed very similar plots to the other two turbulence models (figure 3.7). One could argue that the realizable turbulence model performed slightly worse than the other two models. The lesser iterated realizable model appeared to perform worse as it over predicted wind tunnel observations by a greater margin. Both the short and long run lengths for the realizable K-S produced some of the poorest results in the statistical analysis. The shorter run simulation had the highest RMSE at 1.92 ms"1. The systematic and unsystematic components of the RMSE were equally as poor with the systematic being quite higher than the unsystematic at 1.72 ms"1 and 0.851 ms'1, respectively. The Willmott d score was the lowest in the field at 0.882. The longer run realizable K-e simulation performed better than the shorter one, but still not as well as the other turbulence models. The long run had a RMSE of 1.53 ms"1 with the systematic RMSE being 1.3 ms"' and an unsystematic RMSE of 0.724 ms"1. The Willmott d score for the longer realizable K-S simulation was 0.917. 3.2.4. Combined summary of model performance A Taylor diagram (Taylor, 2000) provides a convenient way of illustrating all of the simulations together in one diagram. In this chapter, it illustrates how well each of the short and long iterated simulations, from each turbulence model, compare to the wind tunnel observations ("Ref.") and between one another. The Taylor diagram in figure 3.18 visually shows what has been stated previously (and in table 3.1), that the longer iterated simulations generally performed 72 better than the shorter ones, and of those the standard and RNG K-e models performed slightly better than the realizable K-£ turbulence model. 3.3. Conclusion Considering that the forest canopy was described as solid object, the results of these simulations look promising. Given the complexity of the wind tunnel domain and the individual tree elements in the canopy, analysis showed that FLUENT and the turbulence models performed well. As stated earlier, three different variations of the K-e turbulence model were used, each simulated until convergence as well as a longer amount of time to a point where residuals had levelled off (no longer decreasing). Better results occurred when convergence criteria were adjusted to values that would allow residuals to level out. The unsteady simulations showed that there was no difference in extending simulation time from 4.5 seconds to 6 seconds. 73 1 3 4 5 6 7 8 9 Horizontal wind velocity (ms1) 10 11 12 ummMm 13 Figure 3.1: Isotachs of horizontal velocity (flowing from left to right) from the shorter iterated standard K-£ turbulence model are presented. The grey vertical lines indicate the ten sampling locations (A to J - refer back to figure 2.21 for explanation) where simulation output can be more directly compared to wind tunnel data (vertical profiles or quantitative analysis). White areas, exterior to the canopy elements, indicate reverse flow (negative horizontal velocity). JmmMlwmmaMM 0 1 2 3 4 5 6 7 8 9 Horizontal wind velocity (ms1) 10 11 12 13 Figure 3.2: Filled isotachs of the converged, longer iterated K-E turbulence model simulation have a slightly smoother profile over the forest unlike the isotach intervals from the shorter one (figure 3.1). This profile compares much better with the isotachs from the wind tunnel (figure 3.5 a). As in figure 3.1, white areas exterior to the canopy elements indicate a reverse horizontal velocity. MMlfflWlUfflMlMfflmimm 0 <3\ 1 2 3 4 5 6 7 8 9 2 2 Turbulent kinetic energy (TKE) (m" s ) 10 imrnmnTwrmtf Figure 3.3: Turbulent kinetic energy from the lower iterated K-E turbulence model simulation. White areas within the canopy indicate that no domain is present and therefore no flow is permitted through the canopy. A large area of TKE is seen over the stands and clearing. This pattern seems unrealistic and is strong evidence that the simulation needs to be iterated longer. ^ | W | ^ a 3 ^ u ^ u < | ^ p | f ^ ^ y K f ^ ^ | S | f l | ^ | ^ ^ & ^ ^ g ^ 0 1 2 3 4 5 6 7 8 9 Turbulent kinetic energy (TKE) (m"V2) 10 'YWYYYYYYYYYYYYYYYYW^k Figure 3.4: Similar to figure 3.3, this figure shows TKE from the longer iterated simulation using the standard K-e turbulence model. The TKE from this simulation is much lower above the clearing and stands, but notably higher above the upwind edge of the second stand. TOOTHS 0 00 •12 0 -12 -4 0 4 Distance (x/h t ) 8 -8 -4 0 4 Distance (x/h t ) 8 Calctke: 0.2 0.4 0.6 0.8 1 1.21.41.61.8 2 2.22.4 -8 1 1 . 5 2 2.5 3 3 . 5 4 4,5 5 5.5 6 6.5 7 7.5 8 12 12 B Figure 3.5: Observations from the wind tunnel were interpolated, plotted and isolines drawn from the data. In the plots above, the most upwind extent (-12 x/ht) would represent location "A" in figure 2.21; the most downwind extent (+14.25 x/ht) would represent location "J". Dashed lines show the extent of the canopy. A) isotach of horizontal velocity show that the flow over the upwind stand is relatively streamlined with a tighter gradient close to the top of the forest. Above the upwind portion of the clearing wind speed increases, this increase becomes more pronounced with increasing height. Ahead of the second stand, there is a sudden decrease in horizontal velocity. The intensity of the decrease above the height of the canopy is mitigated with increasing height. Above the second stand, the speed increases as the flow adjusts to the new surface conditions. Also in the second stand, stem flow can be seen close to the ground. B) Isolines of turbulent kinetic energy calculated from wind tunnel measurements yield an interesting pattern. The area of greatest turbulence values is seen above and downwind of the first stand. There is also tight gradient that bounds the forest perimeters. © .S> p £ 0.60-rg0.45-r So.30- r £0.15t x 0.15 _JMeanu: p £0.60 O0 4 5 ^ So.30-r 0) Scatterplot of Horizontal Windspeed Observed vs. Modelled (Steady std K-E) Scatterplot of Horizontal Windspeed Observed vs. Modelled (Steady std K-E after 21000it) 10.0 jT y _«*• 10.0 S s " / ~- 6.0 • •• •• £•• • • Arm & o.o . •//•? ' A Regression line 1:1 line 1 0.0 2.0 4.0 6.0 8.0 10.0 Measured (ms"1) 0.0 2.0 4.0 6.0 8.0 10.0 Measured (ms"') Figure 3.6: Scatterplots of the shorter iterated (a) and long iterated (b) simulations show inconsistent trends at lower speeds, which were likely from those points extracted within or just above the canopy. The scatterplot with the shorter iterated simulation shows a fair degree of scatter at the higher speeds, while the converged simulation has a tighter profile at the faster speeds. 79 Figure 3.7: Vertical profile plots of horizontal velocity for short and long simulations of all three K-E turbulence models from the domain with the solid tree shaped elements, as well as wind tunnel observations. Positions of each location (a-j) are in relation to the upwind edge of the clearing and can be visually referenced in figure 2.21. 6.073 6.517 5.173 5.173 5.173 5.173 5.173 5.173 std 4.5 s std 6.0 s RNG RNG (long) Realizable Real, (long) 6.175 6.250 6.103 6.095 6.091 5.173 std (long) 5.762 P 5.173 0 std K-6 2.116 2.116 2.116 2.116 2.116 2.116 2.116 2.116 *o 3.100 3.300 3.017 3.079 2.983 2.983 2.965 2.680 P S 278 278 278 278 278 278 278 278 N -1.193 -1.267 -1.150 -1.072 -1.032 -1.048 1.424 1.505 1.396 1.415 1.379 1.381 1.372 1.180 -0.341 -1.003 m b 1.445 1.803 1.291 1.460 1.291 1.281 1.272 0.977 MAE 1.346 1.717 1.231 1.386 1.229 1.225 1.209 0.702 RMSE, 0.724 0.851 0.611 0.723 0.617 0.605 0.611 0.974 RMSEU 1.528 1.917 1.374 1.564 1.375 1.366 1.355 1.200 RMSE 0.917 0.882 0.930 0.913 0.929 0.930 0.931 0.937 d 0.945 0.933 0.959 0.945 0.957 0.959 0.958 0.868 r2 Table 3.1: Statistics for horizontal velocity for both short and long iterated steady state simulations for each ofjhe K-S turbulence models, along with unsteady 4.5 and 6.0 second_simulations using the standard (std) K-8 turbulence model. In the table below, O is the averaged wind speed from all wind tunnel observations. P is the average wind speed output from sampled points from the simulation. The standard deviation for the observed data is noted by s„ and sp for the model output. N is the number of sampled points. The y-intercept is represented by b and m is the slope of the best fit line. MAE stands for Mean Absolute Error. RMSE is the Root Mean-Square Error with RMSES and RSMEU being the systematic and unsystematic components of the RSME, respectively. The Willmott d score, ranging from 0 to 1, is indicated by d, and finally, the correlation coefficient is r2. 00 iiiL?ii.uw«»jj^u.uu , " ' n 1 " 2 m ' . 3 -" m'«- ' ^ ». . .m„. 10 iii-imijiimiLiii-imu.^ij^jjiii 4 5 6 7 8 9 Horizontal wind velocity (ms1) .iji.junu^'ii 11 13 TPW^1 12 ,.ji-Ui—. » u.imm p^n^^j ;|p. Figure 3.8: Isotachs of horizontal velocity. Similar setup to the simulation seen in figure 3.1 with the exception that this output was from the shorter iterated simulation using the RNG K-E turbulence model. ju, 0 00 1 3 4 5 6 7 8 9 Horizontal wind velocity (ms1) 10 11 12 13 Figure 3.9: Isotachs of horizontal velocity, this image is similar to figure 3.8 (RNG K-£ turbulence model used) except that it was simulated for more iterations compared to the previous figure. 0 oo -1^ 1 _ 2 3 4 5 6 7 8 9 Turbulent kinetic energy (TKE) (m"2s2) 10 .mmmmmrrm^ Figure 3.10: Turbulent kinetic energy isolines from the lower iterated RNG K-£ turbulence model simulation. Aside from the turbulence model used, the setup for this simulation was very similar to the simulation used to generate figure 3.3. TTTTTHITlTrn7TTTTm^^ 0 00 2 3 4 5 6 7 8 9 Turbulent kinetic energy (TKE) (m"2s2) 10 mmnrvmnm Figure 3.11: The main difference between the simulation from which this figure was produced from and the simulation from figure 3.10 was the number of iterations. This simulation was also used the RNG K-S turbulence model. Isolines of turbulent kinetic energy are shown. rmrmTTTntTrmrmmTTTTTTTm 1 Scatterplot of Horizontal Windspeed Observed vs. Modelled (Steady RNG K-E) Scatterplot of Horizontal Windspeed Observed vs. Modelled (Steady RNG K-E after 13000it) .# / - •sWi* x / fllfflr jdtier •• /' r f i i r i 2lr^ ^ • \ f i •*' ** ^ _~ 6.0 'w G £ *£'J' • • *f / ' > 1 • 1 > • ' - " rr* Atif 8.0 4tMr / *ar/ 8.0 Jy / 7\ '' »• iai ipprr s / • <*Jir ^Jmr JABS S/ ~3rV / . F*s s / 6.0 4.0 • • •» J v/y ***/ • i / • / 4.0 ,// •SM+Am ¥/ 1 0. y» 2.0 /?* 0.0 2.0 4.0 •— Regression line 1:1 line 6.0 8.0 0.0 10.0 / / ' • * 7%» . ,* s . /. . • 2.0 /' * / Regression line 1:1 line 4.0 8.0 10.0 Measured (ms1) Measured (ms") Figure 3.17: Scatterplots of the short iterated (a) and long iterated (b) simulations using the realizable K-E turbulence model. These two plots exhibit similar trends seen in the other two turbulence models. The two plots are fairly similar, except for a tighter clustering of points at the higher windspeeds. 91 o o K-e model cr> _c>' 1.50 Short Long • Standard RNG Realizable 1 x/ht upwind. The other exception is with the divergence as velocity increases with height. It should be noted that as distance from the edge increases the height at which the divergence occurs at increases, and the degree of difference between wind tunnel and model velocity decreases. This trend continues in the second half of the clearing. At each of these three locations, the standard K-S model seems to over predict horizontal velocity the least. Renormalized group K-£ over predicts only marginally more than the standard model, while realizable over predicts the most, but still not by much. It is interesting to note that the degree to which one model over predicts compared to another is roughly the same, for example the difference in over prediction between the RNG K-S model and the standard K-E model, is very similar to the difference seen between the realizable K-S model and the RNG K-S model. Similarities between wind tunnel observations and model output are apparent at locations F to H (figure 4.15 f-h). Slight divergences are seen between wind tunnel and model velocities, which are exacerbated with the edge from the second stand. This is especially seen at location H. The same differences between K-S models seen at locations C to E are also identifiable in locations F to H. The final two locations (I and J) in the second stand show similar traits seen in the first two locations from the first stand (figure 4.15 i, j). Divergence between the wind tunnel and model profiles can be seen and the point at which divergence occurs depends on the distance from the upwind edge. Divergence begins at an earlier height at the location closer to the edge (I). One aspect the model did not account for as hoped, was stem flow. While a very slight increase in velocity can be seen close to the ground at location I, it was not as apparent as it was in the wind 109 tunnel. 4.2.2.4. Vertical profile descriptions of turbulent kinetic energy Vertical profiles of all the locations clearly show the discrepancy between the turbulent kinetic energy values calculated from wind tunnel observations and what was predicted by the standard K-S turbulence model and to lesser extents, the renormalized group and realizable K-e models. One difference seen between wind tunnel and simulations that relate to all locations, not just at specific locations: wind tunnel TKE values only seem to reduce slightly after reaching their peak values at each location. This differs from the simulated output that sees a large reduction in TKE with increased height after reaching their maximum TKE value. Regardless of surface type (forest or clearing) the models show little difference at ground level in TKE values between clearing and stands. The wind tunnel model had a distinct absence of TKE within most of the stands, with the exception of the very top of the stand. The notable difference was at location I which did have a small TKE value in the stem area close to the ground. The differences between locations A and B in the wind tunnel were somewhat similar, although there was a slight difference between the increase in TKE close to the top of the canopy (figure 4.16 a, b). The output from the model at the same locations show a large increase in TKE, which illustrates that the upwind roughness elements may not have served their intended purpose, to shorten the length of the required forested section. Location C in the wind tunnel showed evidence of the upwind stand (figure 4.16 c), as the TKE profile under 1 z/ht looked very similar to locations A and B. An increase in TKE, upwards from the ground, for both wind tunnel and model locations, were somewhat similar for locations 110 D to H (figure 4.16 d-h). Peak TKE from those locations seemed to increase in height with each subsequent location downwind. This was easy to see in the modelled output, and somewhat harder to determine with the wind tunnel profiles due to the much smaller values, but still present. Locations I and J in the second stand have similar profiles to A and B from the first stand (figure 4.16 i, j and a, b). With the exception of the stem area in location I, the wind tunnel values are very similar between first and second stands. Between K-S models, the realizable K-S model was superior to the standard model in all cases and superior to the renormalized group model in all but one instance (figure 4.16 j). At all locations, the realizable K-S model seemed to adjust more quickly to the changes in TKE compared with the renormalized group model and especially with the standard model. Only a small discrepancy existed between the realizable and renormalized group models; however, a very large difference was seen with the standard K-S model. 4.2.2.5. Quantitative description Quantitative statistical analyses were performed on the output of the three K-e turbulence models. The values used were horizontal wind speed and turbulent kinetic energy. Table 4.3 shows the various statistics derived from the horizontal velocity output of the model and wind tunnel used to validate the models ability to predict horizontal wind speed. Table 4.4 shows the same but for turbulent kinetic energy. Table 4.3 shows that for horizontal wind velocity the standard K-S model output had values that agreed more with the wind tunnel measurements. The main difference was with the r2 values, which were very good for all models, show that the realizable K-S model had less 111 variability than the other two. Table 4.4, unlike the table for horizontal velocity, show that the standard K-E model was the poorest at accurately predicting TKE. The realizable K-E turbulence model performed the best out of three models but still had a low Willmott d score of 0.384, not a good value for model validation, but expected. 4.2.2.6. Discussion of layered tree stand The results did not conclusively point to a single K-S turbulence model being superior for all statistics examined. For horizontal wind speed, the standard K-e model seemed to perform best overall. Examining the vertical profiles and quantitative results of horizontal velocity, it would seem that the standard K-S turbulence model performed best. The standard K-S model seemed to track best when compared to the other two models when it came to the vertical profile plots. The quantitative analysis showed the highest Willmott d score, lowest RMSE value and the best RMSEU to RMSES quotient were all from the standard K-S output. The other two K-E turbulence models did have tighter clustering of points on the scatterplots. This pattern was seen with higher r2 values, the highest being with the realizable K-E model. Additionally, the quantitative analysis, and the scatterplots to a less extent, showed that both the y-intercept and slope of the standard K-S model were closer to 0 and 1, respectively, than the other two models. The seemingly better validation for the standard K-e model did not carry over to the turbulent kinetic energy analysis. In all the analyses performed, the standard K-S model was inferior to the other two K-S models and in some instances, much worse. By including both velocity and TKE, statistics showed that out of the three K-S turbulence models used, the realizable K-E model performed best. It had better statistical scores in most categories: lower 112 RMSE value, higher Willmott d score, better RMSEU to RMSES quotient and better y-intercept. It was very close to having the same slope as the RNG model. The r2 value was the only statistic where the RNG K-e model clearly performed better. 4.2.3. Three dimensional layered blockstands Three dimensional simulations were performed with a three dimensional domain using a layered blockstand style of canopy. Most qualitative analysis will show that the output from the three dimensional simulations performed very similarly to the two dimensional simulations. Output from the simulation will be compared with the wind tunnel data. Differences between the three dimensional and two dimensional output will also be highlighted. 4.2.3.1. Qualitative description of profile Horizontal isotachs of a plane sliced in the middle of the narrow width three dimensional wind tunnel show some differences from the two dimensional simulation. Closely examining the side profile of the three dimensional simulations reveals that the standard K-e turbulence model (figure 4.17) still slightly over predicted the wind speeds (figure 3.5 a) at upper levels sampled in the wind tunnel. A decrease in wind speed is seen above 1 z/ht in the clearing, similar to the findings in the two dimensional domains. Differences between the output from the two and three dimensional domains were more apparent with TKE (figure 4.18). Considering that the two dimensional blockstand simulation greatly over predicted TKE, an over prediction was also expected in the three dimensional domain. Turbulent kinetic energy for the three dimensional simulation was over predicted in most areas in and above the stands by a margin of about 4 - 5 times the values calculated in the 113 wind tunnel. Areas of peak output of TKE in the three dimensional simulation were similar to those areas seen in the two dimensional simulations 4.2.3.2. Scatterplot description Horizontal velocity output from the three dimensional simulation was plotted with wind speed data collected from the wind tunnel (figure 4.19). The scatterplot showed similar patterns seen in the two dimensional simulations: slight under prediction from 2 to 6 ms"1 and slight over prediction at horizontal velocities > 8 ms"1. 4.2.3.3. Velocity Profile and description Vertical profiles of horizontal velocity from the three dimensional simulation were quite similar to the output from the two dimensional simulations (figure 4.20). Where there were differences, output from the three dimensional simulation over predicted wind tunnel velocity slightly less than the output from the two dimensional simulation. 4.2.3.4. Turbulent kinetic energy profile and description Like the vertical profiles of horizontal velocity, the vertical profiles of TKE from the three dimensional domain (figure 4.21) were very similar to those from the two dimensional domain, with one exception: magnitude. At each location the two and three dimensional simulation profiles may have had similar values close to the ground but that difference grew with increasing height up until the maximum TKE value was attained for each specific location (excluding locations I and J). After peak TKE, the values gradually converge closer to each other towards the top of the sampling profile (excluding locations G to J). Locations I and J seem to have the 114 greatest difference at the top of the profile, while locations G and H seem to maintain the difference between the TKE values of the two and three dimensional simulations. 4.2.3.5. Quantitative description Quantitative analysis performed on the output from the three dimensional simulation showed fairly good results (table 4.5), better than the two dimensional output in some instances. Compared to the quantitative analysis from the two dimensional standard K-e simulation, the analysis from the three dimensional simulation indicate statistical values that are more comparable to the wind tunnel observations. Quantitative analysis on turbulent kinetic energy showed poor values (table 4.6). These values were poorer than what was simulated with the two dimensional domain. The poor values were expected given the analysis results of the two dimensional simulations and vertical profiles seen in (figure 4.21). The slope of the regression line was closer to 1 than for the two dimensional simulation. 4.2.3.6. Discussion of the three dimensional layered stand The three dimensional domain with a layered stand, was time consuming to setup and simulate. Compared with the wind tunnel data, the horizontal velocity output from the three dimensional simulation in many cases was predicted quite well, slightly better than the output from the similar two dimensional simulation. Compared to the two dimensional domain, the simulation from the three dimensional domain had a slightly better Willmott d score, moderately better RMSE and a markedly improved RMSEU to RMSES quotient. The over prediction of horizontal wind speed above the stands is still an issue. 115 The accuracy of turbulent kinetic energy continued to be a problem with this turbulence model, or at least the application of it in this setup. The TKE output from the three dimensional simulation was even higher than what was over predicted in the two dimensional porous blockstand domain using the same standard K-S turbulence model. 4.2.4. Layered tree stand (Blockstand) in 3 m tall domain In both the porous tree stands and layered tree stands, it has been noted that both the horizontal velocity and turbulent kinetic energy were over predicted compared to the wind tunnel values, especially higher in the domain. It was thought that this may have been due to the roughness elements and trees themselves "filling up" the space that the air was flowing through and causing the speed to increase to allow the same amount of air to move through the domain. It was believed that if this was the cause then a larger domain volume would therefore make the roughness elements and trees a smaller portion of the overall flow area and thus have a lesser impact on wind speed and turbulence. Qualitative analyses for these simulations were limited to vertical profiles of horizontal velocity and TKE at ten locations. Quantitative statistics were also performed to note the differences between the two domains. 4.2.4.1. Vertical profiles of horizontal velocity between two domain heights Plots of horizontal velocity for the two domains were evaluated. The output from the normal and tall (3 m) domains using the standard K-S turbulence models were used as the basis for the comparison. Wind tunnel data was also included on the plots to show how well each domain performed in comparison with wind tunnel data. 116 The vertical profiles at the ten locations (figure 4.22 a-j) were categorised into four groups: the upwind stand (locations A and B), the first half of the clearing (locations C to E), the second half of the clearing (locations F through H) and the downwind stand (locations I and J). The first two locations (figure 4.22 a, b) from the wind tunnel appear to have been better matched by the output from the taller domain for the upper points sampled, whereas the regular height domain over predicted horizontal velocity. Wind tunnel points within and just above the canopy (up to two tree heights for location B) seem to be better matched with the normal size domain compared to the taller domain. Simulation output for locations C to E, in the first half of the clearing (figure 4.22 c-e), show that the taller domain did not perform as well as the normal sized domain. Only the top quarter of sampling heights from the wind tunnel is where the taller domain performed better against the standard domain with respect to the wind tunnel observations. As distance form edge increases, the horizontal velocity profile become more similar to the regular domain height. The standard domain height simulation also performed better in the second half of the clearing (locations F to H). The wind tunnel profile was more closely followed by the output from the simulation with standard height domain, through most of the sampled heights, especially at locations F and G (figure 4.22 f, g). At location H (figure 4.22 h), the wind tunnel data were mostly situated between the output from the two domain profiles. Above 3 z/ht, the taller domain at these three locations has output more closely associated with wind tunnel observations, but the differences were small. At the last two locations (figure 4.22 i, j), there was a similarity seen between these profiles and those from locations A and B (figure 4.22 a, b). Within the canopy, simulations from both domains performed fairly well to the wind tunnel observations. From the canopy top 117 to 2 z/ht, the standard height domain performed better than the taller domain. Above 2 z/ht, the wind tunnel speed began to slow with height and became more aligned with the profile from the taller domain. 4.2.4.2. Vertical profiles of TKE between the standard and 3 m domain Vertical profiles for TKE were also examined between the standard and 3 m domains. The profiles will show that there was no difference in shape between the two domains, but there was a difference in magnitude (figure 4.23 a-j). Due to the similarities between the two simulated profiles, all locations will be presented together. The vertical profiles between the wind tunnel TKE and simulated TKE were quite different. Within the stand, the observed TKE rises quickly at the top of the canopy then mostly maintains the same degree of turbulence, with some degradation, through the rest of the profile. The simulated TKE had much different characteristics. As with the observed values, there was a quick rise and peak just above the top of the canopy, but at larger magnitudes. Both modelled TKE domains saw large reductions of TKE after their peak TKE values were reached. The clearing (figure 4.23 c-h) was where wind tunnel TKE and modelled TKE were most similar, but the magnitude of TKE was still very great. The difference between the standard height domain and the 3 m tall domain were quite noticeable. In every instance, the simulation from the standard height over predicted TKE more than in the tall domain. These over predictions were greatest in the profile where the TKE values peaked at each location. Close to the ground the two domains had similar TKE values, as well as towards the upper portion of each location. 118 4.2.4.3. Quantitative comparison Quantitative analysis for the 3 m tall domain is presented and the analysis from the standard height domain is posted again to show the differences between the two domains. Table 4.7 shows the analysis for horizontal velocity for the two domains. Both domains performed fairly well, with the standard domain having a better systematic and unsystematic root mean square values of 0.463 and 0.552, respectively, while the 3 m tall domain had values of 0.543 for the systematic RMSE and 0.489 for the unsystematic RMSE. The overall RMSE for the 3m tall domain was slightly larger at 0.730 versus 0.720 for the standard height domain. Willmott d scores were nearly even with the standard height domain fairing a marginally better agreement of 0.976, while the taller domain had a Willmott d score of 0.971. One interesting statistic to point out was the slope (m) of the 3 m tall domain, which was nearly 1. Analysis of turbulent kinetic energy from the two domains compared with the calculated value from the wind tunnel show a larger difference between the two domains. Most of the statistics in Table 4.8 show that the 3 m domain performed better than the standard height domain. Though the TKE statistics were quite poor as a whole, accuracy of K-S turbulence models in canopy settings models are quite often not as good as other indicators like velocity (Foudhill et al., 2005 and Liu et al., 1996). The RMSE for the standard domain was about 1/3 higher than the 3 m tall domain. Willmott d and r2 from the taller domain were both better with values of 0.277 and 0.275, respectively. Alternatively, the standard domain had a Willmott d score of 0.202 and an r2 value of 0.119. 119 4.2.4.4. Discussion of 3m tall domain Model results from the taller domain indicate that the presence of the forest stands may have contributed to the over prediction of certain variables. Though increasing the domain height did help to explain the differences, it did not make for a better domain setup. Reviewing the output from the models, it appears that within and just above the canopy, the domain that was of equal height to the wind tunnel (i.e. the standard height domain) performed as good or better than the domain that was made twice the height of the wind tunnel. Furthermore, from the six locations within the clearing, the output from the standard height domain out performed the taller domain through most of the observed profile. These factors led to better statistical scores for horizontal wind speed. The taller domain did seem to provide better output with respect to turbulent kinetic energy. This was likely wholly due to the decrease in magnitude that was seen in the vertical profiles of TKE (figure 4.23 a-j). 4.2.5. Miscellaneous domain/canopy simulations Additional simulations were performed with different attributes of the mesh and domain. In total the results of three additional simulations will be presented. These were based upon domain layout and mesh attributes along with length of simulation. The additional simulations were done to determine if the roughness elements were necessary instead of a large canopy and if a boundary layer grid system would have any additional benefits to the domain. A longer simulation was also done to ensure that more simulation time had no additional impact on the output. 120 4.2.5.1. Stands only, no roughness elements with and without modified grid This simulation was used to determine if it was necessary to recreate the upwind roughness elements seen in the wind tunnel. In the wind tunnel, those roughness elements were present in order to cut down on the size of the model forest. In creating the numerical domain, roughness elements were just as time consuming as making an entire stand of trees, and it was much easier to forgo the process with a blockstand configuration (figure 2.11). Also included for comparison and grid independence was a finer boundary layer type grid, which was used along with the stand only layout. Table 4.9 shows the statistical analyses of the different simulations including the standard layout. These statistics show that the absence of the bluff bodies and roughness elements result in better agreement between the wind tunnel flow and modelled flow as evidenced by the d score and RMSE value. The domain that had boundary layer grid setup performed very similar to the normal grid, but no clear differences were found. Similarly, it is shown in the table that there were no advantages by simulating the models for a longer period of time. Compared to the standard layout, the y-intercept was much closer to 0 in the all canopy layout and the all canopy layout with the boundary layer grid with -0.168 and -0.265, respectively. The slope (m) was nearly 1 for both all canopy (0.996) and boundary layer grid (1.014) domains. All of the RMSE values including the RMSEU to RMSES quotient were very good. RMSE for the all canopy layout was 0.552 with unsystematic RMSE being 0.518 and systematic RMSE only 0.190. The Willmott d statistic was also very high at 0.983 for all canopy only domains. The only slight advantage that the standard layout had over the canopy only domains was a slightly higher r2 score. While the simulations showed excellent results for horizontal velocity in the all canopy 121 domains, the same could not be said for turbulent kinetic energy. Statistics for turbulent kinetic energy were worse than those from the standard K-S turbulence model, which itself were already considered poor compared to the other K-S turbulence models. The canopy-only domains all performed poorly even with the boundary layer grid settings and the longer simulation time. The only acceptable statistics that could be pulled from these results were from the slope, which had values close to 1, with 1.276 for the normal grid all canopy layout and 1.114 for the boundary layer grid settings. Full statistics for turbulent kinetic energy from these simulations can be seen in table 4.10. 4.2.5.2. Discussion of miscellaneous domain/canopy simulations The results from the simulations that had no roughness elements, only an extended canopy stand were originally promising. Both the "all canopy stand" simulation and the "all canopy stand with boundary layer grid" simulation showed very good results for predicting horizontal wind speed. The statistics presented showed that there was very little difference between incorporating a boundary layer grid system for these simulations. As in section 4.2.2.2.1., simulating the model for 6 seconds made no appreciable difference compared to the model simulated for 4 seconds. Turbulent kinetic energy was more over predicted than in the standard domain. 4.3. Discussion In this chapter, many simulations were performed with various canopy and domain 122 properties. The main objective of all these simulations was to assess, using wind tunnel data, the validity of using a commercial CFD program, FLUENT, to predict airflow through and above a forest, and suggest model configurations that perform best. Qualitative and quantitative analyses were performed on the output from simulations at their equivalent points in the wind tunnel. Wind tunnel points were sampled up to 0.59 m or roughly four tree heights and thus the analyses were performed up to that height as well. Many past studies presented analysis of data up to 2.0 to 2.5 tree heights (eg. Green, 1992; Li et al., 1990 and Foudhil et al., 2005). In order to compare results with the widest range of other studies, the bottom half of the locations sampled were also used (< 2 z/ht), which in some cases, the models simulated appear to agree much better. This also had an impact on the quantitative statistics performed, which were absent from many studies. It was seen in the vertical profiles that at output below two tree heights, horizontal wind speed generally performed very well and comparable to what was observed in the wind tunnel, compared to the output from below 4 tree heights. Discussion for this chapter will first include a presentation of the quantitative analyses at heights below two tree heights(z < 2 ht) which will be presented as the "half profile". These quantitative statistics will be compared to the statistics presented in the results section, which included output from the full height from the sampled locations (z < 4 ht) or the "full profile". Furthermore, qualitative comparison with other studies will mainly be done with output from the half profile (under two tree heights) but where available, with the full profile. A Taylor diagram of all eight simulations can be seen in figure 4.24. This diagram visibly shows how well all simulations performed compared to each other and with the wind tunnel observations. From the diagram, it would appear as though all of the simulations performed 123 quite well, even if some had slightly higher standard deviations. 4.3.1. Porous tree stand Below two tree heights, horizontal wind speed has superior statistical results compared to statistics from the full profile. Table 4.11 shows the statistics from the simulation that had the tree shaped porous canopy with both the full profile and the lower half of the profile. Some statistics did not change much with the elimination of the upper half of the profile, but others, like RMSES, RMSE and the y-intercept, did. This difference shows that the turbulence model and canopy settings of this setup performed fairly well within and just above the canopy. Overall the RNG turbulence model performed better than the standard K-e turbulence model when examining horizontal velocity below two tree heights. Turbulent kinetic energy predicted by the models was not nearly as representative as the horizontal wind velocity was regardless of the profile used (half versus full). Statistics from the lower half of the profile were generally worse than statistics from the entire full profile. Two exceptions were the slope, Willmott d statistic and r2 (table 4.12). The reasons for this could be seen in the vertical profiles from the locations above the stands, where large TKE values were reported above the stands (figure 4.7). These large values of TKE were much greater than those from the wind tunnel. 4.3.2. Layered tree stand (Blockstand) With a four tree height profile, the blockstand canopy seemed to be very similar to the porous tree shaped canopy using the same standard K-E turbulence model. The only clear advantage between the two was the RSMEU to RMSES quotient the blockstand simulation had. 124 When looking at the half height (profile) statistics (tables 4.11 vs. 4.13), the porous tree canopy performed slightly better than the porous blockstand, this was even with the disadvantage of the solid stems used on the lowest portion of the profile, especially at locations downwind of the clearing where stem flow was retarded. The porous tree design may have had the benefit of good reproduction of the wind tunnel trees. In the real world, such dimensions would not have been perfectly reproduced in a homogeneous stand. However, in a heterogeneous stand, using a porous tree shape may have possible real application as each tree could be authored with it's own shape, height and LAD characteristics. While making individual model trees, whether physical or numerical, is quite labour intensive, breaking up the classical homogeneous stands into sections with different heights and permeability have been investigated along edges (Dupont and Brunet, 2008). Three different types of K-S turbulence models were tested using the layered blockstand. The three turbulence models have not been altered nor their constant altered. Other authors of similar studies used the standard K-S turbulence model with various modification to account for the influence the canopies had on the flow. The porous media tool was utilised in the CFD program that allowed the canopy to interact with the flow based on LAD and Cd of the wind tunnel tree canopy. For horizontal wind speed, when looking at the analyses from other studies (eg. Foudhil et al., 2005 and Wilson and Flesch, 1999), output from the standard K-S turbulence model performed best out of the three models used for the full (4 z/ht) profile. Much like the porous tree simulation the half (2 z/ht) profile from the RNG K-S turbulence model simulation had the better statistics. Table 4.13 shows that most statistics improved when the top half of the profile was removed. The realizable K-S turbulence model performed poorest compared to the other two K-S turbulence models. 125 There were differences seen when examining the TKE values between full and half profiles of the three K-S turbulence models. Overall, the same over prediction that was seen in the tree shaped porous canopy was also seen in the block stand. As with the tree shaped model, the block stand models had the same over prediction just above the canopy top, followed by a dramatic decrease. Therefore one would expect that it would also cause poorer statistics, which it did in some instances. From figures 4.16 a-j, it was shown that the standard K-e model performed the poorest, with the greatest over prediction. This was also seen quantitatively when comparing the full and half profiles and with the realizable and RNG K-8 turbulence models. Statistics in table 4.14 show that the lower half of the profile did have worse output compared to the whole profile with respect to RSME and its components, MAE and at the y-intercept. Slope, r2 and the Willmott d statistic each had better scores when only considering the lower half of the profile. These differences may be due to the shape of the profile. Above two tree heights, observed TKE at all ten locations does not decrease very much from it's maximum value, whereas the simulated TKE values quite drastically decrease the further up in the profile one looks, thus the trend between the two are dissimilar. Under two tree heights, the models may have wildly over predicted TKE values, but the trend between observed and modelled were similar. Table 4.14 also shows that for TKE statistics for both the RNG and realizable K-e turbulence models, the bottom half of both profiles agreed much better to the observed values than the bottom half of the standard K-8 model did. Similar patterns between statistics of the full and half height profiles that were seen in the standard K-8 turbulence model were also seen in the realizable and RNG models. Again, better values for slope, Willmott d and r2, where found looking at the lower half of the profiles. Overall, they were better than similar statistics seen in the standard K-S model because that model over predicted and therefore the magnitude of 126 increase was so great and that the increases seen in the realizable and RNG models were more closer to the reality seen in the wind tunnel. Overall, as seen in the statistics of the full profiles and the vertical profiles (figures 4.16 a-j), the realizable K-S turbulence model performed best at heights below two tree heights at all ten sampling locations. 4.3.3. Comparison of flow traits with other studies The infiltration of higher wind speeds close to the ground as distance increased from the leeside of the clearing was seen in figures 4.1 and 4.4 has also been reported in many other studies most notably in Raynor (1971) and most recently in Dupont and Brunet (2008). The over prediction of velocity, particularly above stand height in the simulations performed was also seen in Wilson and Flesch (1999). In their study, a clearing 21.3 tree heights in length was simulated and they reported a slight over prediction of horizontal velocity within the clearing and above stand height. Foudhil et al. (2005) seemed to overcome this problem when they tackled the same scenario using their own modified K-S turbulence model with different model constants. Wind tunnel values of TKE, were calculated from the standard deviation of the wind components. Throughout the vertical profile, there was very little reduction in TKE above two tree heights in the wind tunnel, regardless of where the sampling location was. This is contradictory to other studies which do show similar peaks of TKE between one and two tree heights (eg. Frank and Ruck, 2008; Panferov and Sogachev, 2008; and Foudhil et al., 2005) and then report a noticeable decrease above two tree heights (Frank and Ruck, 2008; and Foudhil et al., 2005), in both their simulations and observations. A trait predicted in all K-S turbulence models and reported in other studies is the 127 development of TKE several tree heights downwind of the edge over the second stand. Both Dupont and Brunet (2008) and Wilson and Flesch (1999) reported that their models predicted increasing TKE values as distance downwind from the edge increased. This is opposite to what was reported in a similar numerical model by Frank and Ruck (2008) whose output showed an immediate rise in TKE following flow over a forest edge from a clearing, with TKE subsiding with increasing distance from the edge. Aside from TKE values from the standard K-S turbulence model, predicted values from the RNG and realizable K-S models compared well with those from other studies. Turbulent kinetic energy values predicted between one and two tree heights above the stands for the RNG and realizable K-e turbulence models were approximately 4 to 5 m2s"2. These values were comparable to those predicted in similar conditions by both Foudhil et al. (2005) (4-5 m2s"2) and Wilson and Flesch (1999) (4-6 mY 2 ). 4.3.4. Exit flows Two canopy types, tree shaped porous and block stand porous, using the standard K-S turbulence models, were more closely examined to compare flow characteristics along forest edges. We can look at forest edge flows in the same way as the exit flow schematic diagrams from Lee (2000), Chen et al. (1995) and Judd et al. (1996) (who were actually examining flow through porous fences) and discuss the output from the simulations in a similar way. These authors all broke the flow down in a number of fluid-mechanically distinct zones. The first zone contains the two upwind sampling locations (locations A and B) each of which show similar wind profiles (figures 4.4 a, b), which are pretty consistent with any forest stand with an equilibrated flow. The region upwind of the edge was called the approach flow, where the 128 various LAD profiles would impact each individual forest differently. Towards the top of the stand horizontal wind speed greatly increases and with the presence of a strong inflection point. Some models also observed a slight increase in the speed close to the ground due to the absence of foliage; this was also the case with the wind tunnel data used, albeit small. This would not be modelled in our simulation as solid stems prevented this in the domain using tree shaped elements due to the height at which our lowest points were where output was collected. The next distinct zone was the called the "quiet zone". This is an area that doesn't differ greatly from the conditions present within the stand (Lee, 2000). Wind speed is at or less than what was recorded in the stand. Both simulated blockstand and porous tree stands had very similar quiet zones, extending about 4ht downwind of the stand (locations C to E). Above and downwind of the quiet zone was the mixing or wake zone. In this zone (locations F and G), streamlines were angled slightly downward (eg. Raynor, 1971; Miller et al., 1991; Liu et al., 1996 and Chen et al., 1995). Due to the presence of the second stand, it was not possible to determine the full extent of this zone as the effects of that stand start to present themselves before the mixing was complete. The literature (Chen et al., 1995) cites that this zone should extend to about 18 tree heights downwind of the edge, but they note that surface roughness can have a large impact on how far this extends. The last zone was the point where the flow became re-equilibrated or a new equilibration of a different surface would be found. Again this would occur at about 18 tree heights (Chen et al., 1995), but can occur sooner with an increased surface roughness. Though the surface roughness downwind of the forest wasn't specified in detail, Li et al. (1990) noted in their simulation, that the streamlines seemed to become parallel again at about 8ht downwind of the forest edge. Wake and readjustment zones will be more closely examined in Chapter 5. 129 4.4. Conclusion Two types of porous canopies were examined. The first canopy type was made up of individual tree shaped elements that were given a uniform leaf area density. The other canopy, which was a more standard representation used in numerical modelling, was rectangular shape with multiple layers, each layer given a different leaf area density to reflect the canopy density of trees and ultimately forest they were modelling. For streamwise (horizontal) wind velocity, graphs, plots and statistics showed that the domain which had the porous tree shaped elements performed fairly well compared to wind tunnel observations and overall was comparable to the more traditional layered blockstand. This technique could be used in future numerical simulations for use with heterogenous stands where stand height, stem density and foliage density are not always uniform. With regards to turbulent kinetic energy prediction, values were generally over predicted compared to those calculated from the wind tunnel. Turbulent kinetic energy being over predicted was not surprising, as K-e turbulence models have been known to over predict these values, and no specialised turbulence closure scheme for the canopy was used. Three variations of the K-e turbulence model were used with both the porous tree elements and layered porous blockstand to determine which model performed closest to observations. Visually and statistically, the standard K-e turbulence model predicted horizontal wind velocity better than the realizable and Renormalized Group models, with the only exception being scatter (r2). Upon examination of the turbulent kinetic energy values, the standard K-e mode performed worse than the other schemes. Though it was far from ideal, the realizable K-e model had the 130 best performance out of the three K-S models at predicting TKE above the stands and clearing. The RNG K-E turbulence model performed adequately in predicting both horizontal wind speed and TKE values. It may have not been the best overall model, but it was never the worst model either. Examining the lower half of the profiles (up to two tree heights) in both the tree stand and blockstand domains, showed that the domains were adequate at predicting horizontal wind speed close to the ground, but consistently over predicted further up the profile (> 2 tree heights). Repeating the simulations using a domain twice (3 m) as tall as the original showed better agreement above the canopy, which indicated that the model may have been over predicting above the canopy in the standard height domain due to the conservation of mass. Better results with the normal domain were seen within the stands and above the clearing, these factors explain why the standard domain had slightly better statistics over the taller domain. The taller domain did fare better when it came to more closely predicting turbulent kinetic energy values. The better results were due to the reduced magnitude of the TKE output. Mixed results were seen with an upwind domain entirely absent of roughness elements and bluff bodies and composed entirely of the porous block stand. Statistics showed that the horizontal wind speed was better predicted with this simpler porous upwind section. The opposite was seen when looking at the values of turbulent kinetic energy. No difference was seen by applying a special boundary layer grid along the bottom of the domain. Differences may not have been picked up at the sampling locations as the lowest sampled points were 0.03 m off the ground, and the boundary layer grid was already expanded out to the normal grid size by that height. No evidence of improvement may have been found at all given that any slight improvement in boundary layer prediction may have been lost in the 131 overall large number of sampling points used in the statistics. The domain with the boundary layer style grid demonstrated that there was no improvement in model performance by increasing the time the fluid was simulated for. All simulations were four seconds long, with the exception of the RNG simulation (4.6 seconds). Most simulations had stabilised after about one second of simulation time. 132 1 2 3 4 5 6 7 8 9 Horizontal wind velocity (ms1) 10 11 12 13 Figure 4.1: Side view of the porous tree shaped stands and clearing showing isotachs (ms1) from the simulation using standard K-E turbulence model with an input velocity of 8.7 ms"1. Vertical lines indicate sampling locations (left to right) A to J (refer to figure 2.21 for more details). 0 iutti 1 2 3 4 5 6 7 8 9 Turbulent kinetic energy (TKE) (nrfV2) 10 Figure 4.2: Side view of porous tree stands and clearing, similar to figure 4.1 except showing isolines of turbulent kinetic energy (m2s2). The standard K-E turbulence model was also used for this simulaiton. i£ 0 Scatterplot of Horizontal Windspeed Observed vs. Modelled (Porous tree - std K-E at 4 sec) 10.0 | , . , . , . • -r- £3/ / idEr/'' '•i$F $f •J '•M-r** A A% J? • Regression line 1:1 line • 0.0 2.0 . 4.0 6.0 8.0 10.0 Measured (ms1) Figure 4.3: Scatterplot distribution of horizontal velocity from all sampling points in both the wind tunnel and porous tree domain simulated using the standard K-E turbulence model with an input velocity of 8.7 ms"1 135 f r - - - - V- Wind tunnel 3.2(1.225) P o r - t r e e i n e +1.255 PoMree-flNG line +1.255 Por-lree-real h i e +1.255 ,••' Por-tree-RNG line -0.5 / / • - - - ^ ./ P o r t r e e line +1.505 Por-tree-RNG line +1.505 - Por-tree-real line +1,505 ; ,S ' ." / PoMrce-RNG b>e +0.13 W i n d tunnel 0 . 8 7 ( + 0 1 3 > /"! 1 / /' 'ind tunnel 11.7 (+1.76) oMreeine+1.76 or-tiee-RNGIine +1.76 ,, Y, - Por-tree line. +3,13 Por-tree-RNG line. +2 Por-tree-real line + 2 / // Wirtdlunnel 4.2 ( . 0 . 6 2 5 ) P a - l i w l n . +0.625 Por-tree-RNG lne_+0.625 Figure 4.4: Vertical profiles of horizontal wind speed from both wind tunnel observations and output from porous tree elements simulation using the three variants of the K-E turbulence models. Locations a and b are within the first stand, upwind of the clearing; c, d and e are from the first half of the clearing; f, g, and h are locations from the second half of the clearing and i and j are locations found in the stand downwind of the clearing. Locations of each plot in the legend are from the upwind edge of the clearing (refer to figure 2.21). Windtmnel6.7(1.00S) Wmdtunnel-12.2(-1.B3) P o r t r e e i n e -1.83 Por-tree-RNG line -1.63 - Por-tree-real line -1.83 •• \ ') jsrSCL • 1 A x\ ' VA ... c Windtunnel-3.3(-0.5) P o M i e e B n e -0.5 PoMree-RNG Sne_-0. !fee-realline..+1.505 •e-RNGNne +0.1; idtunn e IO.87(- ) -0.13) j , Q 0 15 , - 0 45 E_ .. Je ' - .jEs^r- JJ W i n d unnel 2.5 +0.375) Pnr-t • s i n e +0.375 ee-RNG i n Por-t +0.375 ;. ' \ \ \ ; \ WmdtunneM4.2(+2.1 P o M . e e line +2.13 P o r - t w e - R N G line +2. W i l d tunnel 4 2 (+0.625) C.62S Por-tree-RNG lie +0.626 Figure 4.5: Vertical profiles of turbulence kinetic energy calculated from wind tunnel observations and output from the porous tree elements simulation using the standard K-t turbulence model and two variants. Letters in top left of each plot reference the location of that plot as demonstrated in figure 2.21. W i l d tunnel 6.7(1.005) - PoMteefcie +1.005 PoMiee-RNG fcie_+1.005 - Por-tiee-reallhe +1.005 .,» 0 75 r > 5.173 5.173 5.173 standard RNG realizable 5.862 5.666 5.556 P 2.116 2.116 2.116 ^0 2.716 2.657 2.548 P S 286 286 286 N -0.710 -0.748 -0.556 b 1.270 1.240 1.182 m 0.819 0.674 0.589 MAE 0.895 0.708 0.543 RMSE, 0.392 0.425 0.492 RMSEU 0.977 0.825 0.732 RMSE 0.958 0.970 0.975 d 0.979 0.974 0.963 r2 3.862 1.705 1.705 1.705 standard RNG realizable 3.485 5.033 0 K-e P 1.361 1.361 1.361 *0 1.066 0.922 0.626 P S 286 286 286 N 2.018 2.477 3.569 b 0.860 0.812 0.858 m 1.780 2.157 3.327 MAE 1.782 2.160 3.329 RMSES 0.919 0.769 1.250 RMSEU 2.005 2.293 3.556 RMSE 0.364 0.328 0.225 d 0.256 0.305 0.156 r2 _ Table 4.2: Statistics of turbulent kinetic energy for the porous tree stand simulation using three versions of the K-S turbulence model (variables defined oo in table 4.1). 0 K-6 Table 4.1: Statistics of horizontal velocity for the porous tree stand_simulation using three versions of the K-S turbulence model. In the table below, O is the averaged wind speed from all wind tunnel observations. P is the average wind speed output from sampled points from the simulation. The standard deviation for the observed data is noted by s0 and sp for the model output. N is the number of sampled points. The y-intercept is represented by b and m is the slope of the best fit line. MAE stands for Mean Absolute Error. RMSE is the Root Mean-Square Error with RMSES and RSMEU being the systematic and unsystematic components of the RSME, respectively. The Willmott d score, ranging from 0 to 1, is indicated by d, and finally, the correlation coefficient is r2. 1 2 3 4 5 6 7 8 9 Horizontal wind velocity (ms~1) 10 11 13 •^^^•••fflfflWM^BSilffiSSiiipilHiB 12 Figure 4.6: Side profile view of layered blockstands showing isotachs from the standard K-e turbulence model with an input velocity of 8.7 ms"1. Vertical grey lines represent the sampling locations (A-J) as defined in figure 2.21. Stem flow evidence can be seen in the second stand downwind of the forest edge near the ground. 0 o 1 2 3 4 5 6 7 8 9 2 2 Turbulent kinetic energy (TKE) (m" s ) 10 Figure 4.7: Simulation output from the standard K-E turbulence model. Side profile view of the layered canopy domain showing the stands and clearing. Isoiines of turbulent kinetic energy show areas of increased turbulence are found just above the tops of the canopy and close to the leeward end of the stands. TKE above the second stand is greater than what was simulated above the first stand. '.£• 0 1 2 3 4 5 6 7 8 9 Horizontal wind velocity (ms'1) 10 11 12 13 Figure 4.8: Side profile showing isotachs of horizontal velocity from the porous blockstand simulation much like the one in figure 4.6 but instead using the RNG K-£ turbulence model. 0 1 2 3 4 5 6 7 8 Turbulent kinetic energy (TKE) ( m V ) 10 Figure 4.9: The simulation from figure 4.8 except showing isolines of turbulent kinetic energy. TKE values over both stands seem to be fairly equal, with a small patch of TKE >6 m2s"2 over the second stand. - « ^ § K % « • ! , % .* 0 4^ 3 10 12 13 L.i-im JJ.JI... irmmwm-»T--p-l 11 i.--»-i"MB«»w»MW»p^M»p^l^B!^^p^Mpm.L.jjm.u|t.. 4 5 6 7 8 9 Horizontal wind velocity (ms1) •tpuJUi L ujum-in n.ii _ J J > . 1 2 i,. .i..-um^jJ^— Figure 4.10: Similar to figure 4.6, horizontal velocity isotachs from the simulation using the realizable K-£ model in the porous blockstand domain. ««MWI*.- 0 - v .•• : . > • 1 2 3 4 5 6 7 8 9 Turbulent kinetic energy (TKE) (m'V2) 10 Figure 4.11: Turbulent kinetic energy isolines from the layered domain simulated using the realizable K-E turbulence model and an input velocity of 8.7 ms"1. Greater TKE was predicted over the second stand. &••.-••• 0 Scatterplot of Horizontal Windspeed Scatterplot of Horizontal Windspeed Observed vs. Modelled (std K-£ at 4 sec) Observed vs. Modelled (std K-e at 6 sec) Regression line •1:1 line o.o k 2.0 Measured (ms1) 4.0 10.0 Measured (ms"'| Figure 4.12: Scatterplot distribution of horizontal velocity from observed and modelled sampling points for two unsteady simulations run for 4 seconds (a) and 6 seconds (b). Both simulations were identical in every other way using the porous blockstand canopies, the standard K-E turbulence model and an input velocity of 8.7 ms 1 . These plots show virtually no difference in output for the longer iterated simulation. Scatterplot of Horizontal Windspeed Observed vs. Modelled (RNG K-E at 4.6 sec) 10.0 w 6.0 Regression line 1:1 line 10.0 Measured (ms 1 ) Figure 4.13: Scatterplot distribution for observed wind tunnel values versus modelled output. Aside from the use of the RNG K-e turbulence model and the time of simulation (4.6 seconds) all other settings were identical to those using in figure 4.12. 145 Scatterplot of Horizontal W i n d s p e e d a Scatterplot of Horizontal W i n d s p e e d Observed vs. Modelled (realizable K-E at 4 sec) 0.0 2.0 4.0 6.0 8.0 b Observed vs. Modelled (realizable K-E at 4.6 sec) 10.0 Measured (ms"') 0.0 2.0 4.0 6.0 8.0 10.0 Measured (ms') Figure 4.14: Scatterplot distribution of horizontal velocity from observed wind tunnel data and output from the porous blockstand canopy domain using a realizable K-£ turbulence model and an input velocity of 8.7 ms"1. No difference is seen between the four (a) and 4.6 (b) second simulations. 146 ( 2 4 ''> . ." / 6 U(m/s) 8 I 1, • 10 / — - / ." / / ./ ' real line_+1.005 mg4.6line_*1.005 Wind tunnel 6.7(1.005) // -- U(m/s) Wind tun nel-12.2 (-1.83) std line -1 .83 - realfcifi.-1.B3 m g 4.6 line -1.83 f | 0.45 OSO 0.75 0 15 0.30 £ 0.45 0 75 0 " 2 -...__ • / - U (m/s) S//' • " / . W i n d tunna! 8.2 (1.225) std fins +1.255 r e a l i s e +1.255 . lng4.6line_.+1.255 .• U (m/s) ^ ' ' ' // W i n d tunnel-3.3 (-0.5) s t d l i i e -O.S / /,' / / / o' ' / ? 2 " • - ! 000 .' I£ ,.r £ '•r/t ; El - - 14S g j"' — i ___ " " i'S U (m/s) //'' Windtunnel10(+1.505) s t d line +1.505 real lire +1.505 m g 4 . 6 l i n e +1.505 ,• U(m/s) ^ / y WndtunnelO.B7(+0.13) / .' ----- 10 / ' 2 - . ~h; - / - "" *, - •"=• 0 60 n:- o 15 0 4, 0.75 - - / • U (m/s) , .---"< ' stdfine+1.76 realline +1.76 mg4.6liie.+1.76 U (m/s) mg 4.6 Ihe. +0.375 W i n d tunnel 2.5 (+0.375) stdfine+0.375 . „• / / 10 ' •' 1 - • - - - i; I J "o3„ >. 0 4S OSO 0 7S 0 00 0 45 0.0 0 f. >r " 2 -___ •t / 4 s U (m/s) Z/s' ' e "/' •"" W i n d tunnel 14.2 (+2.13) std line +2.13 realline +2.13 U (m/s) ;£'' realline. +0.625 m g 4 . 6 » n e +0.635 W i n d tunnel 4.2 (+0.625) / ' / i f, , 50 / / ' • 2 J -- Figure 4.15: Vertical profiles of horizontal wind speed from the wind tunnel observations and output from three versions of the K-E turbulence model. The three versions of the K-E turbulence models are: the standard (std) model, the realizable (real) model and the renormalized group (RNG) model. The ten plots represent the ten sampling locations (a-j) as shown in figure 2.21. Positions mentioned in the legend are from the upwind edge of the clearing. 0.30 * 0 45 0 60 0 75 •u 0 45 0.75 a d lne_+1.0O5 - rearSrw +1.005 Wmdtunnel-12.2(-1.83) — - std fine_-1.83 - - r e a l f n e -1.83 mg4.66ne_-1.83 f! Wind t u n n e l a . 2 (1.225) s t d l n e +1.255 reallinc_+1.2S5 m g 4 . 6 i n e +1.255 Windtunne!10(+1.505) - stdline +1.505 - real line +1.505 mg 4.6s l i r e +0.13 stdEne +0.13 • • - i -r ^ i t \ \ 1.76 m g 4 . S i ne +1.76 W i n d tun •el 11.7 (+ .76) - " W i n d t u n n e l 14.2 (+2.13) - s t d l i n e +2.13 I h e +2.13 Figure 4.16: Vertical profiles of turbulent kinetic energy. Wind tunnel values were calculated using wind statistics from the wind tunnel observations. Output from the porous blockstand domain includes: standard (std), realizable (real) and renormalized group (RNG) K-Z turbulence models. Each plot is identified by a letter which is also the location where the data and output were sampled from in the wind tunnel or domain, respectively (refer to figure 2.21). A :-- 5.173 5.173 5.173 standard RNG realizable 5.838 5.586 5.448 P 2.116 2.116 2.116 So 2.808 2.714 2.548 P S 286 286 286 N -0.960 -0.954 -0.633 b 1.314 1.264 1.176 m 0.852 0.679 0.568 MAE 0.941 0.695 0.463 RMSES 0.392 0.458 0.552 RMSEU 1.019 0.833 0.720 RMSE 0.958 0.970 0.976 d 0.980 0.971 0.953 r2 0 1.705 1.705 1.705 K-e standard RNG realizable 3.333 0.626 0.626 0.626 5.444 3.920 s0 P 1.046 0.952 1.46 P S 286 286 286 N 1.824 2.391 4.069 b 0.885 0.897 0.806 m 1.628 2.215 3.738 MAE 1.630 2.216 3.740 RMSES 0.887 0.768 1.374 RMSEU 1.856 2.345 3.985 RMSE 0.384 0.323 0.202 d 0.281 0.349 0.119 r2 Table 4.4: Statistics of turbulent kinetic energy for three K-e turbulence model schemes (standard, realizable and RNG) porous blockstand domain (variables defined in table 4.1). 0 K-e Table 4.3: Statistics of horizontal velocity for three K-e turbulence model schemes (standard, realizable and RNG) in the porous blockstand domain (variables defined in table 4.1). 1 2 3 4 5 6 7 8 9 Horizontal wind velocity (ms1) 10 11 12 13 T~$:--.A Figure 4.17: Filled horizontal isotachs from the simulation of the three dimensional porous blockstand domain using the standard K-E turbulence model. Input velocity was 8.7 ms"1. 0 1 2 3 4 5 6 7 8 9 2 2 Turbulent kinetic energy (TKE) (m" s ) 10 11 Figure 4.18: From the same simulation in figure 4.17, isolines of turbulent kinetic energy are displayed. A larger area of TKE was predicted in the downwind stand. 0 Scatterplot of Horizontal Wirtdspeed Qbse*veo vs Modelied (30 std K-Eat 4 sec) 10 C Regression lir.e; I 1.1'line :! Figure 4.19: Scatterplot of horizontal wind speed, observed results measured in the wind tunnel and predicted output from the three dimensional porous blockstand domain simulated using the K-£ turbulence model. 152 (J (m/s) u (mft) u (m/s) u 4«lne_t1JS6 - 3DLaywwf4slkw_-0.5 [ EOrfft") 6 10 e io WindtunmllO^LEOS) 2D Layered 4* lne_*1 JOS 3D Layered 4B In*_+1 JOB B WWtumelO.S7MI.13) 91) Layered 4i lire *0.la - 3D Layered 4* Una 44.13 j ' 1 Wrndtuniwl 11.7 (+1.76) 2DLiyarK!4*lna_+1.7fi Windtamel 2.5 04.375) 3DL*yer»d4alte_-i4.375 - 3DUyem)4Blhe_*a.37S ; 1 - i J . \ 6 10 /hdtumM4-2(+2.13) D Layamd 4a Hr»_f2.13 DLayamd4alM_+2.i3 Wind tunnel A2 (+0X351 3DLn*ensd 3D Blockstand • RNG K-E x All canopy (no upwind elements) A Realizable K-e • Boundary Layer Grid 4s (All canopy) Figure 4.24: Taylor diagram of all ten simulations done in this chapter and how they all relate to one another and to the "perfect simulation" or wind tunnel (Ref.)- 160 3.544 5.173 full half 3.544 5.173 full half 3.544 5.173 half full 0 3.721 5.862 3.546 5.666 3.524 5.556 P 1.600 2.116 1.600 2.116 1.600 2.116 So P 1.934 2.716 1.802 2.657 1.694 2.548 S 150 286 150 286 150 286 TV -0.463 -0.710 -0.344 -0.748 -0.091 -0.556 b 1.180 1.270 1.098 1.240 1.020 1.182 m 0.425 0.819 0.348 0.674 0.360 0.589 MAE 0.338 0.895 0.156 0.708 0.038 0.543 RMSES 0.426 0.392 0.415 0.425 0.462 0.492 RMSEU 0.544 0.977 0.444 0.825 0.463 0.732 RMSE 0.976 0.958 0.983 0.970 0.980 0.975 d 0.951 0.979 0.947 0.974 0.926 0.963 r2 real RNG std K-e half 1.550 1.705 full 1.705 full 1.550 1.550 half half 1.705 full 0 3.842 3.485 4.089 3.862 5.356 5.033 P 0.793 1.361 0.793 1.361 0.793 1.361 *o P 1.225 1.066 1.096 0.922 1.600 0.626 S 150 286 150 286 150 286 TV 2.279 2.018 2.684 2.477 3.837 3.569 b 1.008 0.860 0.906 0.812 0.979 0.858 m 2.292 1.780 2.539 2.157 3.805 3.327 MAE 2.292 1.782 2.540 2.160 3.805 3.329 RMSES 0.929 0.919 0.828 0.769 1.399 1.250 RMSEU 2.473 2.005 2.672 2.293 4.055 3.556 RMSE 0.399 0.364 0.376 0.328 0.268 0.225 d 0.425 0.256 0.430 0.305 0.235 0.156 r2 Table 4.12: Quantitative statistics of turbulent kinetic energy from the full and bottom half of the analysed profile for the porous tree stand simulation using three versions of the K-S turbulence model (variables defined in table 4.1). real RNG std K-e Table 4.11: Quantitative statistics of horizontal velocity from thefull and bottom half of the analysed profile for the porous tree stand simulation using three versions of the K-e turbulence model (variables defined in table 4.1). 3.544 5.173 full half 3.544 5.173 full half 3.544 5.173 half full 0 3.611 5.838 3.405 5.586 3.390 5.448 P 1.600 2.116 1.600 2.116 1.600 2.116 So 1.980 2.808 1.792 2.714 1.602 2.548 P S 150 286 150 286 150 286 N -0.681 -0.960 -0.456 -0.954 -0.011 -0.633 b 1.211 1.314 1.089 1.264 0.960 1.176 m 0.423 0.852 0.367 0.679 0.380 0.568 MAE 0.344 0.941 0.199 0.695 0.167 0.463 RMSES 0.421 0.392 0.425 0.458 0.463 0.552 RMSEU 0.543 1.019 0.469 0.833 0.492 0.720 RMSE 0.977 0.958 0.981 0.970 0.976 0.976 d 0.955 0.980 0.944 0.971 0.916 0.953 r2 real RNG std 1.705 full 1.550 1.550 3.710 3.333 4.142 3.920 1.705 full half 5.770 5.444 1.705 1.550 P O half full half K-e 0.793 0.626 0.793 0.626 0.793 0.626 *o 1.187 1.046 1.121 0.952 1.640 1.46 P S 150 286 150 286 150 286 N 2.085 1.824 2.602 2.391 4.369 4.069 b 1.048 0.885 0.993 0.897 0.904 0.806 m 2.160 1.628 2.592 2.215 4.220 3.738 MAE 2.160 1.630 2.160 2.216 4.220 3.740 RMSES 0.848 0.887 0.798 0.768 1.475 1.374 RMSEU 2.320 1.856 2.712 2.345 4.471 3.985 RMSE 0.419 0.384 0.374 0.323 0.246 0.202 d 0.490 0.281 0.493 0.349 0.191 0.119 r2 Table 4.14: Quantitative statistics of turbulent kinetic energy from the full and bottom half of the analysed profile for the porous blockstand simulation using three versions of the K-e turbulence model (variables defined in table 4.1). real RNG std K-e Table 4.13: Quantitative statistics of horizontal velocity from the full and bottom half of the analysed profile for the porous blockstand simulation using three versions of the K-S turbulence model (variables defined in table 4.1). 5. Effects of gap size on forest edge stress 5.1. Introduction Forest edges are widely known to be the main area where blow down and wind throw occur in forests. This is especially true in cut blocks and along riparian zones within the first three to five years after a harvest (Ruel et al., 2001). Much work in the field is focussed on forest edges and the interaction between the airflow and forest that occurs at this junction and immediately downwind of it (eg. Raynor, 1971; Irvine et al., 1997). This can be divided into two types: exit flows (flows from forests to clearings) and flows into forests. Exit flows have been pretty well studied, with the first major investigation by Raynor (1971). Much of the investigation involved examining how the flow adapted to the new surface. Exit flows and the progression of the air as it moves downwind of a forest edge has been divided into three zones (figure 5.1). The first zone is the quiet zone (Chen etal., 1995; Juddetal., 1996 and Lee, 2000), this zone is immediately downwind of the edge and extends for approximately 4-7 x/hx (Raynor, 1971 and Chen et al., 1995) (where hx is the stand height in the referenced article). The wind profile in this area is very similar to the profile seen before the forest edge. The second zone is the wake (Chen et al., 1995) or mixing (Lee, 2000) zone, as this is the area where the flow is adjusting to the new surface conditions (Judd et al., 1996 and Lee, 2000). This 163 zone has a large range of distances downwind of the forest edge where it can be found. In their simulations, Li et al. (1990), saw evidence from streamlines and pressure output that the edge effects from the stand only persisted for 10 x/hx downwind of the edge. Chen et al. (1995) believed that this zone extended to a distance of 18 x/hx, with their own wind tunnel readings still showing that the horizontal velocity at 22 x/hx was less than the potential value with no upstream obstacles. The final zone was the re-equilibration or readjustment (Lee, 2000 and Chen et al., 1995, respectively) zone. Lee (2000) described the zone as the area where the wind profile is fully adjusted to the new surface, while Chen et al. (1995) describes it as the region far downwind of the forest edge where the vertical diffusivity decreases to less than half the maximum value at heights below hx. Gash (1986) carried out field examinations of a forestheath interface and found by examining the collected wind data, that after 20 x/hx, measurements could not be considered representative of the new surface type. At 25 x/hx, Gash (1986) reported that 90% of the adjustment had been completed, with the presence of instabilities still being reported. For a flow to be completely adjusted, Gash (1986) commented that it may require a fetchof70x/h x . Flows into forests have also been described as distinctive zones, but not to the same level of detail, as flows out of forests. Most of the literature identifies three key zones: an approach zone, adjustment or edge zone and an equilibrium zone. An approach zone is upwind of the forest or canopy edge and is identified by a sudden deceleration of the air flow (Irvine et al., 1997 and Lee, 2000). An adjustment zone is reported to be about 10 x/hx long downwind of the forest edge and is where the flow must adapt to the new surface conditions (Yang et al., 2006 and Lee, 2000). It is the transition area between the edge and where the flow reaches equilibrium. Irvine et al. (1997) discuss the concept of an equilibrium layer that has completely adjusted to the new 164 surface. Flows into forests are where most wind damage can occur, especially within the first few years after the edge was formed (ie. harvesting). This is usually along the edges, where straight line wind blows down trees that are not accustomed to such wind loads (Yang et al., 2006). Conversely, trees a few tree heights downwind of long established, wind firm edges might be more vulnerable to damage than those trees at the edge. Yang et al. (2006) explains that trees along the edge have had to adapt to their new wind regime with higher mean wind speeds, but those downwind of the edge are not wind firm and may actually be more susceptible to extreme wind events. The purpose of this chapter is to examine the downwind edge of a forest clearing to determine if there is an optimal clearing size in which the wind stress on the trees is minimized to help prevent wind damage. 5.2. Setup Eleven domains were created for this study each identical to the blockstand domain used in the previous chapter, the only exception was the length of the domain as dictated by the size of the gaps or clearings. There were 11 clearing sizes used, the different sizes can be found in table 2.2. All simulations were performed using the standard K-S turbulence model. Three wind speeds were examined for this chapter. The speed most thoroughly investigated was the wind speed that was used for the other simulations, 8.7 ms"1; the other two speeds that were simulated were 4.35 ms"1 (half) and 17.4 ms"1 (double). Based on the stand size 165 and the operating conditions for simulations (p=\.225 kgm"3 and/i= 1.7894 x 10"5 kgm's" 1 ), these three wind speeds, 4.35 ras'1,8.7 ms-1 and 17.4 ms"1, have corresponding Reynolds numbers of approximately 44500, 89300 and 178700, respectively. 5.3. Results and discussion A large number of simulations were run for this study. Output of horizontal wind velocity and turbulent kinetic (TKE) energy will be presented along with more quantitative analysis specifically dealing with the downwind edge of the clearing. 5.3.1. Horizontal wind velocity Isotachs show the horizontal wind velocity for the various domains. With the smallest gap, there is almost no change in the wind speed above the clearing (figure 5.2 a) but a small increase in wind speed can be seen just below 1 z/ht. As gap size increases so do the effects on the flow. The 4 x/ht gap size in the simulation did not have any re-circulation (any negative horizontal velocity, or reverse flow, would have shown up as white) within the clearing, in fact no recirculating effects were output downwind of a stand, which agrees with simulations by Panferov and Sogachev (2008). At heights of 0.4 hx, the simulations by Panferov and Sogachev (2008) reported no reverse flow in gaps at 1.7 and 6.1 tree heights in length. Wind tunnel models of forest gaps from Stacey et al. (1994) (2.06, 3.8 and 6.67 x/hx) and Novak et al. (2001) (10.9 x/hx) (figure 3.5 a) did not report reverse flow either. These studies are contrary to Frank and Ruck (2008) and Miller et al. (1991), who both report a small re-circulation cells in their 166 simulations of a 1 x/hx and 2.5 x/hx wide gaps, respectively. For Miller et al. (1991), the reverse flow was predicted below 0.25 z/hx. Wind speed data from field measurements did not support this; although, smoke tracer observations did show evidence of an intermittent rotor (Miller et al., 1991). Most of the studies in the literature deal with small gaps sizes relative to the range of sizes dealt with in the present study (eg. Stacey et al. 1994, Miller et al., 1991 and Gardiner et al. 2005). The next few gap sizes (8 to 14 x/ht) (figure 5.2 b-e) do have classic gap characteristics as greater and greater wind speeds infiltrate further below canopy height. Flow traits in studies with similar gap sizes include Frank and Ruck (2008) (9 x/hx) and Stacey et al. (1994) (6.67 x/hx). As the clearing became larger (20 to 100 x/ht), the properties of the flow began to look like a forest-clearing relationship as greater and greater wind speeds penetrate down to the floor (figure 5.2 f-k). Similar patterns were seen in larger gaps studies including Dupont and Brunet (2008) (21.3 x/hx) and Panferov and Sogachev (2008) (20, 25, 30, 40, 55, 70 x/hx) as well as forestclearing studies by Liu et al. (1996) and Gash (1986). Other features of larger gap sizes include additional wind speed characteristics that occur just above canopy height, both upwind and downwind of the down stream edge. With larger and larger clearing sizes, the air entering the stand decelerates and some of this retarded air is forced upwards causing a slight reduction of wind speed above the canopy (Stacey et al., 1994). This effect becomes more apparent with increasing gap size (figure 5.2 a-k) as greater wind speeds are able to penetrate further down towards the ground. Still above the canopy, but just downwind of the clearing, the other feature noted is the dip in the isotachs caused by the accelerated wind, which has been observed many times in similar canopy edges (eg. Irvine et al, 1997; Stacey et al., 1994 and Gash 1986) and would be expected over any similar obstruction. 167 It is interesting to note that these same features appear regardless of wind speed. Li et al. (1990) noted that with their simulation of flow into a forest, there was an area of higher pressure upwind of the edge and lower pressure area about 4 x/hx downwind of the edge. The four second simulations of both the slower (4.35 m s"1) (figure 5.3 a-e) and faster (17.4 m s"1) (figure 5.4 a-e) simulation speeds are fairly similar to the 8.7 m s"1 used for most of the simulations that were run in chapters 3 and 4, as well as the wind tunnel from figure 3.5 a. Isotachs of the three wind speeds, at all gap sizes, show very similar patterns, although at a different range of speed. The lower speed simulation is the only one that had wind speed return to the upstream speed downwind of the first stand. This occurred sometime after the 60 x/ht (figure 5.3 d) gap but before the 100 x/ht gap (figure 5.3 e), at around 70 x/ht. This would concur with a similar finding by Gash (1986) who found that for a flow to be in equilibrium with the surface a fetch might be required to be at least 70 x/hx. Though in the observations obtained by Gash (1986), at 70 x/hx downwind of the forest-clearing edge, wind speed was still increasing below 1 z/ht. From the faster wind speeds (8.7 ms-1 and 17.4 ms"'), figures 5.1 i-k and 5.4 d and e show that even between 60 x/ht and 100 x/ht, the wind speed has not returned to domain input values. 5.3.2. Turbulent kinetic energy Filled isolines were used to illustrate the distribution of TKE in the simulations. These simulations yielded some interesting insights about TKE at the tops of forest stands where gaps exist. The smaller gap sizes modelled appeared to have some of the largest areas of maximum TKE predicted over the downwind stand. This occurred at all the wind speeds simulated. In the domains with horizontal wind speed set at 8.7 ms"1, the gap sizes that were <60 x/ht, had larger 168 values of TKE over the second stand versus the first stand (figures 5.5 a-h). However, after about 60 x/ht, as the size of gaps increased, the TKE above the downwind stand no longer seemed to be greater (figures 5.5 i-k and 5.7 e). Figures 5.5 i-k show that gap sizes >60 x/ht will have TKE values roughly equal to (figure 5.5 i, j) or even less than (figure 5.5 k) the first stand. This was due to the TKE from the first stand being advected downwind (Dupont and Brunet, 2008). This is analogous to findings by Yang et al. (2006) who reported higher wind moments and thus a greater potential for wind damage were reported downwind of clearings and not along established forest edges. The 4.35 m s"1 (figure 5.6) and 17.4 m s"1 (figure 5.7) simulations shows that an areas average annual (endemic) wind speed will greatly affect how far downwind TKE will be detected. As one may expect, slower wind speeds have turbulence resolved downwind in a shorter length compared to the 8.7 m s"1 and 17.4 m s"' simulations. 5.3.3. Edge stress on downwind stand Stress (F) on the upwind edge of the second stand is determined for each domain (Eq. 2.45). The stress values calculated were then normalised (FRef) using air density and the input velocity from their respective simulations: F Rzf =p{URSff • (5-D These were then plotted to show how stress changed with increasing gap size (figure 5.8). The rate of increase along that edge seems to be logarithmically related to the clearing size. After 60 x/ht, the rate at which stress increases along the edge, notably decreases. With the other two wind speeds used, similar increases were seen in the stress values. The 169 slower model, had nearly identical stress values between the 60 x/ht gap and 100 x/ht gap. This should probably be expected as it was at the 60 x/ht gap size that the input wind speed had almost completely infiltrated towards the ground (figure 5.3 e). The faster wind speed had a similar rate of change; although, it did seem that the rate of increase dropped off quite a bit after about the 30 x/ht gap. Another plot was also made but showing gap transit time instead of gap distance. Gap transit time (in seconds) is the time it takes for wind to travel from the upwind edge to the downwind edge, based upon the input velocity (figure 5.9). Gap transit time is essentially defined as: _ . . gap size Gap transit time = M (5j) Re/ where gap size is the distance between two stands in metres and uRef is the input velocity of the wind tunnel or simulation or the average wind speed over the upwind forest stand. All points regardless of input velocity seem to cluster together along a single logarithmic curve. 5.3.3.1. Logarithmic fits Logarithmic lines of best fit were added to each of the plots (figures 5.8, 5.9). The standard speed (8.7 ms 1 ) graph (figure 5.8) had a logarithmic best fit line that appeared to generally fit the points from the various gap sizes. The lower speed domain (4.35 ms"1) did not fit the logarithmic line as well. This was mainly due to the last two domains (60 x/ht and 100 x/ht) having nearly identical momentum values. The higher speed domain fit the logarithmic profile fairly well, with two points slightly deviating from the line. With gap transit time (figure 5.9), the logarithmic fits were more closely grouped together 170 compared with the fits from figure 5.8. By collapsing all three lines into one curve a logarithmic profile can be plotted (figure 5.10). From this the relative or magnitude difference of total stress on the downwind edge of a clearing can be roughly predicted as: F 1.3827+ 1.0275 log,0 *Re/ f gap size .-_ . (m) /...\\ (5.3) V M Re/ J where Fn is the normalised stress (dimensionless), or magnitude difference (increase or decrease) from a reference stress value, ¥Ref. It should be noted that F„ would be suitable for any gap size or reference velocity in an identical model forest clearing. Also, the logarithmic curve in figure 5.10 should be representative of other forest clearings regardless of their stand composition or normalisation procedure. 5.3.4. Horizontal velocity profile Modelled horizontal velocity was plotted between stands in the 100 x/ht clearing for each input velocity at 1/3 z/ht. Both higher speed input velocities (8.7 and 17.4 ms"') were able to regain 80% of their speeds before the second stand (figure 5.11), while the slower speed simulation (4.35 ms"1) regains over 90% of the input velocity. Unresolved is the sudden apparent change in the lower velocity profile that begins at roughly 7 ht downwind of the upwind edge. This point can also be identified in figure 5.3 e where the 4 ms"1 isotach appears to have a different angle compared to the 3 ms"1 isotach and even other isotachs within the clearing of the other two wind speeds (figure 5.2 k, 5.4 e). When examined as a function of gap transit time, this feature can still be seen (figure 5.12). 171 5.4. Conclusion It was hoped that the stress/momentum values would eventually level out and give forest managers an idea, given the endemic wind speed in an area and site conditions, what the best gap size would be for a forest cut block. The gap sizes used, ended up being quite large, much larger than what other studies looked at, but based on the output, could have been even larger as stress values never completely levelled off as expected. Only the lower wind speed model (4.35 ms"1) showed a levelling off of momentum values on the edge of the second (downwind) stand. In terms of planning cut block size, there was a marked difference depending on which measure is most important, in determining cut block size. If protection from stress due to horizontal wind is a higher priority, then smaller gap sizes will offer greater protection as larger gap sizes allow higher wind speeds to infiltrate closer to the ground. If reducing turbulence downwind of an edge is more important, than opting for larger clearings would be more ideal. 172 Figure 5.1: Exit flow across a two-dimensional forest edge. The approach flow (A) upwind of the forest edge and at the forest edge (A1). Downwind of the stand, the flow field in the open is divided into quiet zone (D), mixing zone, (E) and re-equilibration zone (F) (after Lee, 2000). 173 1 4 5 6 7 8 9 Horizontal wind velocity (ms"1) |M!WpWWPyiUll:§,:, 3 10 11 12 13 ^^^^^I^^^WIS^^^II5J Figure 5.2 a: Horizontal velocity isotachs of airflow over two porous blockstands separated by a 4 x/h, clearing and simulated with an input velocity of 8.7 ms"1 and using the standard K-E turbulence model. 0 -J 1 2 3 4 5 6 7 8 9 Horizontal wind velocity (ms"1) Figure 5.2 b: Horizontal isotachs through a canopy with an 8 x/ht clearing. (UpmwwM mgin 0 10 11 12 13 B -J 0\ 4 6 7 8 9 Horizontal wind velocity (ms") 10 11 12 Figure 5.2 c: Clearing of 10 x/h, showing an increase in horizontal wind speed penetration towards the ground. 3 13 1*2 4 5 6 7 8 9 Horizontal wind velocity (ms1) 10 11 12 13 Figure 5.2 d: Horizontal velocity isotachs showing a 10.9 x/ht clearing, which at full scale (1:200 ht)would be just over 10 ha. 0 D 00 1 2 3 4 5 6 7 8 9 Horizontal wind velocity (ms~1) Figure 5.2 e: Forest gap of 14 x/h„ 3 m/s isotach is well entrained along the ground 0 10 11 12 13 Figure 5.2 f: Clearing representing a 20 tree height gap. 3 UJIIrVWUAPm, 4 5 6 7 8 9 Horizontal wind velocity (ms"1) 10 11 12 13 ^Mpipa^UUJUrL-IIUIkr o 00 10 11 12 13 Figure 5.2 g and h: Increased clearing sizes clearly shows the increased infiltration between the 30 h, (g) and 50 ht (h) gap size. In g, the 6 ms isotach is still above the height of the forest, in h, the isotach has moved down to just above the surface. 4 5 6 7 8 9 Horizontal wind velocity (ms"1) H 00 1 2 3 4 5 6 7 8 9 Horizontal wind velocity (ms~1) 10 11 12 13 K 181 Figure 5.2 i, j and k: Sixty, 80 and 100 ht. I, j and k all show the progression of the higher wind speeds infiltrating further down towards the surface. 0 Figure 5.3 a: Horizontal isotachs flowing over the 4 x/ht gap size from the 4.35 ms 1 (input velocity) simulation also using the porous blockstand canopy and standard K-S turbulence model. 2 3 4 Horizontal wind velocity (ms~1) 00 *• •"• * •'. • : .«..• .¥••;* . * - s Figure 5.3 b: Clearing size of 10 ht. i*: A ri'SK,", m ImT * ^ ;&H ^y® w#. l#:^ m-?s &£.&}•• .*••- • J3»"' **i # # ^ * 4 "*. ;-;£-! -HI 2 3 4 Horizontal wind velocity (ms1) B ,V**..- < , . . 3 1 ._;;.. ... 4^ K- Figure 5.3 c: Clearing size of 30 ht. S- 0 184 2 3 4 Horizontal wind velocity (ms1) 6 2 3 4 Horizontal wind velocity (ms1) D Figure 5.3 d and e: Top image (d) shows a gap length of 60 h„ while the bottom (e) is 100 ht. It appears that the flow has fully adjusted to the new surface in figure e. •>4F2L*&u 0 2 jfi^ffp 4 6 mmmmsmMmmf. IBWapMBffliiWS^WIB^WI^BiJil,! 4ULU jjuwjufc-jij^auiiw-a 8 10 12 14 16 18 20 22 24 26 Horizontal wind velocity (ms1) Itui, Figure 5.4 a: Isotachs of a cross section with a 4 ht gap in a simulation with a 17.4 ms"1 input horizontal velocity using the standard K-S turbulence model and the porous blockstand canopy. '^^^MMPvpfP^op^ 0 00 Figure 5.4 b: Isotachs from the 10 h, gap domain. •WiW-u^iyptfMijjj***- feRfr- •WJjuai.uujjuyjiiyiWMP. 8 10 12 14 16 18 Horizontal wind velocity (ms1) 20 22 26 i4ju.LiLjBij,ijpn 24 ,«|_JlJiUll#lllJUJ».. B ^ Figure 5.4 c: A 30 ht clearing with isotachs. 0 «i 6 8 10 12 14 16 18 Horizontal wind velocity (ms~1) 20 22 24 26 00 8 10 12 14 16 18 Horizontal wind velocity (ms~1) 20 22 24 26 Figure 5.4 d and e: top image (d) shows isotachs of the stands and 60 h, clearing; bottom image (e) shows a 100 ht clearing. 6 D o 1 2 3 4 5 6 7 8 9 Turbulent kinetic energy (TKE) (m'V2) 10 •..f'--- Figure 5.5 a: Filled isolines showing turbulent kinetic energy from the domain with a 4 h, gap between two porous blockstands and an input velocity of 8.7 ms"1, in a simulation using the standard K-£ turbulence model. Higher values of TKE are predicted over the second forest stand, downwind of the clearing. 0 \D 1 2 3 4 5 6 7 8 9 Turbulent kinetic energy (TKE) (m"2s2) Figure 5.5 b: Isolines of TKE in the domain with the 8 h, clearing. 0 10 1 2 3 4 5 6 7 8 9 Turbulent kinetic energy (TKE) (m'V2) Figure 5.5 c: A 10 ht clearing and the resulting TKE isolines in and above the clearing and stands. 0 10 2 3 4 5 6 7 8 9 Turbulent kinetic energy (TKE) (m"2s2) Figure 5.5 d: Filled TKE isolines in the domain with the 10.9 ht clearing. 1 10 4^ 2 3 4 5 6 7 8 9 Turbulent kinetic energy (TKE) ( m V ) 194 Figure 5.5 e: A clearing of 14 ht and the resulting filled isolines of TKE. 1 10 2 3 4 5 6 7 8 9 Turbulent kinetic energy (TKE) (m"2s2) Figure 5.5 f: Clearing of 20 h, with the second stand just out of range of this figure, showing TKE. :'W?'.""! 1 10 1 2 3 4 5 6 7 8 9 Turbulent kinetic energy (TKE) ( m V ) Figure 5.5 g and h: TKE shown in the 30 ht (g) and 50 ht (h) clearings. 0 10 -J VO - - — - — - — — - — — . . ~ — ~ — m ^ - M — j a 10 Figure 5.5 i, j and k: Domains, with gap sizes of 60 h, (i), 80 ht (j) and 100 ht (k), showing TKE. Figures i and j shows that the maximum TKE value over both stands is similar. Figure k shows that the peak TKE value above the second (downwind) stand is lower than what was predicted in the first (upwind) stand. - 1 2 3 4 5 6 7 8 9 2 2 Turbulent kinetic energy (TKE) (m" s" ) 00 2.5 3.0 Figure 5.6 a: Filled isolines of TKE from the simulations with the 4.35 ms'1 input velocity, this figure showing the domain with the 4 h, clearing using the porous blockstands and simulated using the standard K-S turbulence model. Unlike the other simulation with an 8.7 ms'1 input velocity, these two stands have similar maximum values over upwind and downwind stands. 1.0 1.5 2.0 Turbulent kinetic energy (TKE) (m"V2) 0.5 1.5 2.0 2 Turbulent kinetic energy (TKE) (m'V ) 1.0 Figure 5.6 b: TKE in and above two stands and a clearing 10 ht in length. 0 2.5 3.0 K> Figure 5.6 c: Filled isolines showing the distribution of TKE in and adjacent to a 30 ht clearing. o 2.0 1.5 1.0 Turbulent kinetic energy (TKE) (m" s ) 2.5 3.0 0.5 1.0 1.5 2.0 Turbulent kinetic energy (TKE) (m"V2) 2.5 3.0 Figure 5.6 d and e: TKE drops below 0.5 m2s"2 at around 40 ht in both the 60 ht (d) domain and 100 ht (e) domain. Even at these large gap sizes, the upwind and downwind stands have similar peak TKE values. 0 to o 3 6 9 12 15 18 21 24 27 30 Turbulent kinetic energy (TKE) (m"V2) 33 36 39 Figure 5.7 a: Filled isolines of turbulent kinetic energy from the porous blockstand canopy domain with the 4 ht gap. This simulation had the inlet velocity set at 17.4 ms"1 and was simulated using the standard K-8 turbulence model. The peak TKE seen and the extent of it over the second stand is noticeably higher than the first stand. Quite high values of TKE are seen between stands above the clearing. 0 A o to 3 6 9 12 15 18 21 24 27 30 Turbulent kinetic energy (TKE) (m'V2) Figure 5.7 b: TKE from the domain with the 10 h, gap. 0 33 36 39 •.TSrGWW B 3 6 9 12 15 18 21 24 27 30 Turbulent kinetic energy (TKE) (m"V2) g Figure 5.7 c: Gap from the 30 ht domain with isolines of TKE. 0 33 36 39 33 36 39 Figure 5.7 d and e: Turbulent kinetic energy is still present far downwind of the first stand in both 60 ht (d) and 100 ht (e) clearings. > K W » xrnm •*»:*• aan-JtsgiPa 9 12 15 18 21 24 27 30 Turbulent kinetic energy (TKE) (m"2s2) to ON o 1.4 0 . r . . . . 1- . / . 10 i // . —i- . . . ^/ »— . •'•' ' -, ' i . . • •—• 20 | -. . 30 1 _ jr'" | —r . i . -r- J . 1 . 40 L , - L. - r - , . •_ • , . '-• • 1 1 60 • 1 1 I 1 J 70 1 •-, Gap Size (ht) 50 •_ I ~"8l . •-1 1 1— | 1 — . —1 . •] •. . """ " . . T • ' " T ^ * ^ - 1 '- • ' ____ „---' 1 ' ' 1 i 1 1 80 1 1 l_ 90 1 —L 1 - ^ . 1, _l 100 -I L. J , 110 •_.! . 1 120 I • • -- « ^ Noramalised Stress (where uRB, = 4.35 ms"1) • Normalised Stress (where uRel = 8.7 ms'1) ""*, Noramalised Stress (where uM = 17.4 ms"1) N .- 1 Figure 5.8: Normalised stress values on the downwind edge of various gap sizes for three input velocities. 0.2 0.0 0.4 o 0.6 E CD J 0.8 -a m 2> 1.0 CO CO uT u_ 1 1.6 I- 1.8 2.0 Normalised stress for various gap sizes from simulations of three different input velocities / / ?1 1 /•/ / / I . I . I I / Wkf f 1 1 I 7 mm/ s . J * . x 1 yd -J 1 —T 1 1 1 S* *• "M ..*•'" 1 ' 1 '—— 1 I , I . p T i . • i 1 *~ i 1 • ' i 1 -i ] J 1 "| 1 1 -«-T 1 1 1 — ^ ' 1 ' ""•^ Noramalised Stress (where uRef= 4.35 ms"1) • Normalised Stress (where uRef= 8.7 ms'1) " j u Noramalised Stress (where uRef= 17.4 ms"1) 1 • - I Gap Transit Time (s) 0.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 4.0 0.2 0.4 1 ' '' / * 0.6 V- ; I i 1 /A / i 1 Figure 5.9: Normalised stress from each input velocity versus the time taken to cover each gap distance. E •S CD .2 0.8 CD T3 CO 5 i.o 1-2 1.4 1 i O CO CO LL LL LL II 1 |- —. 1 1.6 - 1.8 2.0 Normalised stress for various gap sizes from simulations of three different input velocities i to ooo ii "~% / /• / •/• • I• • / • • / • ^ — • • Gap Transit Time (s) 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 4.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 Normalised stress from various gap sizes and input velocities as a function of transit time Figure 5.10: Normalised stress from each clearing size and each input velocity simulated as a function of time required to flow from start to end of clearing. Normalised Stress to o - 0.8 0.7 0.6 0 10 20 30 ^^^^^^^^^aa^^__i__i__l__i__^_^_i__l_ - ' 0.1 0.0 H 0.2 0.3 o.4 H 40 60 Gap Distance (x/ht) 50 70 80 90 100 Input velocity = 4.35 ms"1 • Input velocity = 8.7 ms"1 , • Input velocity? 17.4ms'1 Figure 5.11: Horizontal velocity profile of the 100 x/h t clearings at 1/3 stand height between stands. £ > o CD k E 'o CD - - 0.9 m £> 0.5 •o CD - 1.0 o - 0.70 0.60 H 0.00 0.0 o.io H 0.20 0.30 0.40 H - 0.80 0.50 - 0.90 0.5 —r1.0 2.0 Gap Transit Time (s) 1.5 2.5 i • 3.0 1 ' 3.5 Input velocity = 4.35 ms"1 • Input velocity = 8.7 ms"1 • Input velocity = 17.4 ms"1 Figure 5.12: Similar to figure 5.11, but with horizontal velocity profile plotted against gap transit time. CD E & "Q .= CO - 1.00 6. Summary of conclusions Two different aspects of forest research were examined. The first involved the use of a commercial CFD program to numerically simulate air flow through and above forests and forest clearings. Three different aspects of CFD simulation were investigated: predicting turbulent kinetic energy using three variations of the K-S turbulence model, using tree shaped entities as the canopy versus using simplified blockstands and how various gap sizes affected airflow. 6.1. Comparison and testing of K-E turbulence models The third chapter of the thesis used solid tree shaped objects with three K-S turbulence models. Three variations of the K-S turbulence model were compared and tested. The standard K-S model was considered very robust and has been use by many other researchers for about two decades. The other two K-S turbulence models, Renormalized Group and Realizable, were ones that had been adapted by determining values of some variables that were constants in the standard version. The purpose was not to refine or enhance an existing turbulence model, but to see the differences between three similar turbulence models that were offered by the program and suggest if one may be better suited for use in airflow studies involving forest landscapes. The use of tree shaped objects was more of a test to determine if simulations using the K-S turbulence model could be successful despite having many relatively complex structures in the 211 numerical domain, and report realistic wind and turbulent values. This was followed up in chapter four which dealt with a porous canopy structure, one that would mimic real trees. The results showed that FLUENT, a computational fluid dynamics program, was able to successfully simulate to convergence tree shaped entities using the K-e turbulence models. Horizontal velocity results were reasonable given the solid structure of the canopy. From these results, analysis appeared to indicate that the standard K-S model yielded better horizontal velocity output than the other two K-E models. 6.2. Comparison and testing of different canopies and domains The fourth chapter dealt with a wide range of issues with many simulations being performed. Carrying over from the third chapter, tree shaped elements had porous settings applied, which were based upon LAD and Cd. Block-like, layered stands with porous settings were also used to examine the three versions of the K-S turbulence model. For both the porous tree and blockstand simulations, horizontal wind velocity was better predicted by the standard K-e turbulence model while the realizable K-S model was better at predicting the turbulent kinetic energy generated by the stands. Between both canopy types, and using the same standard K-S turbulence model, had nearly the same output for horizontal wind velocity. Values of turbulent kinetic energy from the tree shaped elements model were actually superior to those from the layered canopy. However, the superior turbulence model for predicting TKE was the realizable K-S model, which had better statistics with the blockstand domain. Though it should be pointed out that the overall TKE values predicted in these simulations did not compare well with those 212 calculated from wind tunnel statistics. Between the two domain setups, porous trees and blockstands, given a homogeneous stand, there does not seem to be any clear advantage for simulating a domain setup with tree elements over a blockstand. For horizontal velocity, the domain created with twice the height was determined to only perform better for areas well above forested locations. For areas within and just above a stand, as well as within forest clearings the standard domain height performed best. The largest issue was the reduced air speed that was seen throughout the vertical profiles at many locations. The 3 m tall domain regularly under predicted wind conditions leading to a better statistical result for the standard domain. Conversely, the lower TKE values in the taller domain, yielded better TKE values. Throughout all of the simulations, values of turbulent kinetic energy were in poor agreement with wind tunnel values. This was not surprising as it has been known that standard K-£ turbulence models do not predict TKE as well as many would hope. This is why many researchers in this field attempt to tweak models to better account for canopy properties, with some researchers moving to the use of large eddy simulations, which is the next step up on the numerical modelling ladder. When examining vertical profiles of horizontal velocity, it was noted that points higher up above the canopy did not compare well with the data collected in the wind tunnel. Simulations had more accurate results when output under 2 z/ht were examined. Output from models above 2 z/ht tended to over predict horizontal wind speed. By examining this shorter vertical output, the simulation with the porous trees was shown to have better output than the porous blockstand simulation. As with the full profile heights, there was no real noteworthy differences between porous tree elements and blockstands for homogeneous stands. One interesting note, for both 213 canopy types, the RNG K-e turbulence model generally preformed better than the standard K-e model for airflow up to two tree heights. For TKE values, the realizable K-S model performed best out of the three turbulence models, with a slightly better result in the blockstand canopy versus the tree element canopy. For the other domains examined including the three dimensional one, none seemed to improve markedly from the blockstand. Usually slight improvements were seen with respect to horizontal velocity, only to be offset by a worse TKE performance. Similar to how the porous tree elements that would likely see better output in a heterogeneous stand, the three dimensional domain would probably be worth the effort if more dynamic edges or terrain were required. 6.3. Effects of gap size on forest edge stress Horizontal wind velocity and turbulent kinetic energy were examined in and above various sizes of forest gaps. Stress/pressure on the downwind edge of the clearing was also examined. Exposure of progressively higher values of horizontal wind velocity was seen with increasing gap size as the higher speed airflow was able to infiltrate closer to the ground as distance increased. Conversely, downwind turbulent kinetic energy decreased the larger the gap size became. This may have been due to the larger gap sizes allowing for the dissipation of some of the turbulence before the next stand. From all of the gap sizes examined, and from three different input velocities, a logarithmic equation was generated that would roughly predict stress along the downwind edge of a forest clearing, given a stand with similar characteristics. 214 6.4. Reflections and recommendations for future work The two largest contributions that this thesis has made will likely be with numerical model canopy design and the size of gaps being examined. The content of my thesis will be of interest to other researchers and is not likely to have practical use for the professional forester on the ground. Throughout my thesis, I have not come across any other numerical model that has tackled canopy design in a large stand in the same manner as I have done: by using a stand that has been made up of trees resembling real trees. My thesis was seen as an initial step into additional research with more complex forests. The two main types of canopies created and examined in my thesis were the two dimensional porous tree elements stand and the porous blockstand canopy. The domain with the porous blockstand canopy was much easier to construct and simulate than the domain with the porous tree elements, though the two canopy styles validated similarly. Given the amount of time required to construct the porous tree elements and the difficulty in obtaining a successful simulation, I would suggest the porous blockstand style canopy for most uses in a CFD simulation. While porous tree elements would have their benefits in some situations (eg. edge effects downwind of a clearing or stand reserves in a clearing), a porous blockstand upwind forest would probably perform adequately and be easier to implement. Also, if a finer degree of detail is required close to a stand or individual trees are being investigated a porous tree element would likely yield better results. With respect to three dimensional domains versus the two dimensional domain, not much additional construction time is required for the porous blockstand. Simulation time, while longer (depending on the width of the domain) did not appear to cause any great difficulties with convergence. A three 215 dimensional domain with porous tree elements, though not introduced in this thesis, was much more difficult to construct over the two dimensional domain. This complex domain also posed problems with meshing errors that were difficult to find and repair. Successful simulations using the K-S turbulence model were not achieved using a three dimensional porous tree element stand. One issue I would like to understand more would be the reason why turbulent kinetic energy was predicted to be so dissimilar from the wind tunnel. While it is known that K-S turbulence models are not the best for predicting TKE, I was surprised to see them so far off in my simulations and why other studies had more success with them. The next logical step would be to further examine forest edges and see how altering the size, shape and LAD of the trees along the edge might compare with wind tunnel and field work already completed. This work should be implemented using large eddy simulation and therefore three dimensions, as numerical simulations on this scale seem to be progressing more and more to that style of simulation. 216 7. Nomenclature Af = frontal area, m 2 As = surface area, m 2 b = y-intercept Clc, C2c, C3e, C/( = K-S model constants (model dependant) C2 = inertial resistance factor, m"1 Cd - Canopy drag coefficient F = Momentum (Stress) on downwind stand, kgnv's"2 FRef = Stress based upon fluid density and upwind or reference velocity, kgm's Fn = Normalised stress, F/FRef p = external body force ht = height of stand in study, 0.15 m hx = height of stand for (other) referenced stand Hz = timesteps (or samples) per second / = unit tensor I] = first interval on a mesh edge in domain, m LAD = leaf area density, m 2 m"3 (m"1) /„ = last interval on a mesh edge in domain, m m = slope n - number of intervals/nodes along a mesh edge in domain 217 = observed (measured) value = observed mean = static pressure = predicted (model output) value = predicted mean = ratio of mesh node spacing in Gambit = Momentum sink, kgm'Y 2 (Nm~3) = observed standard deviation = predicted standard deviation = time, s = velocity, ms"1 = fluctuating velocity component, ms"1 = mean velocity component, ms"' = mean fluctuating velocity (standard deviation), ms"1 = upwind or reference velocity, ms"1 = streamwise coordinate, m = spanwise coordinate in three dimensional domain, m = vertical coordinate, m = streamwise coordinate in the three dimensional domain, m = vertical coordinate in the three dimensional domain, m = inverse Prandtl number = thermal diffusivity, mV 1 = Kronecker delta 218 £ = dissipation of kinetic energy, m V 3 K - turbulent kinetic energy, mV 2 jj. = dynamic viscosity (air), 1.7894 x 10"5 kgm'V (Pas) jut = turbulent viscosity, kgnv's"1 v = kinematic (eddy) viscosity, mY 1 p = fluid density (air), 1.225 kgm"3 pg = gravitation force aK, at. - Prandtl numbers for turbulent kinetic energy and its dissipation, respectively q> - scalar variable (p' - fluctuating scalar component (p = mean scalar component j - stress tensor X - eddy diffusivity, mY 1 co = angular velocity £X. = mean rotational rate 219 8. Bibliography Amiro, B. D. 1990. Drag coefficients and turbulence spectra within three boreal forest canopies. Boundary-Layer Meteorology 52(3): 227-246. doi: 10.1007/BF00122088 Anderson, Jr., J. D. 1996. Basic philosophy of CFD. In Wendt, J. F. (Ed.), Computational fluid dynamics: An introduction. Springer: New York. p. 7-8. ANSYS, Inc. 2006. FLUENT® Releases 5.0-6.3.26, Help System, User Guide, ANSYS, Inc. Canonsburg, Pennsylvania, USA. Atterson, J.A. 1980. Gambling with gales. Paper presented to the British Association for the Advancement of Science, Section K (Forestry), Salford, UK. p. 10. Belcher, S.E., Jerram, N. and Hunt, J.C.R. 2003. Adjustment of a turbulent boundary layer to a canopy of roughness elements. Journal of Fluid Mechanics. 488: 369-398. doi: 10.1017/S0022112003005019 British Columiba Ministry of Forests. 1992. Annual report of the ministry of forests 19901991. Ministry of Forests, Victoria, British Columbia. British Columbia, Ministry of Forests. 1995. British Columbia Forest Practices Code. Ministry of Forests, Victoria, British Columbia. British Columbia, Ministry of Forests and Range. 2006. The state of British Columbia's forests, 2006. Ministry of Forests and Range, Victoria, British Columbia. CCFM (Canadian Council of Forest Ministers and Natural Resources Canada). 2006. Criteria and indicators of sustainable forest management in Canada: national status 2005. Canadian Council of Forest Ministers: Ottawa, pp.154. Chen, J.M., Black, T.A., and Novak, M.D. and Adams, R.S. 1995. A wind tunnel study of turbulent airflow in forest clearcuts. In Wind and Trees. Coutts, M.P. and Grace, J. (editors). Cambridge University Press, New York, USA. p.71-87. Clark, T.L., Mitchell, S.J., and Novak, M.D. 2007. Three-dimensional simulations and wind-tunnel experiments on airflow over isolated forest stands. Boundary-Layer Meteorology. 125: 487-503. doi: 10.1007/sl0546-007-9198-l Counihan, J. 1969. An improved method of simulating an atmospheric boundary layer in a wind tunnel. Atmospheric Environment. 3:197-214. 220 Doty, K., Tesche, T.W., McNally, D.E., Timin, B. and Mueller, S.F. 2002. Meteorological modeling for the Southern Appalachian Mountains Initiative (SAMI). Tennessee Valley Authority, Knoxville (TN). Available from: http://www.tva.gov/sami/met/eval/SAMI_v2A.pdf [last accessed 27 February, 2009; posted 14 August, 2002] Dupont, S. and Brunet, Y. 2008. Impact of forest edge shape stability: a large-eddy simulation study. Forestry. 81(3): 299-415. Ferziger, J.H. and Peric M. 1999. Computational Methods for Fluid Dynamics. New York: Springer, p. 268-288. Foudhil, H., Brunet, Y., and Caltagirone, J.-P. 2005. A fine-scale K - s model for atmospheric flow over heterogeneous Landscapes. Environmental Fluid Mechanics. 5: 247-265. doi: 10.1007/sl 0652-004-2124-x Fox, D.G. 1981. Judging air quality model performance: A summary of the AMS Workshop on Dispersion Model Performance. Bulletin of the American Meteorological Society. 62: 599-609. doi: 10.1175/1520-0477(1981)062<0599:JAQMP>2.0.CO;2 Frank, C. and Ruck, B. 2008. Numerical study of airflow over forest clearings. Forestry. 81(3): 259-277. doi: 10.1093/forestry/cpn031 Gardiner, B.A., Stacey, G.R., Belcher, R.E. and Wood, C.J. 1997. Field and wind tunnel assessments of respacing and thinning of tree stability. Forestry. 70(3): 233-252. doi: 10.1093/forestry/70.3.233 Gardiner, B., Marshall, B., Achim, A., Belcher, R. and Wood, C. 2005. The stability of different silvicultural systems: a wind-tunnel investigation. Forestry. 78(5): 471-484. doi: 10.1093/forestry/cpi053 Gash, J.H.C. 1986. Observations of turbulence downwind of a forest-heath interface. Boundary-Layer Meteorology. 36(3): 227-237. doi: 10.1007/BF00118661 Green, S.R. 1992. Modelling turbulent airflow in a stand of widely spaced trees. Phoenics journal of computation fluid dynamics and its applications. 5(3): 294-312. Hinze, J.O. 1975. Turbulence. Toronto: McGraw-Hill Publishing Co. 221 Irvine, M.R., Gardiner, B.A., and Hill, M.K. 1997. The evolution of turbulence across a forest edge. Boundary-Layer Meteorology. 84(3): 467-196. doi: 10.1023/A: 1000453031036 Jones, W.P. and Launder, B.E. 1972. The prediction of laminarization with a 2-equation model of turbulence. International Journal of Heath and Mass Transfer. 15(2): 301314. doi:10.1016/0017-9310(72)90076-2 Judd, M.J., Raupach, M.R., and Finnigan, J.J. 1996. A wind tunnel study of turbulent flow around single and multiple windbreaks, part I: Velocity fields. Boundary-Layer Meteorology. 80(1-2): 127-165. doi: 10.1007/BF00119015 Lakehal, D. 1999. Computation of turbulent shear flows over rough-walled circular cylinders. Journal of Wind Engineering and Industrial Aerodynamics. 80(1): 47-68. doi: 10.1016/SO167-6105(98)00122-6 Landsberg, J.J. and Jarvis, P.G. 1973. A numerical investigation of the momentum balance of a spruce forest. Journal of Applied Ecology. 10(2): 645-655. Launder, B.E., Morse, A.P., Rodi, W. and Spalding, D.B. 1972. Free Turbulent Shear Flows Vol. 1, NASA, Langley. p. 361-426 SP 321 . Launder, B.E. and Spalding, D.B. 1974. The numerical computation of turbulent flows. Computer methods in applied mechanics and engineering. 3(2): 269-289. doi: 10.1016/0045-7825(74)90029-2 Lee, X. 2000. Air motion within and above forest vegetation in non-ideal conditions. Forest Ecology and Management. 135:3-18. doi: 10.1016/S0378-1127(00)00294-2 Li, Z., Lin, J.D., and Miller, D.R. 1990. Air flow over and through a forest edge: A steadystate numerical simulation. Boundary-Layer Meteorology. 51(1-2): 179-197. doi: 10.1007/BFOO 120467 Liu, J., Chen, J.M., Black, T.A., and Novak, M.D. 1996. E-s modelling of turbulent air flow downwind of a model forest edge. Boundary-Layer Meteorology. 77(1): 21-44. doi: 10.1007/BF00121857 Mathieu, J. and Scott, J. 2000. An Introduction to Turbulent Flow. New York: Cambridge Univerity Press, p. 4. Meroney, R.N. 1968. Characteristics of wind and turbulence in and above model forests. Journal of Applied Meteorology. 7: 780-788. 222 Miller, D.R., Lin, J.D., and Lu, Z.N. 1991. Air flow across an alpine forest clearing: A model and field measurements. Agricultural and Forest Meteorology. 56: 209-225. Mitchell, S.J. 1995. A synopsis of windthrow in British Columbia: occurrence, implications, assessment and management. In Wind and Trees. Courts, M.P. and Grace, J. (editors). Cambridge University Press, New York, USA. p 448-459. Novak, M.D., Orchansky, A.L., Adams, R., Chen, W., and Ketler, R. 1997. Wind and temperature regimes in the B-5 clearing at the Sicamous Creek silvicultural systems research area: Preliminary results from 1995. In Sicamous Creek Silvicultural Systems Project: Workshop Proceedings. Hollstedt, C. and Vyse, A. (editors). April 24-25, 1996. Kamloops, British Columbia, Canada. Res. Br., B.C. Min. For., Vicotoria, B.C. Work Pap. 24/1997. p. 45-56. Novak, M.D., Warland, J.S., Orchansky, A.L., Ketler, R. And Green, S. 2000. Wind tunnel and field measurements of turbulent flow in forests. Part I: Uniformly thinned stands. Boundary-Layer Meteorology. 95(3): 457-495. doi: 10.1023/A: 1002693625637 Novak, M.D., Orchansky, A.L., Warland, J.S. and Ketler, R. 2001. Wind tunnel modelling of partial cuts and cutblock edges for windthrow. In Windthrow Assessment and Management in British Columbia: Proceedings of the Windthrow Researchers Workshop. Mitchell, S.J. and Rodney, J. (editors). January 31-February 1, 2001. Richmond, British Columbia, Canada, p. 176-192. Panferov, O. and Sogachev, A. 2008. Influence of gap size on wind damage variable in a forest. Agricultural and Forest Meteorology. 148:1869-1881. doi: 10.1016/j.agrformet.2008.06.012 Parish, R. 1997. Age and size structures of the forest at Sicamous Creek. In Sicamous Creek Silvicultural Systems Project: Workshop Proceedings. Hollstedt, C. and Vyse, A. (editors). April 24-25, 1996. Kamloops, British Columbia, Canada. Res. Br., B.C. Min. For., Vicotoria, B.C. Work Pap. 24/1997. p. 16-31. Puttonen, P. and Murphy, B. 1997. Developing silvicultural systems for sustainable forestry in British Columbia. In Sicamous Creek Silvicultural Systems Project: Workshop Proceedings. Hollstedt, C. and Vyse, A. (editors). April 24-25, 1996. Kamloops, British Columbia, Canada. Res. Br., B.C. Min. For., Vicotoria, B.C. Work Pap. 24/1997. p. 1-4. Quine, C.P. 1995. Assessing the risk of wind damage to forests: practice and pitfalls. In Wind and Trees. Courts, M.P. and Grace, J. (editors). Cambridge University Press, New York, USA. p. 379-403. 223 Quine, C.P. and Bell, P.D. 1998. Monitoring of windthrow occurrence and progression in spruce forests in Britain. Forestry. 71(2): 87-97 . doi: 10.1093/forestry/71.2.87-a Raupach, M.R. and Thorn, A.S. 1981. Turbulence in and above plant canopies. Annual Review of Fluid Mechanics. 13: 97-129. doi: 10.1146/annurev.fl.l3.010181.000525 Raynor, G.S. 1971. Wind and Temperature Structure in a Coniferous Forest and a Contiguous Field. Forest Science. 17(3): 351-363. Rollinson, T.J.D. 1987. Thinning control of conifer plantations in Great Britain. Annales des Sciences Forestieres. 44(1): 25-34. doi: 10.1051/forest: 19870103 Rowan, C.A., Mitchell, S.J., and Temesgen, H. 2003. Effectiveness of clearcut edge windfirming treatments in coastal British Columbia: short-term results. Forestry. 76(1): 55-64. doi: 10.1093/forestry/76.1.55 Ruck, B. and Frank, C. 2007. Numerical study of the airflow over clearings. Presented at International Union of Forest Researchers Organisation, International Conference on Wind and Trees (IUFRO 8.01.11), August 5-9, 2007, Vancouver, BC. Rudnicki, M., Mitchell, S.J., and Novak, M.D. 2004. Wind tunnel measurements of crown streamlining and drag relationships for three conifer species. Canadian Journal of Forest Research. 34(3): 666-676. doi: 10.1139/x03-233 Ruel, J.-C, Pin, D., and Cooper, K. 2001. Windthrow in riparian buffer strips: effect of wind exposure, thinning and strip width. Forest Ecology and Management. 143: 105113. doi: 10.1016/S0378-1127(00)00510-7 Shih, T.-H., Liou, W.W., Shabbir, A., Yang, Z. and Zhu, J. 1995. A new K-8 eddy viscosity model for high reynolds number turbulent flows. Computers & Fluids. 24(3): 227-238. doi: 10.1016/0045-7930(94)00032-T Stacey, G.R., Belcher, R.E., Wood, C.J., and Gardiner, B.A. 1994. Wind flows and forces in a model spruce forest. Boundary-Layer Meteorology. 69(3): 311-334. doi: 10.1007/BF00708860 Stathers, R.J., T.P. Rollerson, and S.J. Mitchell. 1994. Windthrow handbook for British Columbia forests. B.C. Ministry of Forests, Victoria, B.C. Working Paper 9401. 224 Taylor, K.E. 2000. Summarizing multiple aspects of model performance in a single diagram. Program for Climate Model Diagnosis and Iintercomparison. University of California, Lawrence Livermore National Laboratory, Livermore, California, USA. Report No. 55. Thorn, A.S. 1971. Momentum absorption by vegetation. The Quarterly Journal of the Royal Meteorological Society. 97:414-428. doi: 10.1002/qj .49709741404 Vyse, A. 1997. The Sicamous Creek silvicultural systems project: How the projects came to be and what it aims to accomplish. In Sicamous Creek Silvicultural Systems Project: Workshop Proceedings. Hollstedt, C. and Vyse, A. (editors). April 24-25, 1996. Kamloops, British Columbia, Canada. Res. Br., B.C. Min. For., Vicotoria, B.C. Work Pap. 24/1997. p. 4-14. Willmott, C.J. 1981. On the validation of models. Physical Geography. 2(2): 184-194. Willmott, C.J. 1982. Some comments on the evaluation of model performance. Bulletin American Meteorological Society. 63(11): 1309-1313. Willmott, C.J. and Wicks, D.E. 1980. An empirical method for the spatial interpolation of monthly precipitation within California. Physical Geography. 1: 59-73. Wilson, J.D. 1988. A Second-Order Closure Model for Flow through Vegetation. Boundary-Layer Meteorology. 42(3): 371-392. doi: 10.1007/BF00121591 Wilson, J.D. and Flesch, T.K. 1999. Wind and remnant tree sway in forest cutblocks. III. a windflow model to diagnose spatial variation. Agricultural and Forest Meteorology. 93(4): 259-282. doi: 10.1016/S0168-1923(98)00121-X Yakhot, V. and Orszag, S.A. 1986. Renormalization group analysis of turbulence. I. Basic theory. Journal of Scientific Computing. 1 (1): 3-51. doi: 10.1007/BF01061452 Yang, B., Shaw, R.H., and Pau U, K.T. 2006. Wind loading on trees across a forest edge: A large eddy simulation. Agricultural and Forest Meteorology. 141: 133-146. doi 10.1016/j.agrformet.2006.09.006 225