Simulating Trees

Why Simulate Trees?

Trees not only convert carbon dioxide into oxygen and eliminate pollution and airborne toxins, but they also significantly influence wind patterns as it passes through their leaves. Trees play a crucial role in enhancing pedestrian wind comfort in urban areas. Here are the key points:

  1. Windbreak Effect: Trees act as natural windbreaks, reducing wind speeds at ground level. Even a single tree can make a noticeable difference, creating a more comfortable environment for pedestrians. Removing all trees from an area can double the wind speeds, significantly decreasing pedestrian comfort.
  2. Canopy Density and Shape: The shape and density of tree canopies are important. Denser canopies with lower porosity provide better wind sheltering.
  3. Aspect Ratio: The aspect ratio (height to width ratio) of a tree canopy affects the size of the sheltered zone behind the tree and the distance required for wind speeds to recover.
  4. Bare Branches in Winter: Even the bare branches of deciduous trees in winter can moderate airflow and reduce wind pressure on buildings.
  5. Sheltering Effect Distance: The sheltering effect of trees extends up to around seven times their height downwind, with the most significant reduction in wind speed occurring within 0.5 times the tree height.

Urban planners and designers should consider tree planting and species selection to maximize pedestrian wind comfort. Strategic placement of trees can serve as effective windbreaks, significantly enhancing outdoor spaces.

Different tree species, characterized by varying shapes, sizes, and leaf types, affect wind velocity to different extents. While wind is not completely stopped by trees, its speed is reduced as it moves through their branches. The key factor here is leaf density: trees with denser foliage offer greater resistance to airflow. For instance, a Silver Birch tree has less air resistance compared to a Chestnut tree due to its lower leaf density.

This ability of trees to allow air to pass through is known as porosity and is quantified using a leaf area index (LAI). The LAI values are determined through experiments.

When trees with high leaf area index values—indicating greater wind resistance—are strategically placed in areas with fast airflow and strong winds, they can significantly improve wind comfort and safety for pedestrians. These trees act as effective windbreaks, creating more pleasant and secure outdoor environments. For more detailed explanations on this topic, refer to the upcoming webinar.

How Do You Simulate a Tree?

Introduction

In Computational Fluid Dynamics (CFD), trees can be modeled as porous media to accurately represent their impact on airflow. This approach considers trees as objects that allow air to pass through while imposing resistance, thus affecting wind speed and turbulence. Modeling trees as porous media is essential for understanding environmental phenomena like windbreaks, pollution dispersion, and microclimate regulation.

Concept of Porous Media

Porous media are materials containing pores (voids). The flow through porous media is characterized by the interaction between the fluid and the solid matrix. In the context of trees:

  • The foliage and branches constitute the solid matrix.
  • The spaces between leaves and branches act as the pores through which air flows.

Key Parameters in Porous Media Modeling

To represent the plant canopy regions, the user must specify the C_d and LAD values.

  1. Leaf Area Density (LAD):
    • Definition: LAD is the leaf area per unit volume of foliage.
    • Role: Determines the porosity of the tree and influences how much air resistance the foliage provides.

LAD is a local characteristic, it is common to make the hypothesis of horizontal homogeneity. However, LAD will depend on the vertical location LAD(z). The empirical relation allows us to link LAD(z) to a global quantity LAI. The values for LAI can be found in the literature.

LAI Values for Woody Plants

$$ LAD = LAI * \frac{h_t - z_m}{h_t - z}^n \exp {n{ 1 - \frac{h_t - z_m}{h_t - z}}} $$

Here are some of the values for the most common species of trees in Europe:

Tree SpeciesLAI
Sycamore2.97
Silver Birch3.24
Quercus5.16
Chestnut5.19
Platanus5.28
  1. Drag Coefficient C_d:
    • Definition: A dimensionless number representing the drag force exerted by the leaves and branches.
    • Role: Used in momentum equations to calculate the resistance offered by the tree to the airflow.

Equations for Modeling Trees as Porous Media

We base our tree model on the OpenFoam atmPlantCanopy sources. To apply the atmPlantCanopy sources, the user must specify the C_d and LAD values in the constant/fvOptions file for the desired cell zones representing the plant canopy regions.

We provide the relation that supposes the foliage density is independent of the vertical direction. This is a simplification to generalize the LAD at the volume of the tree with an average value:

$$ LAD = \frac{LAI}{h} $$

  1. Momentum Sink Term:

    Represents the loss of momentum due to the drag force of the tree.

    $$ S_p = - \alpha \rho C_d LAD |u_0| u $$

    Where:

    • S_p is the source term applied to the momentum equation
    • α is the phase fraction (1 for single-phase flows)
    • ρ is the fluid density
    • C_d is the plant canopy drag coefficient (a user-specified value)
    • LAD is the leaf area density (a user-specified value representing the plant canopy density)
    • |u_0| is the magnitude of the velocity field from the previous iteration
    • u is the current velocity field

The source term S_p acts as a drag force opposing the flow velocity u within the plant canopy region. The magnitude of this drag force is proportional to the leaf area density LAD and the drag coefficient C_d, which characterize the plant canopy geometry and aerodynamic properties.

This source term is then incorporated into the momentum equation during the solution process, effectively reducing the velocity within the plant canopy regions to account for the drag induced by the vegetation.

  1. Turbulence Source Term:

For the turbulence kinetic energy dissipation rate (ε) equation:

$$ S_p = \alpha \rho \left( C_1 - C_2 \right) C_{\text{canopy}} ε $$

with:

$$ C_{canopy} = 12.0 \sqrt{C_\mu} C_d LAD |u_0| $$

Where:

  • S_p: Source term
  • C₁, C₂: Model constants
  • C_μ: Empirical model constant
  • ε: Turbulent kinetic energy dissipation rate
  • C_d: Drag coefficient
  • LAD: Leaf area density
  • |u₀|: Velocity magnitude (previous iteration)
  • α: Phase fraction
  • ρ: Fluid density

The source terms S_p account for the production and dissipation of turbulence within the plant canopy region. The magnitudes of these source terms are proportional to the leaf area density LAD and the drag coefficient C_d, which characterize the plant canopy geometry and aerodynamic properties.

These source terms are then incorporated into the turbulence transport equations during the solution process, effectively modifying the turbulence levels within the plant canopy regions to account for the effects of vegetation.

NUMERICAL IMPLEMENTATION-DYNAMIC GEOMETRIC ANALYSIS: VOXEL COLUMN PROFILING

The core of the methodology lies in a Voxel-Based Canopy Profiling (VCP) approach. This ensures that the Leaf Area Density (LAD) is calculated based on the effective biological thickness of the vegetation, regardless of its placement in the digital domain.

SLOPE-INVARIANT HEIGHT EXTRACTION

The VCP algorithm was implemented to address two critical geographical challenges:

  • Orography & Slopes: In presence of inclined terrain, a standard height calculation would include the elevation gain of the ground, artificially inflating the canopy volume.

  • Multi-Altitude Assets: When a single STL file contains vegetation at different ground levels (e.g., trees in a valley and trees on a hill), a global evaluation fails to provide a consistent physical density.

To resolve this, the algorithm discretizes the STL mesh into a 2D horizontal grid:

  • Voxelization & Resolution: The mesh is discretized using a horizontal grid with a spatial resolution of 0.1m (10cm). This dimension is the optimal balance between capturing canopy structural detail and maintaining computational efficiency for the CFD solver.

  • Peak-to-Peak (PTP) Analysis: For each grid cell, the local height Hi is computed as: $$ Hi= Zmax,i−Zmin,i $$ This method eliminates the “elevation bias,” ensuring the height represents only the canopy thickness and not the ground altitude.

  • Equivalent Height (Heq): The final LAD is derived from the mean of all local PTP heights, providing a consistent physical value even for “hollow” meshes exported from Blender.

NUMERICAL STABILITY AND THE K STABILITY FACTOR

To prevent numerical divergence—a common issue when high-density vegetation zones interact with complex urban wind fields—a dynamic stability limiter (K_stability) was introduced.

  • The Stability Criterion: The K factor limits the maximum permissible LAD based on the characteristic length scale of the vegetation volume (V1/3).

  • Technical Justification: This prevents the momentum source term from becoming disproportionately large relative to the cell size, which would otherwise lead to “Floating Point Exceptions” or unrealistic pressure gradients. The final applied LAD is defined as:

LADfinal = min(HrealLAI, V1/3K)

MATHEMATICAL JUSTIFICATION OF THE STABILITY CRITERION

The vegetation is represented in the Navier-Stokes equations by a source term S_p (momentum sink):

$$ S_p = - \alpha \rho C_d LAD |u_0| u $$

Where Cd is the drag coefficient (set to 0.2).

From a numerical standpoint, when the LAD value is excessively high, the source term becomes non- linearly dominant over the convective and diffusive terms. This leads to an “ill-conditioned” matrix in the linear solver. The K_stability factor acts as a diagonal dominance enhancer. By limiting the LAD based on the characteristic length V1/3, we ensure that the pressure drop across a single computational cell does not exceed the local kinetic energy, preventing numerical oscillations and “checkerboarding” effects in the pressure field.

SCIENTIFIC & REGULATORY COMPLIANCE

The methodology is strictly aligned with international standards and peer-reviewed literature to ensure both physical realism and legal defensibility of the results:

  • VDI 3787 Part 1 & Part 2 (Environmental Meteorology): These standards provide the fundamental framework for determining environmental parameters and controlling urban pollution/climate. Our methodology strictly adheres to VDI 3787 for the parameterization of porous media. This ensures that the momentum sink terms in OpenFOAM accurately reflect real-world plant-atmosphere interactions, emphasizing that for large-scale urban simulations, porous media must be parameterized to avoid “numerical blockage.”

  • COST Action 732 (Quality Assurance of Microscale Meteorological Models): These guidelines recommend that source terms for vegetation be carefully scaled with the grid resolution. By implementing the Voxel-based Canopy Profiling (VCP) approach, ArchiWind satisfies the requirement for grid-resolved canopy morphology. This prevents the numerical overestimation of drag forces typical of static height models and maintains a balanced Cell Peclet Number, ensuring the advection- diffusion stability of the solver.

  • Voxel-Based Volume Modeling (Shihua et al., 2017 & Wang et al.): The transition to a discretized height extraction follows the “Voxel Column” protocols. The algorithm utilizes a 0.1m (10cm) spatial resolution, which is the recognized gold standard in LiDAR Remote Sensing and GIS (ISO/TC 211). This specific resolution allows for the capture of complex 3D canopy structures while filtering sensor noise, as established by the USDA (United States Department of Agriculture) forestry inventory protocols.

  • Aerodynamic Modeling (Osama, 2006): The analytical derivation of the Leaf Area Density (LAD = LAI/H) is supported by the seminal work of Osama (2006) on wind flow simulation through urban greenery. This approach is the industry standard for translating individual tree geometries into porous zones within numerical solvers, ensuring that the resistance (momentum sink) is proportional to the foliage density.

  • Numerical Stability & The K-Factor: To prevent the modeled vegetation from acting as an unintended “numerical wall,” the K-stability factor incorporates the volumetric dimension (V1/3). This scales the resistance to the physical size of the object, maintaining a local momentum balance that prevents the source term from overwhelming the diffusive transport at the specific grid resolution.

IMPACT ON LOCAL TURBULENCE MODELING

By capping the LAD at a scientifically defensible threshold through K, we also ensure the validity of the k−ε (Standard or Realizable) turbulence model. Excessive LAD values can lead to unrealistic production of turbulent kinetic energy (k) and dissipation (ε) at the vegetation interface. The K_stability ensures that the SDR (Source-to-Diffusion Ratio) remains within the limits of the Boussinesq hypothesis, providing a more accurate estimation of the wake turbulence and wind comfort zones behind the trees.

SOLVER CONFIGURATION AND CONVERGENCE

Given the high turbulence and sharp gradients generated by vegetated zones, the SimpleFoam solver in Relaxed mode was utilized.

  • Under-Relaxation Factors: Optimized relaxation factors for pressure (0.3) and velocity (0.7) were applied to ensure smooth convergence.

  • Convergence Reliability: By integrating dynamic geometric correction with a stabilized solver configuration, the simulation achieved residual convergence below 10−4, providing a robust foundation for pedestrian comfort assessment and urban wind flow analysis.

FeatureLegacy Approach (Static)Optimized Approach (Dynamic)
Reference HeightFixed Default (10m)Mesh-Based (e.g., 50m)
LAD AccuracyOverstimated (up to 5x)Physically Consistent
Numeical StabilityHigh risk of divergenceControlled via K_stability
Solver Succes RateLow (Frequent FPE crashes)High (Consistent Convergence)

How to Upload Your Tree 3D Model?

ArchiWind allows adding vegetation under the Trees patch. Trees are a key agent for mitigating wind exposure locally. They will act as a porous medium and naturally dissipate excessive wind energy. The optimal placement for trees can be assessed through the comparison feature. Tree geometries are to be expressed as volumetric shapes. Adding trees with foliage represented as surfaces or aggregating tiny elements will not be recognized properly.

Valid Tree Example

Setting tree porosity

Dissipation of wind’s kinetic energy depends heavily on how porous the foliage is. To improve the accuracy of your simulation, you can specify the porosity level of the tree canopy. This ranges from Low Dense to Very Dense. The denser the foliage, the more it slows down the wind.

If you’re not sure which tree species will be planted, it’s recommended to choose the “Low Dense” option. This provides the most conservative estimate, ensuring that the simulation doesn’t unrealistically reduce wind speeds.

The table below links common tree species to their Leaf Area Index (LAI) and the porosity setting:

Tree SpeciesLAIPorosity
Eucalyptus1.0770% (Low Dense)
Sycamore2.9750% (Medium Dense)
Silver Birch3.2430% (Dense)
Platanus5.2820% (Very Dense)