A Study on Single Tree Modeling for Wind- Induced Fluid–Structure Interaction Exploring the Level of Detail for Tree Representation Master’s thesis in Applied Mechanics, MSc (MPAME) John Paul Perumadan Mannummal Department of Industrial and Materials Science CHALMERS UNIVERSITY OF TECHNOLOGY Gothenburg, Sweden 2025 Master’s thesis in Applied Mechanics, MSc (MPAME) 2025 A Study on Single Tree Modeling for Wind-Induced Fluid–Structure Interaction Exploring the Level of Detail for Tree Representation John Paul Perumadan Mannummal Department of Industrial and Materials Science Chalmers University of Technology Gothenburg, Sweden 2025 A Study on Single Tree Modeling for Wind-Induced Fluid–Structure Interaction: Exploring the Level of Detail for Tree Representation Master’s thesis in Applied Mechanics, MSc (MPAME) John Paul Perumadan Mannummal © John Paul Perumadan Mannummal, 2025. Supervisor: Franziska Hunger, franziska.hunger@fcc.chalmers.se & Gustav Kettil, gustav.kettil@fcc.chalmers.se Examiner: Knut Andreas Meyer, knut.andreas.meyer@chalmers.se Master’s thesis in Applied Mechanics, MSc (MPAME) 2025 Department of Industrial and Materials Science Chalmers University of Technology SE-412 96 Gothenburg Telephone +46 31 772 1000 Cover: A 12 branch partitioned tree crown model. Typeset in LATEX Gothenburg, Sweden 2025 iv A Study on Single Tree Modeling for Wind-Induced Fluid–Structure Interaction: Exploring the Level of Detail for Tree Representation Master’s thesis in Applied Mechanics, MSc (MPAME) John Paul Perumadan Mannummal Department of Industrial and Materials Science Chalmers University of Technology Abstract Trees and green spaces are vital for sustainable urban living, yet tools to guide opti- mal tree placement and assess wind related risks remain limited. This study presents a numerical modeling framework for simulating wind-induced fluid-structure inter- action (FSI) on a single tree. Three tree models of increasing level of detail (LoD) are developed, and detailed FEM and CFD parameter studies are conducted to iden- tify key mechanical and aerodynamic properties. Rather than relying on a detailed CAD model, the tree geometry is defined through spatial functions, allowing for the assignment of material and permeability properties to stem and crown regions within a single mesh used by both solvers. The immersed boundary method (IBM) is employed to capture tree motion in the CFD simulations. Through FEM parameter studies, a set of material properties for the artificial crown is proposed, identifying that a stiffness on the order of 102 Pa is necessary to observe wind-induced defor- mations. The study also evaluated how modeling the stem as a permeable object rather than a solid body affects the aerodynamic forces and wake properties. It was found that a permeable stem yields aerodynamic forces and velocity fields compa- rable to those of a solid stem. However for both types of stem, the aerodynamic influence of the stem is confined to its immediate vicinity. Beyond the crown’s radial extent, the presence or absence of a stem showed negligible impact on the flow field. The results further highlight the dominant role of crown permeability in determin- ing aerodynamic loading and wake characteristics, underscoring the importance of measurement-based permeability inputs for crown. Overall, the study demonstrates the feasibility and current limitations of the proposed modeling approach, providing a foundation for future FSI based tree–wind interaction studies. Keywords: Tree Modeling, Tree-Crown Modeling, Fully Permeable Tree, Level of Detail, Wind-Tree Interaction, Fluid Structure Interaction, CFD, FEM v Acknowledgements I would like to sincerely thank my supervisors and examiner for their immense and unwavering guidance throughout this project. Their commitment to quality and dedicated involvement made this a true team effort. This thesis is part of the DUT project Urban ElemenTREE with financial sup- port from the Swedish Research Council for Sustainable Development Formas under the grant 2024-01221. This work is also part of the Digital Twin Cities Centre supported by Sweden’s Innovation Agency Vinnova under Grant No. 2024-03904. John Paul Perumadan Mannummal Gothenburg, August 2025 vii Acronyms ABL Atmospheric Boundary Layer. BT Branched Tree. CFD Computational Fluid Dynamics. DBH Diameter at Breast Height. FCC Fraunhofer-Chalmers Centre. FEM Finite Element Method. FSI Fluid–Structure Interaction. GS-FSI Gauss-Seidel FSI. IB Immersed Boundary. IBM Immersed Boundary Method. IBOFlow Immersed Boundary Octree Flow Solver. LAD Leaf Area Density. LaStFEM FEM Structural Solver. LoD Level of Detail. MRE Mean Relative Error. PAD Plant Area Density. PAD-crown Plant Area Density of Crown Cell. PAD-wood Plant Area Density of Wood Cell. PT Partitioned Tree. RANS Reynolds-Averaged Navier-Stokes. ST Simple Tree. TLS Terrestrial Laser Scanning. VIVs Vortex Induced Vibrations. ix List of Acronyms x Contents Nomenclature xiii List of Figures xv List of Tables xix 1 Introduction 1 1.1 Aim & Goals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 Limitations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2 Theory 5 2.1 Wind-Tree Interaction . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.1.1 Tree Nomenclature . . . . . . . . . . . . . . . . . . . . . . . . 5 2.1.2 Plant Area Density . . . . . . . . . . . . . . . . . . . . . . . . 6 2.1.3 Scales of Motion in Trees . . . . . . . . . . . . . . . . . . . . . 6 2.1.4 European Aspen Tree (Populus tremuloides) . . . . . . . . . . 7 2.2 Finite Element Method . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.2.1 Small Strain Elasticity . . . . . . . . . . . . . . . . . . . . . . 8 2.2.2 Governing Equations in Structural Mechanics . . . . . . . . . 9 2.2.3 Weak Formulation . . . . . . . . . . . . . . . . . . . . . . . . 9 2.2.4 Finite Element Discretization and Equilibrium Solver . . . . . 10 2.3 Computational Fluid Dynamics . . . . . . . . . . . . . . . . . . . . . 12 2.3.1 Immersed Boundary Method . . . . . . . . . . . . . . . . . . . 12 2.3.2 Governing Equations in Fluid Dynamics . . . . . . . . . . . . 13 2.3.2.1 k-ϵ Realizable Turbulence Model . . . . . . . . . . . 13 2.3.3 Source Terms from Permeability Model Used in Study . . . . 14 2.3.4 Atmospheric Boundary Layer Flow . . . . . . . . . . . . . . . 15 2.4 Fluid Structure Interaction . . . . . . . . . . . . . . . . . . . . . . . . 15 2.4.1 FSI Coupling Algorithms . . . . . . . . . . . . . . . . . . . . . 15 2.4.1.1 Monolithic & Staggered Approach . . . . . . . . . . 15 2.4.1.2 Strong (semi-monolithic) & Weak (explicit) Coupling 16 2.4.2 Gauss-Seidel FSI Solver . . . . . . . . . . . . . . . . . . . . . 16 3 Methodology 19 3.1 Tree Representation Geometries . . . . . . . . . . . . . . . . . . . . . 19 3.1.1 Simple Tree . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 3.1.2 Branched Tree . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 xi Contents 3.1.3 Partitioned Tree . . . . . . . . . . . . . . . . . . . . . . . . . . 22 3.2 Mesh Pre-Processing Steps . . . . . . . . . . . . . . . . . . . . . . . . 22 3.2.1 Single Object Mesh Generation . . . . . . . . . . . . . . . . . 23 3.2.1.1 Marking of Wood and Crown Elements . . . . . . . . 23 3.2.2 Local Mesh Refinement Along Stem & Branches . . . . . . . . 24 3.2.3 Partitioning of Tree Crown . . . . . . . . . . . . . . . . . . . . 24 3.2.4 Mesh Connectivity Check . . . . . . . . . . . . . . . . . . . . 26 3.2.5 Surface Mesh Extraction for IB-tree Local Refinement . . . . . 27 3.3 FEM Simulations in LaStFEM . . . . . . . . . . . . . . . . . . . . . . 28 3.3.1 FEM Simulation Setup . . . . . . . . . . . . . . . . . . . . . . 29 3.3.1.1 Material Properties for Artificial Crown Material . . 29 3.3.1.2 FEM Boundary Conditions . . . . . . . . . . . . . . 29 3.3.2 FEM Mesh Study . . . . . . . . . . . . . . . . . . . . . . . . . 29 3.3.2.1 Mesh Study Using ST . . . . . . . . . . . . . . . . . 30 3.3.2.2 Mesh Study Using BT . . . . . . . . . . . . . . . . . 32 3.3.3 FEM Parameter Study . . . . . . . . . . . . . . . . . . . . . . 32 3.3.4 Results & Discussions from FEM Parameter Study . . . . . . 32 3.3.4.1 Stiffness Ratio . . . . . . . . . . . . . . . . . . . . . 32 3.3.4.2 Stiffness Ratio Field . . . . . . . . . . . . . . . . . . 34 3.3.4.3 Poisson Ratio . . . . . . . . . . . . . . . . . . . . . . 36 3.3.4.4 Density Ratio . . . . . . . . . . . . . . . . . . . . . . 37 3.3.5 Base Case Deformations in ST, BT & PT . . . . . . . . . . . 38 3.4 CFD Simulations in IBOFlow . . . . . . . . . . . . . . . . . . . . . . 39 3.4.1 CFD Simulation Setup . . . . . . . . . . . . . . . . . . . . . . 40 3.4.2 CFD Mesh Study . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.4.3 CFD Parameter Study . . . . . . . . . . . . . . . . . . . . . . 45 3.4.4 Results & Discussions from CFD Parameter Study . . . . . . 45 3.4.4.1 PAD for Wood Cells . . . . . . . . . . . . . . . . . . 45 3.4.4.2 PAD for Crown Cells . . . . . . . . . . . . . . . . . . 48 3.4.4.3 Relative Influence of Stem & Crown . . . . . . . . . 50 4 Wind - Tree Interaction FSI Simulation Results 55 4.1 Discussion on Parameter Study Results . . . . . . . . . . . . . . . . . 55 4.2 One Way Coupled Simulations . . . . . . . . . . . . . . . . . . . . . . 57 4.3 Two Way Coupled Simulations . . . . . . . . . . . . . . . . . . . . . . 60 4.3.1 Two Way Coupled Simulation Setup . . . . . . . . . . . . . . 60 4.3.2 Results & Discussions . . . . . . . . . . . . . . . . . . . . . . 60 5 Conclusion 63 5.1 Summary of Methodology . . . . . . . . . . . . . . . . . . . . . . . . 63 5.2 Summary of Results . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 5.3 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 xii Nomenclature General APAD : Plant area density (PAD) [m2/m3] Stiffness_ratio : Ratio of Young’s modulus of crown material to that of the stem material Density_ratio : Ratio of density of crown material to that of the stem material FEM Theory E : Young’s modulus [Pa] ν : Poisson’s ratio σ : Stress vector in Voigt notation [Pa] ϵ : Strain vector in Voigt notation ∇̃ : Matrix differential operator [1/m] b : Body force vector [N/m3] u : Exact local displacement vector [m] δu : Virtual displacement vector [m] x : Position vector [m] Ui : Nodal displacement at node i [m] U : Global nodal displacement vector [m] δuh : Virtual displacement field (FEM approx.) [m] δU : Virtual nodal displacement vector [m] B : Strain-displacement matrix [1/m] D : Constitutive matrix [Pa] Ke : Element stiffness matrix [N/m] F e : Element force vector [N] K : Global stiffness matrix [N/m] F : Global force vector [N] C0 : Space of continuous functions Ω : Domain t : Traction vector [Pa] xiii Contents Γt, Γu : Neumann / Dirichlet boundaries nnode : Number of nodes per element Ni(x) : Shape function for node i N (x) : Vector of shape functions ai, bi, ci, di : Coefficients of linear shape function Ni Ωe : Domain of element e Γe t : Neumann boundary of element e CFD Theory u : Local velocity vector [m/s] p : Local pressure [Pa] ρ : Fluid density [kg/m3] ν : Kinematic viscosity [m2/s] νt : Turbulent (eddy) viscosity [m2/s] k : Turbulent kinetic energy [m2/s2] ϵ : Turbulent dissipation rate [m2/s3] σk : Turbulent Prandtl number for k transport c1ϵ, c2ϵ : Model constants Cµ : Eddy viscosity model constant A0, As : Model constants Ωij : Rotation rate tensor [1/s] Su : Momentum source term [N/m3] Sϵ : Turbulent dissipation source term [m2/s4] Cd : Drag coefficient βp, βd : Model coefficients Cϵ4, Cϵ5 : Model constants z : Height above ground [m] u∗ ABL : Friction velocity in atmospheric boundary layer [m/s] κ : von Karman constant Href : Reference height [m] uref : Velocity at reference height [m/s] z0 : Aerodynamic roughness length [m] xiv List of Figures 2.1 Tree nomenclature used in the study. . . . . . . . . . . . . . . . . . . 5 2.2 PAD distribution with leaves (left) and without leaves (right) esti- mated for single tree using TLS and varying grid resolution, taken from [8]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.3 Open grown Swedish columnar Aspen (left) and forest grown (right) American Aspen tree [38]. . . . . . . . . . . . . . . . . . . . . . . . . 7 2.4 Cartesian grid with immersed boundary showing different cell types. . 12 2.5 Working of a loosely coupled partitioned Gauss-Seidel FSI solver. . . 16 3.1 Tree representations. (a) Simple tree ; (b) Branched tree ; (c) Parti- tioned tree. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 3.2 Overall form and scale of the tree representation models (a) simple- tree overlaid on reference tree geometry; (b) reference forest grown Aspen tree. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.3 Variation of stem diameter with height for the chosen taper rate. . . . 22 3.4 Schematic diagram of exaggerated stem region and second branch region showing parameters involved in the geometric criteria defining tree-skeleton region. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.5 Wood elements (marked in red) resolved within the volume mesh by the element-locating-algorithm following local refinement of input mesh using a threshold maximum edge length of 3 cm. (right) Zoomed in view of resolved main-stem showing the resultant rough surface . . 25 3.6 Partitioned crown elements in PT. (a) Viewed in YZ plane (b) Viewed in XZ plane . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.7 Resolved stem before and after refinement and connectivity check. (a) before refinement and connectivity check (b) after refinement and connectivity check . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.8 Extracted tree-skeleton surface (marked in white) inserted into carte- sian mesh to perform local refinement. Viewed in X-Z plane. . . . . . 27 3.9 PAD distribution in IB-tree defining wood cells with PAD = 50 m2·m−3 and crown cells with PAD = 0.1 m2·m−3. (a) ST PAD dis- tribution (b) BT PAD distribution . . . . . . . . . . . . . . . . . . . 28 3.10 Comparison of meshes for different max element size for stem elements. 30 3.11 Comparison of meshes for different max element size for crown elements. 31 xv List of Figures 3.12 Deformation of crown viewed in X − Y plane passing through center of the crown. (a) stiffness ratio = 10−6 (b) stiffness ratio = 10−7 (c) stiffness ratio = 10−8 . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 3.13 Equivalent strain plots of crown viewed in X − Y plane passing through center of the crown. (a) stiffness ratio = 10−6 (b) stiffness ratio = 10−7 (c) stiffness ratio = 10−8 . . . . . . . . . . . . . . . . . 33 3.14 Equivalent strain plot of crown section viewed in X-Z plane passing through center of the crown for stiffness ratio = 10−8, density ratio = 10−8 and Poisson ratio = 0.44 . . . . . . . . . . . . . . . . . . . . . 34 3.15 Characteristic stiffness of elements across the discretized crown vol- ume consisting of 50 levels (a) ST stiffness field (b) BT stiffness field . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 3.16 Equivalent strain plots in deformed state of elements across the dis- cretized crown volume consisting of 50 levels (a) ST stiffness field (b) BT stiffness field . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.17 Equivalent strain plots of crown viewed in X − Y plane passing through center of the crown. (a) Poisson’s ratio = 0.44 (b) Poisson’s ratio = 0.0 (c) Poisson’s ratio = -0.2. . . . . . . . . . . . . . . . . . . 36 3.18 Deformation under self-weight of the crown viewed in the X–Z plane for the estimated density ratio of 0.00162, (stiffness ratio = 10−8, Poisson ratio = 0.0). . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.19 Pure FEM base case deformations in ST, BT, and PT, viewed in the X–Z plane. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 3.20 Pure FEM base case deformations in ST, BT, and PT, viewed in the X–Y plane. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 3.21 Overall dimensions of computational domain in CFD simulations with boundary conditions used. . . . . . . . . . . . . . . . . . . . . . . . . 40 3.22 Velocity, turbulent kinetic energy, dissipation rate inlet profiles gener- ated according to standard ABL wall functions.(z0 = 0.0018, u∗ ABL = 0.73 m/s, κ = 0.42). . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.23 Velocity contour plots for the benchmark cases using ST solid stem and global mesh parameters of Mesh M3 . . . . . . . . . . . . . . . . 42 3.24 Velocity contour plots for different mesh resolutions of IB-tree with PAD for wood cells = 50 m2·m−3 and PAD of crown cells equal to 0.1 m2·m−3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 3.25 Pressure contour plots different mesh resolutions of IB-tree with PAD for wood cells = 50 m2·m−3 and PAD of crown cells equal to 0.1 m2· m−3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 3.26 Through-flow velocity and streamwise pressure force component vari- ations from PAD-wood parameter study results . . . . . . . . . . . . 46 3.27 Wake velocity profiles for solid stem and PAD based permeable stem (PAD-wood = 50 m2·m−3 & 90 m2·m−3) plotted along stations down- stream the stem.(Plots are based on mesh M4) . . . . . . . . . . . . . 47 3.28 Wake turbulent kinetic energy profiles for solid stem and PAD based permeable stem (PAD-wood = 50 m2·m−3 & 90 m2·m−3) plotted along stations downstream the stem. (Plots are based on mesh M4) . . . . 48 xvi List of Figures 3.29 Average through-flow velocity and streamwise pressure force compo- nent variations from PAD-crown parameter study results . . . . . . . 49 3.30 Wake velocity profiles for different PAD-crown values, plotted along stations downstream the stem. . . . . . . . . . . . . . . . . . . . . . . 50 3.31 Velocity contour plots for different PAD-crown magnitudes, PAD- wood = 50 m2·m−3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 3.32 Pressure contour plots for different PAD-crown magnitudes, PAD- wood = 50 m2·m−3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 3.33 Velocity contour plot comparison between homogeneous ST, solid- stem ST, ST and BT. . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 3.34 Wake velocity profiles for different tree cases discussed in 3.9. . . . . . 54 4.1 Deformation plots in X-Z and X-Y plane for ST & BT from one-way coupled simulations. . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 4.2 Equivalent stress plots in X-Z and X-Y plane for BT and ST stem from one-way coupled simulations. . . . . . . . . . . . . . . . . . . . . 59 4.3 Deformation plots in X-Y plane passing through crown center from PT one-way coupled simulations. . . . . . . . . . . . . . . . . . . . . 60 4.4 Deformation plots in X-Z plane passing though center of the crown of BT from two-way coupled simulations. . . . . . . . . . . . . . . . . 61 4.5 Time series of pressure force components in streamwise and flow- normal directions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 xvii List of Figures xviii List of Tables 3.1 Material properties of Aspen wood used in the simulations. . . . . . . 29 3.2 Initial set of meshes comparing threshold maximum element size for wood elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 3.3 Second set of meshes comparing maximum element size of the crown elements for max edge length of stem equal to 0.05 m . . . . . . . . . 31 3.4 Second set of meshes comparing maximum element size of the crown elements using BT for max edge length stem equal to 0.05 m . . . . . 32 3.5 Mesh study using solid-object stem of ST . . . . . . . . . . . . . . . . 43 3.6 Mesh study using fluid-object IB-tree of ST with PAD-wood = 50 m2·m−3, PAD-crown = 0.1 m2·m−3 . . . . . . . . . . . . . . . . . . . 45 3.7 Results from parameter study of PAD value for wood cells using ST stem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 3.8 Results from parameter study of PAD value for crown cells for fixed PAD-wood of 50 m2·m−3 . . . . . . . . . . . . . . . . . . . . . . . . 49 3.9 Forces and effective PAD magnitudes of compared IB-trees . . . . . . 52 xix List of Tables xx 1 Introduction According to data published by the United Nations (UN), the urban population is steadily increasing and is expected to exceed 50 percent of the global population by 2030 [23]. In response, urban expansion is also happening rapidly. In this context, ensuring the sustainability of urban living spaces is both critical and challenging. The World Health Organization (WHO) recommended per-capita area of urban green space for the physical and mental well-being of a person is 9 m2 [31]. With cities responsible for 70 percent of global CO2 emissions [18], it is also critical to maintain sufficient vegetation to serve as a CO2 sink. Inadequate vegetation and poor planning of urban green spaces can negatively impact the urban microclimate [35, 28]. Realizing the vital role of vegetation in sustainable urban development, efforts are required to restore lost green spaces, manage existing ones effectively, and find a balance between green and built spaces during urban expansion. The importance of urban trees in this endeavor cannot be overstated. However, urban trees face many limitations related to the scarcity of space, high costs, and the often-challenging conditions reducing the longevity of trees. In Sweden, wind is a major factor influencing tree establishment, growth, and po- tential mechanical failure. However, our practical understanding of the interplay between wind and trees remains extremely limited. Consequently, there is limited information on optimal positioning when planting new trees. Further, our assess- ment of the risk for wind-induced tree failures, due to wind-throw (uprooting) or wind breakage (branch or trunk damage), often falls short. Therefore, urban decision makers, planners, and other stakeholders are in need of tools supporting evidence based urban greening strategies. To investigate the vulnerability of existing trees to the expansion of built spaces and to support the planning of new tree placement, Computational Fluid Dynam- ics (CFD) simulations in conjunction with Finite Element Method (FEM) based tree models can be employed. Different approaches to modelling trees, with varying Level of Detail (LoD), exist in the literature. In the context of urban microclimate simulations the required LoD of individual trees has been a topic of discussion over the recent years [44, 10, 11, 47]. With expanding access to methods like Terrestrial Laser Scanning (TLS) which em- ploys technologies like Light Detection and Ranging (LiDAR), the definition of LoD in individual trees are changing. An extensive review of existing definitions of LoD for trees is done by Xu et al. in their work [44]. According to Xu et al., a set of 1 1. Introduction tree parameters are used in defining the LoD for trees in the context of urban mi- croclimate simulations. These include tree dimension, canopy geometry, Leaf Area Density (LAD) and leaf type. In general, existing models for individual trees have primarily been employed to in- vestigate airflow patterns around trees and their effects on the urban microclimate [13, 1, 8, 29, 9, 5], as well as to study the mechanical response of trees to wind loading [14, 45, 33, 2]. Among the works focused on the former category, the LoD varied from 2D to volumetric representations, with partial or no consideration of the tree skeleton. Most of these models used permeability-based approaches involv- ing source term modifications to account for through-flow variations in tree crown. They typically relied on common assumptions such as constant aerodynamic and permeability properties like drag coefficient and LAD. In the latter category, the LoD ranged from simplified mechanical representations for FEM, such as cantilever beam models for the stem, to FEM models based on highly detailed tree skeletons reconstructed from TLS point cloud data. However, as pointed out by Xu et al [44]., both excessive geometric detail and over- simplification that neglects key physical characteristics can render existing models impractical and inefficient for urban microclimate simulations. The literature review highlights the absence of a robust and computationally efficient tree representation model suitable for such applications. From the literature review, it was also evident that, hardly any literature exists on modelling tree motion, considering a permeable crown. A recent work explore the bending of a simplified tree under simplified fluid conditions [2]. Yet, as trees adapt to the prevailing wind conditions, streamlining effects, crown reconfiguration and movement are important measures to reduce the force on the tree. Experiments of single tree measurements under different wind conditions have shown these effects [46, 22]. A numerical model that captures these effects enables further studies on forces exerted on trees in different urban environments, risk estimation of trees and might improve the prediction of the flow field downstream of the trees, being an important characteristic for urban planning. 1.1 Aim & Goals Develop a parametric numerical model of a single tree that balances level of detail and physical accuracy, suitable for analyzing its structural and aerodynamic responses to wind using Fluid–Structure Interaction (FSI) simulations. To achieve this, a FSI based modeling approach will be adopted for the proposed tree models. Emphasis is placed on streamlining the model generation pipeline to facilitate scalability of the study across different tree species and geometries. The goal will be accomplished through several sub-objectives that guide the devel- opment and integration of both structural and fluid dynamic components. 2 1. Introduction • Modelling of the tree as a single unified mesh, thereby eliminating the need for spatial coupling between the stem and crown in the FEM and CFD solvers. • Developing an FEM model capable of capturing the key dynamics of both the tree skeleton and the tree crown. • Conducting a FEM parameter study to identify and quantify the artificial material properties most relevant to accurately modeling the tree crown. • Developing a CFD model capable of capturing variation of flow properties over and through the tree-crown. • Conducting a CFD parameter study to identify appropriate permeability prop- erties for the permeable stem and crown within the single-object tree repre- sentation. • Analysis of the developed tree models using physical reasoning of structural and aerodynamic responses in FSI simulations. The LoD for the proposed tree representations includes the following aspects: • Dimension: True-to-scale tree geometry. • Tree skeleton: A tapered cylindrical cantilever stem with tapered cylindrical first-order branches. • Tree crown: Represented as a spherical volume with internal partitions; through-flow is modeled using source-term modification permeability models. • Plant Area Density (PAD): Defined as spatial fields within the crown volume. 1.2 Limitations The following aspects which are relevant in an actual tree are neglected or approxi- mated in this study: • Soil mechanics and related root-plate effects of the tree are neglected. • Instead of fully resolving the flow through the tree crown, permeability models are used to model variation of through flow properties. • The specific modelling approach approximates the solid behavior of stem using a fluid object with very low permeability. Further in the report, Chapter 2 discusses the relevant supporting theory for the structural and fluid dynamics concepts used in the development of the tree models. Chapter 3 details the development process, covering the FEM and CFD components of the framework. In Chapter 4, the results and discussion of the FSI simulations conducted using the developed tree models are presented. The report concludes with a discussion of the overall framework, key insights obtained, and suggestions for future work in Chapter 5. 3 1. Introduction 4 2 Theory 2.1 Wind-Tree Interaction In this section, a basic introduction to the bio-mechanical description of a tree, definition of PAD, different scales present in tree motion and relevant peculiarities of forest grown Aspen tree that led to the choice of it for the study are given. 2.1.1 Tree Nomenclature Figure 2.1: Tree nomenclature used in the study. Structurally, the root system anchors the tree and makes it stable under wind and other external loading. Structure of the tree extends from the root system into the stem and then into first order and higher order branches. Depending on the nature of external stimuli like wind direction, exposure to sunlight and presence of obstacles, the morphology of the tree structure will vary in order to adapt to such stimuli. Open grown trees are subjected to relatively higher wind loads, and this causes such trees to have thicker stem than forest grown trees which tend to have 5 2. Theory long and narrow stem. Also features like buttressing at the roots are predominant in open grown trees. Similar morphological evolution of trees is collectively termed thigmomorphogenesis [6]. Bole height is measured along the main-stem from the ground to the starting point (node) of the first branch. The crown comprises of the leaves, twigs and higher order branches above the bole height. The collective mass of leaves and twigs are often called foliage. The largest branches are called first order branches. And the branches that grows laterally from the first order branches are termed second order and so on. A discussion of the growth characteristics of a tree is presented by Franklin. J and Mercker. D in their work [16]. A pictorial representation of tree nomenclature is given in Figure 2.1. 2.1.2 Plant Area Density Figure 2.2: PAD distribution with leaves (left) and without leaves (right) estimated for single tree using TLS and varying grid resolution, taken from [8]. Plant Area Density (PAD) is defined as the sum of one-sided plant surface area (windward facing area) in unit volume. PAD has the units m2·m−3. In contrast to LAD in which only one-sided leaf area is considered, the PAD includes one-sided area of all material parts including leaves, branches and the stem. Measurements of PAD can be made using TLS method. A detailed study measuring PAD of single tree using TLS is done by Dellwik et al. in their work [8]. The measured PAD distributions for a single tree, taken from their work is given in Figure 2.2 for reference. 2.1.3 Scales of Motion in Trees The wind induced dynamic response of a tree occurs at different scales. This multi- scale response helps the trees to dissipate the kinetic energy transferred by the wind via aerodynamic loads. These scales of motion include: 6 2. Theory • Large scale bending and sway motion of the entire tree. • Intermediate scale buffeting motion of branches. • Small scale fluttering motion of leaves. In general, these motions get initiated when the aerodynamic force disturbs the static equilibrium position of the tree corresponding to no-wind condition. For steady wind conditions, the tree will try to re-adjust itself to a modified equilibrium position. With fluctuating wind flow, these motions can persist longer with different amplitudes. With complex architecture of a tree, higher number of modes and degrees of freedom can exist in the tree [15], making it often difficult to find an equilibrium state even in low wind conditions. Also, high wind conditions can result in high amplitude oscillations of branches or the entire tree by self excitation or resonance by Vortex Induced Vibrations (VIVs) [36]. Such motions can lead to structural failures like stem breakage or uprooting as well as create airflow patterns downstream. While the sway motions add to the canopy scale eddies, the buffeting of branches and motion of the leaves (collectively called foliage motion) create small scale eddies in the flow field. Within the foliage motion, the individual leaf motion predominates over effects due to motion of branches. Foliage motion makes the aerodynamic characteristics of the crown like drag coefficient and PAD, a function of time. This change of properties induced by foliage motion is termed reconfiguration of the crown. A clear discussion of the foliage dynamics and its effects on the wind is discussed in detail by Tadrist et al. in their work [36]. 2.1.4 European Aspen Tree (Populus tremuloides) Figure 2.3: Open grown Swedish columnar Aspen (left) and forest grown (right) American Aspen tree [38]. Aspen (Populus tremuloides) is a large, fast-growing broad-leaf tree. They are char- acterized by a very narrow growing stem. Based on the growing conditions, PAD distributions of trees will vary. Forest grown Aspen trees seems to have a slender form with small spherical crown at the top of the narrow stem, while open grown 7 2. Theory trees like those in urban setting usually have short stem and cylindrical crown shapes. Reference images showing the overall geometry of the columnar Swedish Aspen and a forest grown American Aspen tree are shown in Figure 2.3. Its unique ability of regrowth into multiple new trees from its own root system makes it attractive to the wood industry. A detailed discussion about the role of Aspen trees in serving bio-diversity, carbon sequestration, wood industry and other commercial needs is presented by Rogers et al. in their review article [30]. The choice of forest grown Aspen trees for the study is primarily due to its simplicity in structure and overall form. 2.2 Finite Element Method This section outlines the fundamental concepts relevant to the FEM simulations conducted in the study. A brief overview of the small-strain linear elasticity model which was used to describe the material behavior of the tree elements is provided. Further, the weak formulation leading to finite element discretization of the gov- erning equations for a linear elastic solid and a mathematical description of static equilibrium analysis in FEM are also presented. Supporting theory for FEM discussed in this report are based on concepts explained in the book Introduction to the Finite Element Method [24]. 2.2.1 Small Strain Elasticity When strains and rotations are small, and for materials that exhibit no strain- softening, a small-strain theory based on kinematics, equations of motion, and con- stitutive laws is commonly used. The basic assumptions of purely mechanical be- havior, reversibility, and path-independence imply the absence of energy dissipation during deformation. The energy expended in deforming the material is stored as strain energy, which is fully recoverable upon unloading. The regime of linear elastic behavior is typically confined to small strain magnitudes, which is why the small-strain elasticity theory is used to model such materials. For a linear elastic material, stress is related to strain via the constitutive relation known as Hooke’s law. For a general three-dimensional isotropic elastic material, this relation can be written in matrix (Voigt) notation as: σ = D ϵ (2.1) D = E (1 + ν)(1 − 2ν)  1 − ν ν ν 0 0 0 ν 1 − ν ν 0 0 0 ν ν 1 − ν 0 0 0 0 0 0 1−2ν 2 0 0 0 0 0 0 1−2ν 2 0 0 0 0 0 0 1−2ν 2  8 2. Theory Here, σ is the 6-component stress vector, ϵ is the 6-component strain vector, and D is the constitutive matrix containing the elastic moduli. E is the Young’s modulus and ν is the Poisson’s ratio. 2.2.2 Governing Equations in Structural Mechanics In the context of small-strain linear elasticity, the governing equations in their strong form consist of the equilibrium equation (or conservation of linear momentum), the constitutive relation linking internal stresses to strains, and the kinematic relation expressing strain in terms of displacement. These equations form the foundation of the finite element formulation and are given as follows: 1. Equilibrium Equation (Cauchy Momentum Balance) ∇̃T σ + b = 0 (2.2) where b is the body force vector per unit volume. ∇̃ is the matrix differential operator. 2. Constitutive Equation (Hooke’s Law): σ = Dϵ (2.3) 3. Kinematic Equation (Strain-Displacement Relation): ε = ∇̃u (2.4) where u is the displacement vector. Governing equations 2.2, 2.3 and 2.4 together with the respective Dirichlet and Neumann boundary conditions, constitute the strong formulation (i.e., the original set of differential equations) of the linear elasticity problem. The strong formulation must be satisfied point-wise throughout the entire domain and therefore defines the exact state of the system. 2.2.3 Weak Formulation While the strong form requires the solution to satisfy the differential equations point-wise and demands higher order continuity and differentiability, this is often too restrictive for complex geometries and discontinuous material properties. To overcome these limitations, we introduce the weak form, which is based on a variational principle. In the context of linear elasticity, this variational principle corresponds to the principle of virtual work. According to this principle, the total virtual work done by internal stresses is equal to the total virtual work done by external forces for all kinematically admissible virtual displacements. Instead of enforcing the equilibrium equation point-wise, the weak form ensures that the equation is satisfied in an integral sense over the domain. This approach allows us to relax the differentiability requirements of the trial solution. 9 2. Theory By integrating the equilibrium equation against a set of suitably chosen test func- tions (typically piecewise continuous functions from a C0 function space, such as linear Lagrange polynomials), the requirement for the displacement field to be twice differentiable is reduced to only once differentiable. This makes the method well- suited for piecewise polynomial approximations commonly used in the finite element method. Furthermore, by applying the Gauss divergence theorem, the spatial derivatives act- ing on the stress field are transferred to the test functions. This avoids the need to directly differentiate the stress components, which is particularly beneficial when stresses are derived from displacements via constitutive laws. The weak form, derived using the principle of virtual work, is obtained from the strong form equation (2.2) as follows: We begin by multiplying the equilibrium equation by a virtual displacement δu ∈ V (test function space), such that δu = 0 on the Dirichlet boundary Γu (where actual displacements are prescribed). We then integrate the resulting expression over the domain Ω: ∫ Ω δuT (∇ · σ) dΩ + ∫ Ω δuT b dΩ = 0 (2.5) Applying the divergence theorem (i.e., integration by parts) to the first term allows the derivative to be transferred from the stress tensor σ to the virtual displacement δu, reducing the differentiability requirements on the solution. This yields:: − ∫ Ω ε(δu)T σ dΩ + ∫ Γ δuT (σ · n) dΓ + ∫ Ω δuT b dΩ = 0 (2.6) Using the traction boundary condition σ · n = t̄ on the Neumann boundary Γt, and noting that virtual displacement δu = 0 on the Dirichlet boundary Γu, the boundary integral reduces to: − ∫ Ω ε(δu)T σ dΩ + ∫ Γt δuT t̄ dΓ + ∫ Ω δuT b dΩ = 0 (2.7) Substituting the constitutive relation σ = Dε(u) into the expression yields the weak form of the linear elasticity problem:∫ Ω ε(δu)T Dε(u) dΩ = ∫ Ω δuT b dΩ + ∫ Γt δuT t̄ dΓ (2.8) Rearranging the weak form reveals that it corresponds to the first variation of the total potential energy functional, which must vanish for equilibrium. Thus, the weak form can also be interpreted as the variational statement of the principle of minimum potential energy, which holds for linear elastic materials. 2.2.4 Finite Element Discretization and Equilibrium Solver To solve the weak form numerically using the finite element method, the unknown continuous displacement field u is approximated using a finite-dimensional basis function defined over the discretized domain. In this study, the domain is discretized 10 2. Theory using four-node linear tetrahedral elements. Within each element, the displacement field is interpolated using piecewise linear shape functions associated with the ele- ment’s nodes. Accordingly, the displacement field u and virtual displacement δu are approximated as: uh(x) = nnode∑ i=1 Ni(x)Ui = N (x)U (2.9) δuh(x) = nnode∑ i=1 Ni(x)δUi = N (x)δU (2.10) where U are the nodal displacement vectors, and N (x) is the vector of shape func- tions. For a linear tetrahedral element, the gradient of Ni is constant over the element and each Ni has the form: Ni(x, y, z) = ai + bix + ciy + diz The strain-displacement matrix B relates nodal displacements to strains via: ε(uh) = BU , ε(δuh) = BδU Substituting these finite-dimensional approximations into the weak form results in a matrix expression at the element level: δUT (∫ Ωe BT DB dΩ ) U = δUT (∫ Ωe NT b dΩ + ∫ Γe t NT t̄ dΓ ) (2.11) From the above equation, the element stiffness matrix and element force vector are defined as: Ke = ∫ Ωe BT DB dΩ (element stiffness matrix) (2.12) F e = ∫ Ωe NT b dΩ + ∫ Γe t NT t̄ dΓ (element force vector) (2.13) To form the global system, all element matrices and vectors are assembled based on the global numbering of nodes and degrees of freedom: KU = F (2.14) where: • K = ∑ e Ke is the global stiffness matrix, • U is the global displacement vector (unknowns), • F = ∑ e F e is the global force vector. For static and linear elasticity problems, the resulting sparse system is solved im- plicitly without time stepping, using a direct or iterative linear solver. Solving this system provides the displacement field that satisfies mechanical equilibrium under the given loading and boundary conditions. 11 2. Theory 2.3 Computational Fluid Dynamics This section provides an overview of the methods and models used in CFD simula- tions from the study. The Immersed Boundary Method (IBM), as implemented in the Immersed Boundary Octree Flow Solver (IBOFlow) solver developed at Fraunhofer-Chalmers Centre (FCC) [19], is first introduced to describe how ge- ometries are handled without the need for conformal meshing. Then, the governing equations solved by the steady Reynolds-Averaged Navier-Stokes (RANS) solver are presented along with the associated turbulence model. The permeability model used to represent the variation of flow properties within the tree crown is also discussed. Finally, the characteristics of the atmospheric boundary layer flow are outlined, as they are used to define the inlet boundary conditions for the CFD simulations in the study. 2.3.1 Immersed Boundary Method Immersed boundary method uses a fixed cartesian grid. The boundary (surface mesh) of the object to be analyzed is inserted into this cartesian grid and thus the need for a rigorous body-fitted fluid mesh is avoided. Mesh refinement can be easily applied to the surrounding cartesian grid to accurately capture the geometry and flow features near the immersed object. In applications involving motion of objects in a fluid domain and in situations where large displacements in fluid domain is plausible, IBM can be computationally efficient and affordable compared to a adaptive meshing approach. Figure 2.4: Cartesian grid with immersed boundary showing different cell types. Once the boundaries are immersed into the cartesian grid, appropriate boundary conditions using forcing functions or source terms are incorporated into the cells lying in the vicinity of the boundary. Different Immersed Boundary (IB) methods exists in literature [27, 39]. In the IBOFlow solver used in the study, a mirroring 12 2. Theory IB method [19] is used in which, the velocity is implicitly mirrored across the IB. To impose the no-slip boundary condition at the IB, the velocity at the mirroring cell is set such that a linear interpolation between the mirroring cell and the known velocity at the mirrored point yields exactly the local velocity of the IB. The fictitious velocities inside the solid are excluded from the continuity equation to ensure that no mass flux occurs across the IB. The different cell types used to define cells in IB method is shown in Figure 2.4. 2.3.2 Governing Equations in Fluid Dynamics For the study, a steady-state RANS CFD solver is used. Incompressible unsteady formulations of governing equations according to the model are discussed below [40]. RANS equations are solved in a segregated manner, with pressure velocity coupling using SIMPLEC method. 1. The continuity equation continuity equation or mass balance equation can be written in vector notation as: ∇u⃗ = 0 (2.15) Where u⃗ is the local velocity vector. 2. The momentum equation The momentum equation also called as Navier Stokes equation or momentum balance equation can be written in vector form as: ∂u⃗ ∂t + u⃗ · ∇u⃗ = −∇p ρ + ∇ · ( (ν + νt)∇u⃗ ) (2.16) where p is the local pressure, ρ is the fluid density, ν and νt are kinematic viscosity and turbulent viscosity respectively. 2.3.2.1 k-ϵ Realizable Turbulence Model To solve for mean flow and turbulence, RANS equations are solved in conjunction with the k-ϵ Realizable eddy viscosity turbulence model. The closure problem arising in RANS equations is addressed using the Boussinesq approximation and solving for the transport of turbulent kinetic energy and turbulent dissipation rate. The local turbulent kinetic energy is found by solving the transport equation for turbulent kinetic energy or the k-equation. And the realizable k-equation is given as: 3. The k-equation ∂k ∂t + u⃗ · ∇k = ∇ · ( (ν + νt σk )∇k ) + νtS 2 − ϵ (2.17) where S = √ 2SijSij, Sij is the strain rate tensor quantifying deformation of fluid volume, ϵ is the turbulent dissipation rate and σk is the turbulent Prandtl number. 13 2. Theory The local turbulent dissipation rate is found by solving the transport equation for ϵ and is given as: 4. The ϵ equation ∂ϵ ∂t + u⃗ · ∇ϵ = ∇ · ( (ν + νt σk )∇ϵ ) + c1ϵS − C2ϵ ϵ2 k + √ νϵ (2.18) νt = Cµ k2 ϵ (2.19) Cµ = 1 A0 + As kU∗ ϵ (2.20) where k is the local turbulent kinetic energy, η = S k ϵ ,c1ϵ = max ( 0.43, η η+5 ) , c1ϵ = 1.92, U∗ = √ SijSij + ΩijΩij, A0 = 4.04, As = √ 6 cos ϕ,W = SijSjkSki√ (SijSij)3 , ϕ = 1 3 arccos( √ 6W ). 2.3.3 Source Terms from Permeability Model Used in Study To approximate the variation in flow properties inside the tree crown, permeability model with source terms found in literature [10] are used. Su⃗ = −ρCdAP ADu⃗|u⃗| (2.21) Sϵ = ρCdAP AD ϵ k ( Cϵ4βp|u⃗|3 − Cϵ5βd|u⃗|k ) (2.22) where: ρ : density of fluid Cd : drag coefficient AP AD : plant area density βp : fraction of mean kinetic energy converted into turbulent kinetic energy u⃗ : local velocity vector βd : dimensionless coefficient for short circuiting of turbulent cascade k : local turbulent kinetic energy ϵ : local dissipation rate Cϵ4, Cϵ5 : modal constants The coefficients βp, βd, Cϵ4 and Cϵ5 used in the study are 1,1,0 and -2.556 respectively based on study by Dellwik et al.[8]. The above source terms are solved as additional terms in momentum, turbulent dissipation rate transport equations, at the internal cells of IB object. Since the momentum and turbulent source terms are all functions of PAD, by setting a zero PAD value in a region, the permeability model can be switched off in those cells. 14 2. Theory 2.3.4 Atmospheric Boundary Layer Flow The lower part of Atmospheric Boundary Layer (ABL) extending from the ground to an elevation of 200 m is characterized to be significant for studying different phenomena in urban microclimate, pollutant dispersion and deposition by wind and many more. In order to approximate the ABL flow, wall functions for velocity and turbulent flow properties are used in CFD. These wall functions are derived based on a number of assumptions like: • Horizontal homogeneity and absence of streamwise gradients. • Equilibrium between mean flow and turbulent properties with turbulence in- duced by ground roughness. • Local equilibrium exists between turbulent production and dissipation. Based on the above assumptions, following wall functions for fully developed ABL profiles of velocity, turbulent kinetic energy and turbulent dissipation rate are used [4]: U(z) = u∗ ABL κ ln(z + z0 z0 ) (2.23) k(z) = u∗2 ABL√ Cµ (2.24) ϵ(z) = u∗3 ABL κ(z + z0) (2.25) where: u∗ ABL : ABL friction velocity κ : von Karman constant (0.40–0.42) Cµ : turbulent viscosity coefficient from standard k-ϵ turbulence model (0.09) z : height measured from ground z0 : aerodynamic roughness length 2.4 Fluid Structure Interaction This section introduces the staggered FSI approach adopted in the study. The distinction between weak and strong coupling strategies is briefly discussed to clarify the level of interaction between the fluid and structural solvers. Finally, the Gauss- Seidel FSI (GS-FSI) solver used for the two-way coupled simulations is outlined. 2.4.1 FSI Coupling Algorithms 2.4.1.1 Monolithic & Staggered Approach In the monolithic approach, the governing equations from fluid and solid domains are solved as a combined system of equations. While in a staggered (partitioned) approach, separate solvers for fluid and solid variables are used. Here, the interaction between physics is ensured by the use of a coupling algorithm between solvers. The 15 2. Theory fundamental advantage of the staggered approach lies in the ability to use optimized solvers for the respective physics. 2.4.1.2 Strong (semi-monolithic) & Weak (explicit) Coupling The coupling between physics is said to be strong if the case demands synchroniza- tion of all the system variables during every time-step [41, 43]. Cases involving large deformations and small time scales, in which both domains influence each other are examples of strong coupling. In applications where time scale of deformation is large or the feedback to flow is insignificant, weak or loosely coupled approach can be used. In weak coupling, both solvers progress one time-step in each iteration and no synchronization of variables between solvers are made. Based on the nature of the coupling, the choice of coupling scheme is made. Figure 2.5: Working of a loosely coupled partitioned Gauss-Seidel FSI solver. 2.4.2 Gauss-Seidel FSI Solver In this study, a GS-FSI solver is used to couple the structural and flow solvers. The basic working of the GS-FSI solver in a general two way coupled simulation is described in the Figure 2.5. The solver works in a sequential manner transfer- ring aerodynamic forces and deformed tree geometry in each FSI iteration. Internal iterations are also included within each FSI iteration to ensure synchronization of variables between CFD and FEM and strong coupling [41]. During each FSI iter- ation, equilibrium of structural and flow properties is monitored using convergence of respective residuals. Once these residuals are within allowed tolerances for each solver, iterations are stopped. In the case of one-way coupling, the feedback from structural solver to the flow solver is absent. During the first iteration, flow solver calculates the steady state flow field and forces on the tree using the undeformed tree geometry. Using the 16 2. Theory solved pressure field, the structural solver solves for the static equilibrium geometry of the tree. 17 2. Theory 18 3 Methodology As discussed in the beginning of the report, the objective of the study is to develop a single tree model suitable for investigating wind-tree interaction using FSI sim- ulations. The proposed tree models and the FSI algorithm should avoid detailed modeling of tree morphology but instead at the same time be able to resolve the wind-tree interaction physics with sufficient realism. To achieve this, three tree ge- ometries of increasing level of detail are used to study different modelling aspects of the tree-model. Using the basic tree geometry, the relevant material and aero- dynamic parameters of the model are first established. These parameters are then progressively fine tuned for geometries with increased complexity. In general, three main aspects are considered in the investigated tree models. These are: 1. Overall form of the tree - influencing its external aerodynamics. 2. Branching architecture (tree-skeleton) - which defines the mass, stiffness and damping characteristics of tree motion dynamics. 3. Variation of through flow properties in the tree crown - which affect the wake field dynamics. The motivation and steps taken to choose different parameters of the tree-models, according to the above considerations and the planned objective are elaborated in the following sections. 3.1 Tree Representation Geometries The three tree geometries are named as Simple Tree (ST), Branched Tree (BT) and Partitioned Tree (PT) respectively. In contrast to tree models having a solid stem- object and fluid crown-object, the proposed tree-models in common were defined as a single-object in the structural-solver and flow-solver. This approach ensures ease of implementation of the moving tree in the coupled FSI setup by eliminating the need for spatial coupling between stem and crown. The discussed tree geometries are described in the following subsections starting with the ST. A graphical representation of the geometries are shown in Figure 3.1. Throughout the discussion, a right-handed cartesian coordinate system is followed with the tree at the origin, tree height along the positive Z axis and wind being uni-directional along the positive X axis. 19 3. Methodology (a) (b) (c) Figure 3.1: Tree representations. (a) Simple tree ; (b) Branched tree ; (c) Parti- tioned tree. 3.1.1 Simple Tree The overall form of the ST corresponded with the real-scale reference forest grown Aspen tree. It is defined using a tapered cylindrical cantilever stem with a contin- uous spherical crown. The ST is set to resolve the external aerodynamics of the tree as well as the primary bending mode under aerodynamic load. The material properties of the cantilever stem are taken as the properties of European Aspen wood from the reference article [12]. Within the ST geometry, the stem extends from root defined with a root diameter, to the tip of the crown with a defined stem taper rate. The minimum cross-section diameter of the stem is constrained at the tip using a limiter. Further the crown ra- dius and crown-center position are set to match the overall tree form. The diameter variation of the stem for the chosen taper rate is plotted in Figure 3.3. The magni- tude of taper rate used in the study is also given in the Figure 3.3. The Diameter at Breast Height (DBH) for the stem is measured as 37.825 cm. 20 3. Methodology (a) (b) Figure 3.2: Overall form and scale of the tree representation models (a) simple-tree overlaid on reference tree geometry; (b) reference forest grown Aspen tree. 3.1.2 Branched Tree In order to resolve the dynamics of the crown and the entire tree, the stiffness and mass-induced damping characteristics introduced by the branches must be consid- ered. Consequently, the first-order branches of the tree are included along different levels of the main stem in BT. The branches are also considered to be tapered cylinders with taper rate similar to that of the main stem. Node diameter of the branches are considered to be half the diameter of the main stem measured at the node height of the branch. The orientation of the branches are parameterized using angles alpha (defined between axis of a branch and the positive z axis) and theta (defined between axis of a branch and the positive x axis). For the BT in the study, alpha is considered to be 45◦ and theta to be 0◦, 90◦, 180◦ and 270◦ respectively for the four branches. To position the branches along different levels of the main stem, node positions of the branches are also incorporated as a parameter. The crown is defined as in ST. 21 3. Methodology Figure 3.3: Variation of stem diameter with height for the chosen taper rate. 3.1.3 Partitioned Tree The approximation of a continuous spherical crown restricts the degrees of freedom of the individual branches. In reality, motion of the branches are nearly independent due to morphological evolution and adaptive responses of the tree [6]. This is also vital to the tree response since, a significant amount of wind energy intercepted by the tree crown is dissipated through buffeting motion of branches [32]. Therefore, replacing the continuous spherical crown which is restrictive to the buffeting motion of the branches, partitions are introduced into the crown in PT. The objective of having partitions in the crown is to allow a branch and its associated crown region to move independent of other branches and the main stem. A detailed discussion of partitioning procedure is given in the section 3.2.3. The overall shape of the crown is retained as spherical with the same position and diameter as in the cases of ST & BT . 3.2 Mesh Pre-Processing Steps Having proposed the LoD, the next step was to numerically represent the different tree geometries. A single-object discretization of the tree is implemented in both FEM and CFD solvers. In this section, the steps which led to this single-object tree mesh representation are discussed. 22 3. Methodology 3.2.1 Single Object Mesh Generation The tree-mesh generation process starts with an input three dimensional, uniform and unstructured base volume-mesh with tetrahedron elements. The base mesh is defined using the overall form (volume) of the ST and was created using ANSA pre-processing software. 3.2.1.1 Marking of Wood and Crown Elements With the single-object mesh approach, the elements with properties of wood and those of crown have to be marked. This is done with the help of an algorithm (element-locating-algorithm) to locate and mark wood elements or cells based on a geometric criteria. In summary, the geometric criteria defines the volume of the reference tree-skeleton. Figure 3.4: Schematic diagram of exaggerated stem region and second branch region showing parameters involved in the geometric criteria defining tree-skeleton region. The geometric criteria is parameterized using the root diameter of the stem, the taper rate, number of branches, position of branch nodes, diameter of the branches and their orientation angles. The element-locating-algorithm scans over the base volume mesh and evaluates if the position of the center of a tetrahedron is within the tree-skeleton volume defined by the geometric criteria. If the center of the 23 3. Methodology element is within the tree-skeleton region, then the element is marked as a wood- element and the corresponding material properties are assigned. The element-locating-algorithm using the geometric criteria works as follows. The position vectors of the center of each element with respect to the branch nodes and the root-center are found. The radial distance to the center of an element from a branch-axis or the main-stem-axis is calculated as the magnitude of cross-product between the position vectors and the unit direction vector of the branch and stem axes. The geometric criteria checks if the radial distances fall within the radius of the stem or any of the branches in order to pick the wood elements. The radius of the main stem or of the branches are defined as a function of height and taper rate. The parameters involved in the geometric criteria are schematically represented in the Figure 3.4. 3.2.2 Local Mesh Refinement Along Stem & Branches In order to accurately resolve the stem and branches of the tree, as well as to ensure computational efficiency, local refinement of the input base-mesh along the tree- skeleton is required. This is achieved by incorporating mesh-refinement routines into the pre-processing steps. The previously defined element-locating-algorithm is modified to work as the definition of a refinement zone surrounding the stem and the branches. The local mesh-refinement level is set using a threshold for the maximum edge length of a tetrahedron element. For those elements which fall into the refinement zone and have element edges longer than the given threshold, the corresponding elements are split into two by adding a new vertex at the midpoint of the longest edge. This also splits all other elements sharing the edge. The splitting is repeated until the threshold element size is reached. The appropriate threshold is found from the mesh study. The wood elements marked using the element-locating-algorithm within the single- object tree-mesh, post local refinement is shown in Figure 3.5. 3.2.3 Partitioning of Tree Crown To introduce partitions into the crown, the crown-element-set is further divided into five regions (sets) surrounding the main stem and four branches. The elements are sorted in to these five regions based on their proximity to the nearest branch or the main stem. The radial distances calculated in the geometric criteria are used for this purpose. Based on the evaluated minimum of radial distances calculated between center of the elements and axes of the stem and branches, the respective region where the element belong is identified. Accordingly, region ids are also assigned to every crown-element in the base mesh. Further, routines are added to identify neighbor elements with dissimilar region ids and thus the location of partitions. Along the shared faces of the element pairs identified for partition, shared nodes are initially duplicated. To cut the topology along these faces, remapping of the nodes are done so that the involved element with the largest node index is assigned to use 24 3. Methodology Figure 3.5: Wood elements (marked in red) resolved within the volume mesh by the element-locating-algorithm following local refinement of input mesh using a threshold maximum edge length of 3 cm. (right) Zoomed in view of resolved main- stem showing the resultant rough surface Figure 3.6: Partitioned crown elements in PT. (a) Viewed in YZ plane (b) Viewed in XZ plane the duplicated nodes. Even though the duplicated and previously shared nodes are geometrically coincident, the remapping cuts the topology between elements. Flags are added to the partitioning routines to make sure that the partitions doesn’t cut through the regions of wood elements. By partitioning, the continuous spherical surface of the crown is modified to include the surface boundaries along the partitions. No contact constraints are defined along these new boundaries. Due to which,the regions are allowed to overlap into each other during deformation of the crown. Section views of the partitioned crown 25 3. Methodology surrounding the four branches are shown in Figure 3.6. 3.2.4 Mesh Connectivity Check Prior to the grid independence check of the single-object tree meshes, a primary mesh performance check was also done. For this, a comparison between the aggre- gate of extracted wood element volume and the reference volume of tree-skeleton is made. For coarse meshes, the extracted wood-element volume differed significantly from the reference volume. Also the entire length of the stem inside the crown was not resolved since the threshold element size was bigger than the cross-section diameter at a particular height. This caused the elements to fail the geometric criteria. With increasing local refinement level, the percentage of extracted length and vol- ume of the stem improved. For a particular threshold size, entire length of the stem is resolved, however the extracted volume of stem was still lower than the reference volume. A close inspection of the extracted stem elements showed discontinuities in Figure 3.7: Resolved stem before and after refinement and connectivity check. (a) before refinement and connectivity check (b) after refinement and connectivity check the stem-mesh and the missing elements were the reason for a lower magnitude of extracted stem volume. At the upper section of the extracted stem where the cross-section diameter is small, the elements which were partially within the stem region, also failed the geometric criteria. The position of similar mesh discontinuities in the upper part of the stem varied randomly for different meshes. The discontinuities can result in occurrence of unrealistic stress and strain concentrations in the stem, and this must be avoided. In order to eliminate the presence of mesh discontinuities in the stem, a mesh con- nectivity check was added as a post local refinement step. The connectivity check utilized node positions of the extracted wood-elements to evaluate connectivity with its neighbor elements. Once a discontinuity is located, the extracted wood elements 26 3. Methodology above the point of discontinuity are removed from the wood-element set and is trans- ferred to the crown-element set. The re-definition of wood elements which lie above the point of discontinuity to crown-element effectively reduces the overall height and volume of the stem. Hence, it is important to have sufficient local refinement proportional to the resolution of tree-skeleton that is desired. In the study the required refinement level is identified from the mesh study. The resolved stem of the ST before and after mesh refinement are shown in Figure 3.7. Figure 3.8: Extracted tree-skeleton surface (marked in white) inserted into carte- sian mesh to perform local refinement. Viewed in X-Z plane. 3.2.5 Surface Mesh Extraction for IB-tree Local Refinement As discussed in the introduction, we are not relying on a CAD geometry for the tree-skeleton in the study for both solvers. Also, we want to investigate the use of a fully permeable tree within the IBOFlow solver. Hence an alternative way to locate the tree skeleton and spatially define the permeability properties for stem and crown in the IBOFlow solver is used. 27 3. Methodology Figure 3.9: PAD distribution in IB-tree defining wood cells with PAD = 50 m2·m−3 and crown cells with PAD = 0.1 m2·m−3. (a) ST PAD distribution (b) BT PAD distribution In the IBOFlow solver, local refinement of the tree-skeleton surface is guided by the surface mesh extracted from the FEM solver. This mesh is incorporated into the CFD computational domain as an immersed body with a free-boundary condition that remains invisible to the incoming flow. The invisible surface mesh defines the required level of local refinement, which is crucial for accurately resolving the tree-skeleton’s thickness distribution and minimizing numerical diffusion during the assignment of permeability properties to the wood cells. Figures 3.8 shows the extracted tree-skeleton surface used for local refinement. The resulting single-object ST and BT IB-trees with spatial PAD fields in IBOFlow solver are shown in Figure 3.9. 3.3 FEM Simulations in LaStFEM Within the FEM Structural Solver (LaStFEM), our objective is to simulate the me- chanical response of the tree under both static and aerodynamic loading conditions. More broadly, we aim to capture the structural behavior of real trees subjected to wind by examining the parameters that define the tree models. A critical step in achieving this is the identification of suitable material properties for the artificial crown material used in the tree-models. 28 3. Methodology 3.3.1 FEM Simulation Setup The basic setup used for the FEM simulations, investigating artificial material pa- rameters for the tree-crown in the study are discussed in this section. 3.3.1.1 Material Properties for Artificial Crown Material Based on the small-strain elasticity model used in the study, it is necessary to define the stiffness, density, and Poisson’s ratio for the artificial crown material. Accordingly, the stiffness and density of the crown volume are specified relative to the material properties of the stem using stiffness and density ratios. The appropriate Poisson’s ratio, along with suitable stiffness and density ratios for the crown material, are determined through detailed parameter studies. Table 3.1: Material properties of Aspen wood used in the simulations. Material Property Value Young’s Modulus 11 × 109 Pa Density 440 kg m−3 Poisson’s Ratio 0.44 The material properties of Aspen wood taken from the reference article [12] are listed in the Table 3.1. The stiffness and density ratios for the artificial crown material are defined as follows: stiffness_ratio = Ecrown/Estem (3.1) density_ratio = ρcrown/ρstem (3.2) 3.3.1.2 FEM Boundary Conditions For the FEM parameter studies, the aerodynamic loads on the tree under normal wind conditions were approximated as a uniform pressure force acting on the sur- face nodes over the windward face (surface facing the wind direction) of the tree. Assuming a steady wind speed of 10 m/s and a drag coefficient of 0.5, the average pressure on the crown surface is estimated as half of dynamic pressure equal to 30 Pa. This is an overestimate of actual pressure for a wind speed of 10 m/s, since we neglect the reduction in drag coefficient due to through flow. The root plate was defined as a fixed support. 3.3.2 FEM Mesh Study Prior to parameter studies, the appropriate level of refinement required by the stem and crown regions of the single-object mesh must be identified. This is done using a mesh study for ST and BT geometries. The base mesh is parameterized using two variables. These are: • maximum element size for crown elements 29 3. Methodology • maximum element size for wood elements The size of the surface element is set as the minimum diameter limiter of the cross- section of the stem (0.025 m) and the volumetric growth rate is kept as 20%. Keeping these two parameters fixed, the maximum element sizes for crown and wood elements were varied to obtain different meshes for comparison. 3.3.2.1 Mesh Study Using ST Two set of meshes were evaluated with ST geometry. In the initial set, for a base mesh with fixed and uniform crown element size of 0.1 m, the maximum element size of wood elements were iterated. The previously defined local refinement operation was used to set the threshold maximum wood element size. Four meshes with threshold size as 0.1 m, 0.05 m, 0.03 m, and 0.02 m were generated. Table 3.2: Initial set of meshes comparing threshold maximum element size for wood elements Mesh Element Count (Million) Max Edge Length Stem (m) Extracted Stem Vol (%) Extracted Crown Vol (%) Bending Moment at Root (kNm) Effective Stem Height (m) M1 0.53 0.10 93.55 100.59 -23.85 18.94 M2 1.18 0.05 94.07 100.59 -24.88 21.47 M3 2.89 0.03 94.08 100.59 -24.88 21.96 M4 8.34 0.02 94.09 100.59 -24.89 22.33 (a) Element count vs max edge length (b) Bending moment vs max edge length Figure 3.10: Comparison of meshes for different max element size for stem ele- ments. Since the effective length of stem may vary with refinement level, the displacement over the crown can also vary significantly. Hence to monitor the convergence of mesh properties, the bending moment of the tree (effective load on the tree) measured at its root is used. 30 3. Methodology Table 3.3: Second set of meshes comparing maximum element size of the crown elements for max edge length of stem equal to 0.05 m Mesh Element Count (Million) Max Edge Length Crown (m) Extracted Stem Vol (%) Extracted Crown Vol (%) Bending Moment at Root (kNm) Stem Tip Displace- ment (m) M5 0.70 0.20 94.08 100.59 -27.619 1.562 M6 1.57 0.15 99.39 100.00 -27.651 1.300 M7 2.92 0.10 99.39 100.00 -27.651 1.300 Variation in bending moment seems to converge after mesh M2. And visual com- parison between M2 and M3 also showed that the thickness distribution of effective stem identified by the geometric criteria seems to follow the reference stem. The (a) Element count vs max edge length. (b) Tip displacement vs max edge length. Figure 3.11: Comparison of meshes for different max element size for crown ele- ments. comparison of effective extracted volume and actual tree volume until the point of discontinuity is also reported in the table showing initial mesh study results Table 3.2. From the initial mesh study results, mesh M2 with max edge length 0.05 m was selected for further studies. In the second set of meshes, the maximum element size possible for the crown ele- ments is determined. The converged maximum wood element size from the initial set was used in these cases. Three meshes with maximum crown element sizes 0.2 m, 0.15 m and 0.1 m were compared. The tip displacements of the effective stem were monitored for mesh convergence. The results from the second set are given in Table 3.3. Based on the convergence in tip displacement and root bending moment, mesh M6 with max edge length for crown elements equal to 0.15 m and max edge length for stem element equal to 0.05 m was selected as the final mesh for FEM parameter study. 31 3. Methodology Table 3.4: Second set of meshes comparing maximum element size of the crown elements using BT for max edge length stem equal to 0.05 m Mesh Element Count (Million) Max Edge Length Crown (m) Extracted Stem Vol (%) Extracted Crown Vol (%) Bending Moment at Root (kNm) Stem Tip Displace- ment (m) M5 0.83 0.20 95.50 100.59 -27.612 1.562 M6 1.70 0.15 100.80 100.00 -27.646 1.300 M7 3.03 0.10 100.80 100.00 -27.649 1.299 3.3.2.2 Mesh Study Using BT Using the defined refinement routines and mesh connectivity checks, the two set of cases were re-evaluated for the BT case. The results underlines the convergence of the selected mesh parameters as in M6. The results for BT are shown in Table 3.4. 3.3.3 FEM Parameter Study To re-iterate, the two-fold objective of the FEM parameter studies are: 1. Identify the relevant material parameters for the artificial crown material. 2. Find the bounds for these material parameters that can generate the expected physical character of a real tree. Throughout the parameter study simulations, the ST geometry is used and material properties of the wood elements are kept fixed. 3.3.4 Results & Discussions from FEM Parameter Study 3.3.4.1 Stiffness Ratio The primary response of the tree crown under wind load, happens as the deformation of the outer region of the crown. The local deformations in the crown surface results in overall streamlining of it. The stiffness properties of the simplified spherical crown should be in a range that facilitates this streamlining effect. Therefore, identifying appropriate stiffness ratio for the crown elements was the initial objective in the FEM parameter studies. A uniform stiffness ratio was used across the crown volume initially. The order of magnitude of the stiffness ratio was iterated between 10−1 and 10−8, while keeping density ratio and Poisson’s ratio fixed at magnitudes 10−8 and 0.44 respectively. The choice of density ratio to be 10−8 was to neglect the contribution of weight induced deformation into the streamlining. Also Poisson’s ratio was set to match the Poisson’s ratio of wood to keep the crown nearly incompressible during streamlining. The Poisson’s ratio was not set to 0.5 to avoid singularity of having an infinite 32 3. Methodology Figure 3.12: Deformation of crown viewed in X − Y plane passing through center of the crown. (a) stiffness ratio = 10−6 (b) stiffness ratio = 10−7 (c) stiffness ratio = 10−8 Figure 3.13: Equivalent strain plots of crown viewed in X − Y plane passing through center of the crown. (a) stiffness ratio = 10−6 (b) stiffness ratio = 10−7 (c) stiffness ratio = 10−8 bulk modulus. The observed deformation across the crown surface as well as the equivalent strain fields are shown in Figures 3.12 and 3.13 respectively. With decreasing stiffness ratio and effective stiffness of crown elements, local de- formations started to appear on the crown surface. For stiffness ratio 10−7 the streamlining effect starts to be evident and is significant for 10−8 case. Thus its inferred that the stiffness ratio has to be lowered to the order of 10−8, in order for the crown surface to show significant deformation under normal wind conditions. It is also important to remark that, with increasing resolution of the stem where we include the higher order branches, the resulting crown will be composed more of foliage and voids. The change in composition of the crown must be considered in the choice of stiffness ratio and other properties for the BT. Based on the equivalent strain plots in Figures 3.13 and 3.14, very high strain lev- els are observed in crown elements near the stem. This is attributed to the high accumulated pressure acting on these regions and must be mitigated. To reduce overall strain levels in the crown while still capturing the deformation of the crown 33 3. Methodology Figure 3.14: Equivalent strain plot of crown section viewed in X-Z plane passing through center of the crown for stiffness ratio = 10−8, density ratio = 10−8 and Poisson ratio = 0.44 surface, a non-uniform stiffness ratio is required. This is implemented through a spatially varying (field) definition of the stiffness ratio, which will be discussed in the following section. 3.3.4.2 Stiffness Ratio Field Figure 3.15: Characteristic stiffness of elements across the discretized crown vol- ume consisting of 50 levels (a) ST stiffness field (b) BT stiffness field Like we discussed in the beginning of section 3.3.4.1, compared to the deformation of the crown elements along the crown surface, crown volume close to the stem are less deformed. It is observed in a real tree that the crown volume close to the stem 34 3. Methodology and at the core of the crown, follow the stem during its displacement rather than undergoing significant deformation. Hence it is reasonable to assume that the crown elements close to the stem are stiffer than those at the surface in order to be able to resist deformation. The degree of deformation in crown elements is observed to increase as we move radially outward across the crown. Accordingly, a spatial stiffness ratio field was defined to assign stiffness to crown elements based on their distance from the nearest stem or branch. The crown volume was discretized into multiple levels, with the zeroth level being closest to a stem or branch, and each level assigned with a unique id. The stiffness ratio was then varied across these levels according to the minimum of evaluated distances to the stem or branch axes. (The specific discretization approach was because of a current limitation in LaStFEM restricting element-wise property definition.) The growth of the stiffness ratio from its value at the crown surface to the value at the zeroth level was defined to follow an inverse relationship with the radial distance from the stem axis. Assuming a cylindrical iso-surface around the central stem axis, internal stresses were approximated as pressure applied over the surface of this cylinder. Under this assumption, the internal stress increases inversely with radial distance from the stem axis. To counterbalance this increase, the stiffness ratio is scaled proportionally from its base value at the crown surface. For cylindrical iso-surface assumption: Piso = Faero 2πrisol (3.3) For spherical iso-surface assumption: Piso = Faero 4πr2 iso (3.4) For more complex tree geometries with extensive branching, a spherical iso-surface approximation can be more appropriate than a cylindrical one. In such cases, the growth of the stiffness ratio can be considered to follow an inverse-square relationship with radial distance. The gradual decrease in characteristic stiffness of the crown elements with respect to radius, helps to avoid high strain levels across the volume of the crown and to localize large deformations along the outer regions of the crown. The equivalent strain in crown elements across the different levels achieved using the field definition of stiffness ratio are shown in Figure 3.16. In the plots, moderately high but localized strain is observed at elements close to the tip of stem and branches. However the average equivalent strain values across the volume range between 5 to 10 percent. By defining stiffness ratio as function of the proximity to nearest stem or branch, a more robust distribution of stiffness is achieved in the crown. Depending on the composition of the crown, the order of decay as well as magnitude of stiffness ratio at the surface of the crown can be easily modified with this approach. Hence instead of uniform stiffness ratio definition, hereafter the simulations were carried out using stiffness fields. 35 3. Methodology Figure 3.16: Equivalent strain plots in deformed state of elements across the discretized crown volume consisting of 50 levels (a) ST stiffness field (b) BT stiffness field 3.3.4.3 Poisson Ratio The second key parameter affecting the overall deformation of the crown is the Poisson’s ratio. In the set of cases discussed in section 3.3.4.1, we used a Poisson ratio of 0.44. This makes the crown nearly incompressible for wind loads and this is not the case for an actual tree. The morphology of a tree crown allow it to undergo compression to certain limit. The presence of voids and the independent nature of motion of the branches, makes this compression to occur without lateral expansion or contraction. Resolving this behavior of the crown was the primary motivation to investigate the influence of Poisson’s ratio in crown deformation. Figure 3.17: Equivalent strain plots of crown viewed in X − Y plane passing through center of the crown. (a) Poisson’s ratio = 0.44 (b) Poisson’s ratio = 0.0 (c) Poisson’s ratio = -0.2. Keeping the stiffness ratio and density ratio at 10−8, the Poisson’s ratio of the crown elements was iterated between 0.44 and -0.2. An element with Poisson ratio close 36 3. Methodology to 0.5, behaves like a rubber while a negative Poisson ratio would mean lateral contraction under compression (expansion under tension). However Poisson ratio of magnitude 0.0 allows compression without lateral expansion or contraction. This is similar to response observed in porous materials as well as in certain honeycomb structures [25, 7]. The equivalent strain plot of the deformed crown cross-section for Poisson’s ratio 0.44,0.0,-0.2 are shown in Figure 3.17. For decreasing Poisson’s ratio, the lateral expansion (bulging) of the crown decreased. And for Poisson’s ratios 0.0 and below, some amount of compression followed the streamlining deformation. The area of cross-section of the crown is also observed to be smaller in the case of zero Poisson’s ratio than in the case of 0.44. Between Poisson’s ratio 0.0 and negative values, variation in deformation was minimal. This might be because of the nature of the loading which is uniform across the crown surface. The Poisson’s ratio 0.0 was selected for future cases since it aligned with the expected behavior in actual trees. Even though the average strain levels in the crown is higher than accepted strain levels for small strain elasticity model, the choice of Poisson’s ratio zero for the crown allows to decouple the strains and this can improve stability of the FE solver. The accuracy of resolving local deformation in the crown using small strain elasticity model has to be investigated. However, the objectives of streamlining of the crown is achieved with the selected Poisson’s ratio and stiffness distribution. This will ensure realistic variation of aerodynamic load in flow simulations. Since the strain levels in the stem are within the limit of 0.05%, the primary bending response of the tree is accurately resolved using the selected material models. 3.3.4.4 Density Ratio To accurately predict the overall response of the crown, a realistic representation of crown mass is essential. After establishing the stiffness and Poisson’s ratio for the crown elements, a homogenized bulk density for the crown was initially estimated. Average total mass and foliage mass fraction (leaves and twigs) for an aspen tree were taken from reference [17]. For an average total fresh weight of 636 kg, the fo- liage mass fraction was taken as 0.3, giving a crown mass of approximately 190.8 kg. Given a crown volume of 268 m3 (from the ST crown geometry), this corresponds to a bulk density of 0.712 kg/m3 and a resulting density ratio of 0.00162. For the ST model, the volume of wood is approximately 0.96 m3. Using a reference wood density of 440 kg/m3 for aspen, the total wood mass sums to 613.2 kg. This closely matches the average tree mass referenced in the study, supporting the accuracy of the estimate. Since the density ratio is derived solely from the foliage mass fraction, the value of 0.00162 is used as a reference for the parameter study. Using this estimated bulk density and corresponding density ratio, the crown’s deformation by self-weight was simulated. The resulting deformation is shown in Figure 3.18. The results reveal significant deformation due to self-weight. The crown, initially spherical, deforms into an oval shape, and the stem and branches exhibit noticeable bending. This contradicts the assumption that, in the absence of wind, the crown 37 3. Methodology Figure 3.18: Deformation under self-weight of the crown viewed in the X–Z plane for the estimated density ratio of 0.00162, (stiffness ratio = 10−8, Poisson ratio = 0.0). maintains a spherical equilibrium geometry. The approximated tree skeleton is un- able to support the estimated crown mass under gravity alone. It was therefore concluded that a realistic density ratio must incorporate actual thickness variations across the tree skeleton. To preserve the spherical no-wind ge- ometry and ensure that all deformation arises solely from aerodynamic loading, the crown density ratio was excluded from the simulation. Instead, a very low density ratio of 10−8 was assigned to the crown, effectively rendering it massless. To still account for the crown’s mass and its damping contribution, the estimated crown mass was redistributed across the stem and branches. This assumes a bal- anced mass distribution, justified by the hypothesis that trees optimize their struc- tural mass distribution during growth [26]. 3.3.5 Base Case Deformations in ST, BT & PT To conclude the parameter studies in LaStFEM, pure FEM simulations were per- formed on ST, BT, and PT using a consistent set of established parameters. The common parameters applied across all three cases include: • Stiffness ratio = 10−8 • Poisson ratio = 0.0 • Density ratio = 10−8 A comparison of the deformation plots reveals that, relative to ST, the crown in BT undergoes less compression. This indicates that the added structural stiffness provided by the branches plays a significant role in resisting deformation. Within the PT, an asymmetric deformation is observed. Compared to the windward partition zone of the crown, the lateral and top partitions experience greater dis- 38 3. Methodology placement due to pressure loading. The windward zone undergoes less deformation, primarily because of the reduced exposed area and the corresponding decrease in bending moment, which together diminish the effective wind load acting on that region. Also as expected, the partitions are facilitating overlap of the crown regions. Across all three cases, the crown exhibits streamlining accompanied by primary bending deformation. By incorporating the surrounding pressure field from CFD simulations, the one-way coupled FSI simulations will provide a more realistic rep- resentation of the crown’s deformation behavior. (a) ST (b) BT (c) PT Figure 3.19: Pure FEM base case deformations in ST, BT, and PT, viewed in the X–Z plane. 3.4 CFD Simulations in IBOFlow In the IBOFlow solver, the flow field over and through the tree is computed. From the flow field, the aerodynamic loads acting on the tree can be found. Similar to the single-object tree representation in structural solver, a single IB-object tree with internal cells having varying permeability properties across wood and crown regions is defined in the flow solver. This is in contrast to the approach where a permeable fluid-object crown is spatially coupled to a solid-object stem. Instead of fully resolving the flow through the tree crown, permeability model based on source term modification are used to induce variation in flow properties. In a 39 3. Methodology (a) ST (b) BT (c) PT Figure 3.20: Pure FEM base case deformations in ST, BT, and PT, viewed in the X–Y plane. similar way, the source terms are modified in the wood cells in order to obstruct the flow through them. Using a relatively higher value of PAD than the realistic PAD values used in the crown, the incoming flow is obstructed in the wood cells. Appropriate PAD values for both wood and crown regions are identified through parameter studies. 3.4.1 CFD Simulation Setup Figure 3.21: Overall dimensions of computational domain in CFD simulations with boundary conditions used. The basic setup used for CFD simulations in the IBOFlow solver are discussed in this section. The computational domain is defined according to the best practice guidelines discussed by Blocken et al.[3]. The inlet conditions follow the definition of standard ABL flow [20, 42]. The ground is defined as a rough wall with aerodynamic roughness height chosen as 0.0018 m (flat terrain with no obstacles). The average wind speed at the crown level is approximately 12 m/s, according to the generated inlet velocity profile. The sides and top boundaries of the domain are defined using 40 3. Methodology symmetry boundary conditions. At the outlet, a pressure outlet condition is used. The overall dimensions of computational domain is given in figure 3.21. The log-law inlet velocity profile, inlet profiles for turbulent kinetic energy and dis- sipation rate generated according to the standard ABL flow wall functions discussed in 2.3.4 are given in Figure 3.22. Figure 3.22: Velocity, turbulent kinetic energy, dissipation rate inlet profiles gen- erated according to standard ABL wall functions.(z0 = 0.0018, u∗ ABL = 0.73 m/s, κ = 0.42). To model the through flow in crown, the permeability model discussed by Sogachev et al. [34], with coefficients proposed for a solitary tree by Dellwik et al. [8] are used. Based on the measurements of drag coefficient by Dellwik et al, the drag coefficient in the source term for crown cells was set to be fixed and equal to 0.25. For the wood cells, drag coefficient was set equal to one. This makes the source terms to be a function of PAD alone. Previous validation studies (not part of this work) have shown good performance of the model working with k-ϵ turbulence model. Hence the simulations in the study are conducted using it. 3.4.2 CFD Mesh Study For the mesh study in IBOFlow, two set of cases were conducted. The initial set of meshes focused on identifying the global mesh parameters for the exterior cells of the IB-tree as well as getting a reference for the forces over tree stem. The mesh parameters included global refinement levels and local wake refinements for the tree-object. To this purpose, a benchmark case using solid object stem of ST was conducted. For different mesh refinement levels, the overall pressure force, calculated by integrating over the surface of the IB-tree-stem were monitored. The force comparison for different mesh levels from the initial set of cases are shown in Table 3.5. 41 3. Methodology (a) Velocity field for ST solid stem reference case (b) Pressure field for ST solid stem reference case Figure 3.23: Velocity contour plots for the benchmark cases using ST solid stem and global mesh parameters of Mesh M3 The velocity and pressure fields plotted along the X−Z plane passing through center of the solid stem for ST are shown in Figure 3.23. From the plots, a peculiar pattern is observed in the wake of the tree. This pattern resembles the 2D snapshot of the wake field, for a tapered cylindrical body kept in cross flow [21, 37]. The reducing amplitude across the height corresponds to the reducing characteristic length scale (cross-section diameter) of vortices shed by the cylinder. The overall pressure force calculated over the stem seems to converge from M2. Also, the global mesh parameters and the level of wake refinement used for mesh M2 seems to resolve the mean wake field with sufficient accuracy and therefore, it is followed in further cases. 42 3. Methodology Table 3.5: Mesh study using solid-object stem of ST Mesh Wake Cell Size (m) Streamwise Force Com- ponent ST (kN) M1 0.031 0.4959 M2 0.016 0.4731 M3 0.008 0.4737 (a) Velocity field for BT using mesh M4 (b) Velocity field for BT using mesh M5 (c) Velocity field for BT using mesh M6 Figure 3.24: Velocity contour plots for different mesh resolutions of IB-tree with PAD for wood cells = 50 m2·m−3 and PAD of crown cells equal to 0.1 m2·m−3. The second set of meshes were aimed at identifying the refinement levels required for the internal cells of the IB-tree as wells for resolving the branches of BT. Here actual PAD based fully permeable IB-tree with PAD of the stem equal to 50 m2· m−3 and PAD of crown cells equal to 0.1 m2· m−3 is used. Three meshes with varying 43 3. Methodology (a) Pressure field for BT using mesh M4 (b) Pressure field for BT using mesh M5 (c) Pressure field for BT using mesh M6 Figure 3.25: Pressure contour plots different mesh resolutions of IB-tree with PAD for wood cells = 50 m2·m−3 and PAD of crown cells equal to 0.1 m2· m−3 local refinement level for the tree-skeleton surface was generated and compared. The overall pressure force acting on the IB-tree surface measured for these three meshes are shown in Table 3.6. The corresponding velocity and pressure fields are shown in Figures 3.24 and 3.25. Analyzing the result plots, wake field of the main-stem and the branches are well resolved with mesh M6. This is evident with a comparison of wake near lower part of the stem as well as directly behind the branches. The forces obtained for the ST stem indicate the onset of convergence starting from mesh M5. Therefore mesh refinement properties of M6 is followed in CFD parameter study meshes. 44 3. Methodology Table 3.6: Mesh study using fluid-object IB-tree of ST with PAD-wood = 50 m2·m−3, PAD-crown = 0.1 m2·m−3 Mesh Internal Cell Size (mm) Streamwise Force Component ST (kN) Streamwise Force Component BT (kN) M4 1.95 1.157 1.739 M5 0.97 1.220 1.895 M6 0.65 1.250 2.026 3.4.3 CFD Parameter Study As discussed before, the objective in CFD parameter studies is to identify the ap- propriate level of permeability defined by PAD, for wood and crown cells. Thus the key parameters investigated in IBOFlow parameter study simulations are, • PAD for wood cells (PAD-wood) • PAD for crown cells (PAD-crown) The PAD for the wood cells should allow to mimic the obstruction of flow caused by a solid stem while PAD for crown cells should resolve the through flow properties in the crown. Therefore, two fold simulations are conducted here. 3.4.4 Results & Discussions from CFD Parameter Study 3.4.4.1 PAD for Wood Cells Table 3.7: Results from parameter study of PAD value for wood cells using ST stem PAD-wood (m2·m−3) Velocity Inside Stem (m/s) Streamwise Force Component (kN) 10 4.86 0.326 30 2.97 0.422 50 2.30 0.459 70 1.94 0.478 90 1.54 0.470 In the first set of cases, Plant Area Density of Wood Cell (PAD-wood) is established. For this, PAD value for the crown cells is set to zero. Using ST fluid-object stem and keeping the drag coefficient equal to one, PAD-wood in the permeability model was iterated in the cases. For different PAD values, the overall pressure force acting over the permeable ST-stem was monitored. The iterated PAD-wood values and the corresponding forces on the permeable stem are presented in Table 3.7. The variation of forces and velocity measured inside the stem at a distance of 1 m from the ground are plotted in the Figure 3.26. 45 3. Methodology (a) Through-flow velocity variation (b) Pressure force component variation Figure 3.26: Through-flow velocity and streamwise pressure force component variations from PAD-wood parameter study results Comparing the results for solid stem and the fluid-object permeable stem, permeable stem does have some amount of through flow in it. The velocity of this through flow is decreasing with increasing PAD value. Comparing the forces obtained for solid ST-stem with permeable ST-stem, for PAD- wood values equal to 50, 70 and 90 m2·m−3, the force computed over permeable stem are approximately 97.0%, 101.1% and 99.4% respectively of the force over solid-stem. This indicates that the approach of using high PAD values to mimic a solid stem is effective in reproducing the forces with sufficient accuracy. Although the pressure force converges with increasing PAD, the slightly higher force magnitudes compared to the expected maximum of the solid stem may be attributed to a more uniform pressure field around the permeable stem. Unlike the solid stem, which exhibits a wavy pattern in the wake due to vortex shedding, the permeable stem produces a more stable and symmetric flow. This can result in a higher average pressure in the wake region, thereby contributing to the increased overall pressure force. A second measure of correspondence between the solid stem and permeable stem is done by comparison of the wake fields. The wake velocity profile plotted for the solid stem and permeable stem along different stations in the wake of the tree is shown in Figure 3.27. The root radius of the stem is 0.2 m and radius of the crown is 4 m. The probing stations 1,2,3 and 4 are at locations X = 0.25 m, 0.5 m, 5.0 m and 12.0 m respectively, within the X-Z plane passing through center of the stem. The Mean Relative Error (MRE) of respective profiles with respect to the solid-stem profile are also shown in the Figure 3.27. If we observe the profiles across the four stations, we can see that the profiles closely match over the length of the stem. The maximum deviation of approximately 5% is observed at the station - 1 (X = 0.25 m). Further downstream at all the sta- tions MRE is below 2%. Beyond the radial extent of the crown, the variations are insignificant. A similar comparison of wake turbulent kinetic energy profiles are shown in Figure 46 3. Methodology (a) Station - 1 (X = 0.25 m) (b) Station - 2 (X = 0.5 m) (c) Station - 3 (X = 5 m) (d) Station - 4 (X = 12 m) Figure 3.27: Wake velocity profiles for solid stem and PAD based permeable stem (PAD-wood = 50 m2·m−3 & 90 m2·m−3) plotted along stations downstream the stem.(Plots are based on mesh M4) 3.28. In contrast to the velocity profiles, the turbulent kinetic energy profiles show significant variation from the solid reference profile, at stations 1 & 2. This suggests that the source term for turbulent dissipation rate requires further tuning. With high PAD-wood values, the turbulent dissipation rate is scaled significantly across the volume of the stem. With increasing dissipation rate, it was expected that an under-esti