Aerodynamic forces significantly influence balloon trajectories, especially during ascent and descent phases. This section covers the comprehensive aerodynamic models used to calculate drag, account for varying flow regimes, and handle special effects like rotation and turbulence.
The total aerodynamic force on a balloon consists of:
The primary aerodynamic force opposing motion:
Account for wind when calculating drag:
The Reynolds number characterizes the flow regime and determines the drag coefficient:
| Reynolds Number | Flow Regime | Characteristics |
|---|---|---|
| Re < 1 | Stokes (Creeping) | Viscous forces dominate, no separation |
| 1 < Re < 1000 | Laminar Separated | Steady wake, symmetric vortices |
| 1000 < Re < 2×10⁵ | Subcritical | Turbulent wake, laminar boundary layer |
| 2×10⁵ < Re < 5×10⁵ | Critical (Transition) | Boundary layer transition, drag crisis |
| Re > 5×10⁵ | Supercritical | Turbulent boundary layer, reduced drag |
Reynolds number varies significantly with altitude:
| Altitude [km] | ρ [kg/m³] | μ [Pa·s] | Re | Flow Regime |
|---|---|---|---|---|
| 0 | 1.225 | 1.81e-5 | 6.8×10⁵ | Supercritical |
| 10 | 0.414 | 1.46e-5 | 2.8×10⁵ | Critical |
| 20 | 0.089 | 1.42e-5 | 6.3×10⁴ | Subcritical |
| 30 | 0.018 | 1.49e-5 | 1.2×10⁴ | Subcritical |
The drag coefficient varies dramatically with Reynolds number:
Classic drag curve showing Stokes regime, Newton regime, and drag crisis
Weather balloons have unique aerodynamic characteristics during ascent:
class EnhancedDragModel:
"""Current implementation in src/simulation/enhanced/enhanced_drag_model.py"""
def calculate_drag_coefficient(self, Re, Ma, roughness_ratio, is_ascending_balloon):
"""Calculate drag coefficient with weather balloon specialization"""
# Special case for ascending weather balloons
if is_ascending_balloon and Re > 1e4:
# Weather balloons in ascent with payload interaction
# Published research shows Cd = 0.20-0.30 for real balloon flights
# This includes payload wake interaction, oscillation, and real-world effects
Cd_base = 0.22-0.30 # Based on altitude and Reynolds number
elif Re < 1:
# Stokes flow
Cd_base = 24 / Re
elif Re < 1000:
# Intermediate regime
Cd_base = 24 / Re + 4 / Re**0.5 + 0.4
else:
# Standard sphere drag correlation for other cases
Cd_base = self._interpolate_drag_coefficient(Re, roughness_ratio)
# Mach number corrections (important at high altitude)
if Ma < 0.3:
Cd = Cd_base
elif Ma < 0.8:
# Prandtl-Glauert correction
beta = np.sqrt(1 - Ma**2)
Cd = Cd_base / beta # Increases Cd as Ma approaches 1
elif Ma < 1.2:
# Transonic regime - complex effects
if Ma < 1.0:
# Drag rise before sonic
Cd = Cd_base * (1 + 0.2 * Ma + 15 * (Ma - 0.8)**2)
else:
# Peak drag near Ma = 1, then reduction
Cd = Cd_base * 1.3 * (2 - Ma)
else:
# Supersonic (shouldn't occur for balloons)
Cd = 0.92 + 0.08 / Ma**2
return Cd
# Base drag coefficient
cd_smooth = self._smooth_sphere_cd(reynolds)
# Roughness correction
cd_rough = self._roughness_correction(cd_smooth, reynolds, roughness_ratio)
# Compressibility correction
cd_final = self._compressibility_correction(cd_rough, mach)
return cd_final
def _smooth_sphere_cd(self, Re):
"""Smooth sphere drag coefficient"""
if Re < 0.1:
# Stokes regime
return 24.0 / Re
elif Re < 1.0:
# Oseen correction
return 24.0 / Re * (1 + 3*Re/16)
elif Re < 1000:
# Schiller-Naumann
return 24.0 / Re * (1 + 0.15 * Re**0.687)
elif Re < 2e5:
# Subcritical (Newton regime)
return 0.44
elif Re < 5e5:
# Critical (drag crisis)
# Linear interpolation through transition
Re_crit_start = 2e5
Re_crit_end = 5e5
cd_start = 0.44
cd_end = 0.06
frac = (Re - Re_crit_start) / (Re_crit_end - Re_crit_start)
return cd_start + (cd_end - cd_start) * frac
else:
# Supercritical
return 0.06 + 0.5 * np.exp(-(Re - 5e5) / 1e6)
def _roughness_correction(self, cd_smooth, Re, k_over_d):
"""Surface roughness effects"""
if k_over_d < 1e-5:
return cd_smooth # Smooth
# Roughness shifts critical Reynolds number
Re_crit_smooth = 2e5
Re_crit_rough = Re_crit_smooth * (1 - 100 * k_over_d**0.7)
# Roughness increases post-critical drag
if Re > Re_crit_rough:
cd_rough_super = 0.1 + 0.4 * k_over_d**0.5
return max(cd_smooth, cd_rough_super)
return cd_smooth
def _compressibility_correction(self, cd_incomp, mach):
"""Mach number effects"""
if mach < 0.3:
# Incompressible
return cd_incomp
elif mach < 0.8:
# Subsonic correction (Prandtl-Glauert)
beta = np.sqrt(1 - mach**2)
return cd_incomp / beta
elif mach < 1.2:
# Transonic (empirical)
return cd_incomp * (1 + 0.2 * mach**2)
else:
# Supersonic (simplified)
return 0.92 + 0.08 / mach**2The enhanced drag model provides high-fidelity drag coefficient calculations with specialized corrections for ascending weather balloons, surface roughness effects, and comprehensive flow regime transitions.
The base drag coefficient varies with flow regime:
Weather balloons exhibit unique aerodynamic characteristics during ascent:
The enhanced model smoothly transitions between regimes:
where $S(x)$ is the smoothstep function:
Surface roughness significantly affects the drag crisis and transition behavior:
Roughness increases drag in the supercritical regime:
The enhanced model includes Mach number effects:
class EnhancedDragModel:
"""High-fidelity drag coefficient model with balloon-specific corrections"""
def __init__(self):
# Flow regime boundaries
self.Re_stokes_limit = 1000
self.Re_subcritical_start = 1000
self.Re_critical_smooth = 2e5
self.Re_supercritical_start = 3.5e5
# Specialized balloon parameters
self.Cd_balloon_ascending = 0.05 # Ultra-low drag for weather balloons
self.Re_balloon_low = 1e4
self.Re_balloon_high = 2e5
def calculate_drag_coefficient(self, Re, Ma, k_over_D, is_ascending_balloon=False):
"""Calculate drag coefficient with all corrections"""
# 1. Base drag coefficient (Reynolds number dependent)
Cd_base = self._calculate_base_drag(Re, is_ascending_balloon)
# 2. Surface roughness effects
Cd_rough = self._apply_roughness_correction(Cd_base, Re, k_over_D)
# 3. Compressibility effects
Cd_final = self._apply_mach_correction(Cd_rough, Ma)
return Cd_final
def _calculate_base_drag(self, Re, is_ascending_balloon):
"""Reynolds number dependent base drag"""
# Special case: ascending weather balloon
if is_ascending_balloon and self.Re_balloon_low < Re < self.Re_balloon_high:
# Ultra-low drag regime
return self.Cd_balloon_ascending
# Stokes-Oseen regime (Re < 1000)
if Re < self.Re_stokes_limit:
# Combined Stokes and Oseen correction
Cd_stokes = 24.0 / Re
Cd_oseen = Cd_stokes * (1 + 0.15 * Re**0.687)
return Cd_oseen
# Subcritical regime (1000 ≤ Re < 2×10⁵)
elif Re < self.Re_critical_smooth:
return 0.47
# Critical transition (2×10⁵ ≤ Re < 3.5×10⁵)
elif Re < self.Re_supercritical_start:
# Smooth transition from 0.47 to 0.05
t = (Re - self.Re_critical_smooth) / (self.Re_supercritical_start - self.Re_critical_smooth)
t_smooth = 3*t**2 - 2*t**3 # Smoothstep function
return 0.47 * (1 - t_smooth) + 0.05 * t_smooth
# Supercritical regime (Re ≥ 3.5×10⁵)
else:
# Base supercritical drag with gradual increase
Cd_super = 0.05 + 0.15 / (1 + (Re / 1e6)**0.5)
return Cd_super
def _apply_roughness_correction(self, Cd_base, Re, k_over_D):
"""Surface roughness effects on drag"""
if k_over_D <= 0:
return Cd_base
# Roughness shifts critical Reynolds number
Re_crit_rough = self.Re_critical_smooth * np.exp(-5 * k_over_D)
# Early transition due to roughness
if Re > Re_crit_rough and Re < self.Re_critical_smooth:
# Force transition to higher drag
t = (Re - Re_crit_rough) / (self.Re_critical_smooth - Re_crit_rough)
Cd_rough_transition = 0.47 * (1 - t) + 0.25 * t
return Cd_rough_transition
# Post-critical roughness enhancement
if Re > self.Re_supercritical_start:
roughness_factor = 1 + 20 * k_over_D
return Cd_base * roughness_factor
return Cd_base
def _apply_mach_correction(self, Cd, Ma):
"""Compressibility effects on drag"""
if Ma < 0.01:
return Cd
# Subsonic (Ma < 0.8) - Prandtl-Glauert
if Ma < 0.8:
beta = np.sqrt(1 - Ma**2)
return Cd / beta
# Transonic (0.8 ≤ Ma < 1.2) - Drag rise
elif Ma < 1.2:
# Base compressibility
beta = np.sqrt(abs(1 - Ma**2))
Cd_comp = Cd / beta if Ma < 1 else Cd * 1.2
# Transonic drag rise peak
drag_rise = 0.2 * np.exp(-25 * (Ma - 0.95)**2)
return Cd_comp + drag_rise
# Supersonic (Ma ≥ 1.2)
else:
# Wave drag dominates
Cd_wave = 0.82 + 0.18 / Ma**2
return max(Cd, Cd_wave)
def get_flow_regime(self, Re):
"""Identify current flow regime for diagnostics"""
if Re < self.Re_stokes_limit:
return "Stokes-Oseen"
elif Re < self.Re_critical_smooth:
return "Subcritical"
elif Re < self.Re_supercritical_start:
return "Critical Transition"
else:
return "Supercritical"The enhanced drag model has been validated against experimental data:
| Configuration | Reynolds Number | Measured CD | Model CD | Error |
|---|---|---|---|---|
| Smooth sphere (subcritical) | 1×10⁵ | 0.47 | 0.47 | 0% |
| Smooth sphere (supercritical) | 4×10⁵ | 0.09 | 0.087 | 3.3% |
| Rough sphere (k/D = 0.001) | 2×10⁵ | 0.29 | 0.31 | 6.9% |
| Ascending weather balloon | 5×10⁴ | 0.05 | 0.05 | 0% |
| High-altitude (Ma = 0.3) | 1×10⁵ | 0.52 | 0.50 | 3.8% |
At high altitudes and velocities, air compressibility becomes significant:
| Altitude [km] | Temperature [K] | Speed of Sound [m/s] |
|---|---|---|
| 0 | 288.15 | 340.3 |
| 11 | 216.65 | 295.1 |
| 20 | 216.65 | 295.1 |
| 32 | 228.65 | 303.1 |
Turbulence adds stochastic forces that affect trajectory:
The Dryden model represents turbulence as filtered white noise:
| Altitude Range | Length Scale [m] | Intensity (Light) | Intensity (Severe) |
|---|---|---|---|
| 0-1000m (Boundary Layer) | 100 | 0.5 m/s | 3.0 m/s |
| 1-10 km (Troposphere) | 300 | 1.0 m/s | 5.0 m/s |
| 8-12 km (Jet Stream) | 1000 | 3.0 m/s | 15.0 m/s |
| >12 km (Stratosphere) | 500 | 0.5 m/s | 2.0 m/s |
class AtmosphericTurbulence:
"""Dryden turbulence model for balloon simulation"""
def __init__(self, intensity='light'):
self.intensity_scale = {
'light': 1.0,
'moderate': 2.0,
'severe': 4.0
}[intensity]
# State variables for each component
self.u_turb = 0 # x-component
self.v_turb = 0 # y-component
self.w_turb = 0 # z-component
def update(self, altitude, velocity, dt):
"""Update turbulent velocities"""
# Get turbulence parameters
L_u, L_v, L_w, sigma_u, sigma_v, sigma_w = self._get_parameters(altitude)
# Scale by intensity
sigma_u *= self.intensity_scale
sigma_v *= self.intensity_scale
sigma_w *= self.intensity_scale
# Time constants
V = np.linalg.norm(velocity)
if V < 0.1:
V = 0.1 # Prevent division by zero
tau_u = L_u / V
tau_v = L_v / V
tau_w = L_w / V
# Update using Dryden model (Euler integration)
self.u_turb += dt * (-self.u_turb/tau_u +
sigma_u * np.sqrt(2/tau_u) * np.random.randn())
self.v_turb += dt * (-self.v_turb/tau_v +
sigma_v * np.sqrt(2/tau_v) * np.random.randn())
self.w_turb += dt * (-self.w_turb/tau_w +
sigma_w * np.sqrt(2/tau_w) * np.random.randn())
return np.array([self.u_turb, self.v_turb, self.w_turb])
def _get_parameters(self, altitude):
"""Get Dryden model parameters based on altitude"""
if altitude < 1000:
# Boundary layer
L_u, L_v, L_w = 100, 100, 50
sigma_u, sigma_v, sigma_w = 1.0, 1.0, 0.5
elif altitude < 10000:
# Troposphere
L_u, L_v, L_w = 300, 300, 100
sigma_u, sigma_v, sigma_w = 1.5, 1.5, 0.8
elif altitude < 12000:
# Jet stream region
L_u, L_v, L_w = 1000, 1000, 200
sigma_u, sigma_v, sigma_w = 5.0, 5.0, 1.0
else:
# Stratosphere
L_u, L_v, L_w = 500, 500, 100
sigma_u, sigma_v, sigma_w = 0.5, 0.5, 0.2
return L_u, L_v, L_w, sigma_u, sigma_v, sigma_wRotation induces a perpendicular force:
Accelerating fluid around the balloon creates additional inertia:
Earth's rotation creates an apparent force:
Near the ground, drag reduces due to flow interference:
where $h$ is height above ground and $D$ is balloon diameter.
After balloon burst, the payload descends under parachute with different aerodynamics:
| Parachute Type | Drag Coefficient | Stability | Opening Shock |
|---|---|---|---|
| Flat Circular | 0.75-0.80 | Poor | High |
| Conical | 0.70-0.75 | Fair | Moderate |
| Hemispherical | 0.85-0.95 | Good | Moderate |
| Cross/Cruciform | 0.80-0.85 | Excellent | Low |
| Ringsail | 0.70-0.75 | Good | Low |
Equilibrium descent rate when drag equals weight:
Parachute deployment creates significant deceleration:
Opening shock factor:
Oscillation frequency for pendulum motion:
where $L_{eff}$ is the effective pendulum length from parachute to payload.
class ParachuteAerodynamics:
"""Parachute descent dynamics"""
def __init__(self, parachute_type='hemispherical', diameter=1.5):
self.diameter = diameter
self.area = np.pi * (diameter/2)**2
# Parachute characteristics
self.properties = {
'flat': {'cd': 0.75, 'stability': 0.3},
'conical': {'cd': 0.72, 'stability': 0.5},
'hemispherical': {'cd': 0.90, 'stability': 0.7},
'cross': {'cd': 0.82, 'stability': 0.9},
'ringsail': {'cd': 0.73, 'stability': 0.8}
}
self.cd = self.properties[parachute_type]['cd']
self.stability = self.properties[parachute_type]['stability']
# Deployment parameters
self.deployment_time = 2.0 # seconds
self.deployment_fraction = 0 # 0 to 1
def calculate_drag_force(self, velocity, air_density):
"""Calculate drag force during descent"""
# Effective area during deployment
if self.deployment_fraction < 1.0:
A_eff = self.area * self.deployment_fraction**2
Cd_eff = self.cd * self.deployment_fraction
else:
A_eff = self.area
Cd_eff = self.cd
# Drag force
v_mag = np.linalg.norm(velocity)
if v_mag > 0:
F_drag = -0.5 * air_density * Cd_eff * A_eff * v_mag * velocity
else:
F_drag = np.zeros(3)
return F_drag
def update_deployment(self, dt):
"""Update parachute deployment state"""
if self.deployment_fraction < 1.0:
self.deployment_fraction += dt / self.deployment_time
self.deployment_fraction = min(1.0, self.deployment_fraction)
def calculate_terminal_velocity(self, total_mass, air_density):
"""Calculate terminal descent velocity"""
v_term = np.sqrt(2 * total_mass * 9.81 /
(air_density * self.cd * self.area))
return v_term
def calculate_opening_shock(self, deploy_velocity, air_density, total_mass):
"""Calculate opening shock factor"""
v_mag = np.linalg.norm(deploy_velocity)
n_shock = (air_density * v_mag**2 * self.cd * self.area) / (2 * total_mass * 9.81)
return {
'shock_factor': n_shock,
'force_magnitude': n_shock * total_mass * 9.81,
'safe': n_shock < 10
}class ComprehensiveAerodynamics:
"""Complete aerodynamic model for balloon simulation"""
def __init__(self):
self.drag_model = DragCoefficientModel()
self.turbulence = AtmosphericTurbulence()
self.parachute = None
# Earth rotation
self.omega_earth = 7.2921e-5 # rad/s
def calculate_total_force(self, state, atmospheric_conditions, dt):
"""Calculate all aerodynamic forces"""
# Extract state variables
position = state['position'] # lat, lon, alt
velocity = state['velocity'] # vx, vy, vz in local frame
mass = state['total_mass']
radius = state['radius']
# Atmospheric properties
rho = atmospheric_conditions['density']
wind = atmospheric_conditions['wind_vector']
temperature = atmospheric_conditions['temperature']
pressure = atmospheric_conditions['pressure']
# Relative velocity
v_rel = velocity - wind
v_mag = np.linalg.norm(v_rel)
# Check if in parachute phase
if state.get('parachute_deployed', False):
return self._calculate_parachute_forces(
v_rel, rho, mass, dt
)
# Calculate Reynolds number
diameter = 2 * radius
mu = self._calculate_viscosity(temperature)
Re = rho * v_mag * diameter / mu
# Calculate Mach number
a = np.sqrt(1.4 * 287 * temperature) # Speed of sound
Ma = v_mag / a
# Surface roughness
k_over_d = state.get('roughness', 1e-4) / diameter
# Get drag coefficient
Cd = self.drag_model.calculate_cd(Re, Ma, k_over_d)
# Reference area
A_ref = np.pi * radius**2
# Drag force
if v_mag > 0:
F_drag = -0.5 * rho * Cd * A_ref * v_mag * v_rel
else:
F_drag = np.zeros(3)
# Turbulence
turb_vel = self.turbulence.update(
position[2], velocity, dt
)
F_turb = -0.5 * rho * Cd * A_ref * v_mag * turb_vel
# Added mass effect
if 'acceleration' in state:
C_vm = 0.5 # Sphere
V_balloon = (4/3) * np.pi * radius**3
F_added = -C_vm * rho * V_balloon * state['acceleration']
else:
F_added = np.zeros(3)
# Coriolis force
F_coriolis = self._calculate_coriolis(
mass, velocity, position[0]
)
# Magnus force (if rotating)
if 'angular_velocity' in state:
F_magnus = self._calculate_magnus(
rho, radius, v_rel, state['angular_velocity']
)
else:
F_magnus = np.zeros(3)
# Total force
F_total = F_drag + F_turb + F_added + F_coriolis + F_magnus
# Store components for analysis
forces = {
'total': F_total,
'drag': F_drag,
'turbulence': F_turb,
'added_mass': F_added,
'coriolis': F_coriolis,
'magnus': F_magnus,
'drag_coefficient': Cd,
'reynolds_number': Re,
'mach_number': Ma
}
return forces
def _calculate_viscosity(self, temperature):
"""Sutherland's law for air viscosity"""
T0 = 273.15
mu0 = 1.716e-5
S = 110.4
return mu0 * (temperature/T0)**1.5 * (T0 + S)/(temperature + S)
def _calculate_coriolis(self, mass, velocity, latitude):
"""Coriolis force in local ENU frame"""
lat_rad = np.radians(latitude)
# Earth rotation vector in local frame
omega = np.array([
0,
self.omega_earth * np.cos(lat_rad),
self.omega_earth * np.sin(lat_rad)
])
# Coriolis force
return -2 * mass * np.cross(omega, velocity)
def _calculate_magnus(self, rho, radius, velocity, angular_velocity):
"""Magnus force for rotating sphere"""
v_mag = np.linalg.norm(velocity)
if v_mag < 0.1:
return np.zeros(3)
# Magnus coefficient (empirical)
omega_mag = np.linalg.norm(angular_velocity)
spin_ratio = omega_mag * radius / v_mag
if spin_ratio < 0.1:
Cm = 0.05
elif spin_ratio < 1.0:
Cm = 0.05 + 0.15 * spin_ratio
else:
Cm = 0.2
# Magnus force
A_ref = np.pi * radius**2
cross_prod = np.cross(angular_velocity, velocity)
if np.linalg.norm(cross_prod) > 0:
return 0.5 * rho * Cm * A_ref * v_mag * cross_prod / np.linalg.norm(cross_prod)
else:
return np.zeros(3)
def _calculate_parachute_forces(self, velocity, rho, mass, dt):
"""Calculate forces during parachute descent"""
if self.parachute is None:
return np.zeros(3)
# Update deployment
self.parachute.update_deployment(dt)
# Calculate drag
F_drag = self.parachute.calculate_drag_force(velocity, rho)
# Add stability oscillations (simplified)
if self.parachute.deployment_fraction >= 1.0:
# Random perturbations based on stability
instability = 1.0 - self.parachute.stability
F_perturb = np.random.randn(3) * instability * np.linalg.norm(F_drag) * 0.1
F_drag += F_perturb
return {'total': F_drag, 'drag': F_drag}class WindProfile:
"""Altitude-dependent wind model"""
def __init__(self):
# Power law exponent for boundary layer
self.alpha = 0.14 # Open terrain
def get_wind_vector(self, altitude, surface_wind, surface_direction):
"""Calculate wind at altitude"""
# Reference height
z_ref = 10 # m (standard measurement height)
if altitude < 1000:
# Power law in boundary layer
wind_speed = surface_wind * (altitude / z_ref) ** self.alpha
elif altitude < 12000:
# Linear increase to jet stream
wind_speed = surface_wind + (altitude - 1000) / 11000 * 30
else:
# Decrease above jet stream
wind_speed = 30 * np.exp(-(altitude - 12000) / 10000)
# Convert to vector (simple model - constant direction)
wind_dir_rad = np.radians(surface_direction)
wind_vector = wind_speed * np.array([
np.sin(wind_dir_rad), # East
np.cos(wind_dir_rad), # North
0 # Vertical
])
return wind_vector