WindShear.hpp

Go to the documentation of this file.

// Wind Shear Profiles
// Altitude-dependent wind profiles for boundary layer and low-level wind shear
#pragma once
 
#include <vulcan/wind/WindTypes.hpp>
 
namespace vulcan::wind_shear {
 
// ============================================================================
// Physical Constants
// ============================================================================

inline constexpr double VON_KARMAN_CONSTANT = 0.41;
 
// ============================================================================
// Linear Wind Shear
// ============================================================================

template <typename Scalar>
Scalar linear(const Scalar &altitude, double base_wind, double base_altitude,
              double shear_rate) {
    return base_wind + shear_rate * (altitude - base_altitude);
}

template <typename Scalar>
wind::WindVector<Scalar>
linear_vector(const Scalar &altitude, const wind::WindVector<double> &base_wind,
              double base_altitude, double shear_rate) {
    // Compute base wind speed (numeric)
    double base_speed = std::sqrt(base_wind.north * base_wind.north +
                                  base_wind.east * base_wind.east +
                                  base_wind.down * base_wind.down);
 
    // Avoid division by zero for calm conditions
    if (base_speed < 1e-10) {
        return wind::WindVector<Scalar>{.north = altitude * Scalar(0),
                                        .east = altitude * Scalar(0),
                                        .down = altitude * Scalar(0)};
    }
 
    // Scale factor based on altitude
    Scalar new_speed = linear(altitude, base_speed, base_altitude, shear_rate);
    Scalar scale = new_speed / base_speed;
 
    return wind::WindVector<Scalar>{.north = scale * base_wind.north,
                                    .east = scale * base_wind.east,
                                    .down = scale * base_wind.down};
}
 
// ============================================================================
// Power Law Wind Shear (Atmospheric Boundary Layer)
// ============================================================================

template <typename Scalar>
Scalar power_law(const Scalar &altitude, double ref_wind,
                 double ref_altitude = 10.0, double exponent = 1.0 / 7.0) {
    // V = V_ref * (h / h_ref)^alpha
    // Use janus::pow for symbolic compatibility
    return ref_wind * janus::pow(altitude / ref_altitude, exponent);
}

template <typename Scalar>
wind::WindVector<Scalar>
power_law_vector(const Scalar &altitude,
                 const wind::WindVector<double> &base_wind,
                 double ref_altitude = 10.0, double exponent = 1.0 / 7.0) {
    // Compute base wind speed (numeric)
    double base_speed = std::sqrt(base_wind.north * base_wind.north +
                                  base_wind.east * base_wind.east +
                                  base_wind.down * base_wind.down);
 
    if (base_speed < 1e-10) {
        return wind::WindVector<Scalar>{.north = altitude * Scalar(0),
                                        .east = altitude * Scalar(0),
                                        .down = altitude * Scalar(0)};
    }
 
    // Scale factor from power law
    Scalar scale = janus::pow(altitude / ref_altitude, exponent);
 
    return wind::WindVector<Scalar>{.north = scale * base_wind.north,
                                    .east = scale * base_wind.east,
                                    .down = scale * base_wind.down};
}
 
// ============================================================================
// Logarithmic Wind Profile
// ============================================================================

template <typename Scalar>
Scalar logarithmic(const Scalar &altitude, double friction_velocity,
                   double roughness_length, double displacement = 0.0) {
    // V = (u* / κ) * ln((h - d) / z_0)
    Scalar effective_height = altitude - displacement;
    return (friction_velocity / VON_KARMAN_CONSTANT) *
           janus::log(effective_height / roughness_length);
}

inline double friction_velocity_from_ref(double ref_wind, double ref_altitude,
                                         double roughness_length,
                                         double displacement = 0.0) {
    double effective_height = ref_altitude - displacement;
    return VON_KARMAN_CONSTANT * ref_wind /
           std::log(effective_height / roughness_length);
}
 
// ============================================================================
// Common Roughness Lengths
// ============================================================================
 
namespace roughness {
inline constexpr double OPEN_WATER = 0.0002; // [m]
inline constexpr double OPEN_TERRAIN = 0.03; // [m] - grass, few trees
inline constexpr double RURAL = 0.1;         // [m] - scattered buildings
inline constexpr double SUBURBAN = 0.5;      // [m] - low-rise urban
inline constexpr double URBAN = 1.0;         // [m] - city center
} // namespace roughness
 
// ============================================================================
// Common Power Law Exponents
// ============================================================================
 
namespace exponent {
inline constexpr double UNSTABLE =
    1.0 / 7.0;                               // ~0.143 - daytime, strong heating
inline constexpr double NEUTRAL = 1.0 / 4.0; // 0.25 - overcast, moderate wind
inline constexpr double STABLE = 1.0 / 3.0;  // ~0.333 - nighttime, light wind
} // namespace exponent
 
} // namespace vulcan::wind_shear