From 71865692595f8a45dc21d038e54c0e0a2bbe85be Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludvig=20B=C3=B6klin?= Date: Sat, 24 Apr 2021 18:11:18 +0200 Subject: [PATCH] Reduce abstraction for lift calculation; remove RigidWings struct --- common/src/comp/fluid_dynamics.rs | 260 ++++++++++++------------------ common/src/comp/mod.rs | 2 +- common/src/states/glide.rs | 11 +- common/systems/src/phys.rs | 17 +- 4 files changed, 114 insertions(+), 176 deletions(-) diff --git a/common/src/comp/fluid_dynamics.rs b/common/src/comp/fluid_dynamics.rs index 23589ef285..6afd527fcf 100644 --- a/common/src/comp/fluid_dynamics.rs +++ b/common/src/comp/fluid_dynamics.rs @@ -1,6 +1,6 @@ use super::{ body::{object, Body}, - Density, Ori, Vel, + CharacterState, Density, Ori, Vel, }; use crate::{ consts::{AIR_DENSITY, WATER_DENSITY}, @@ -90,10 +90,9 @@ impl Default for Fluid { impl Body { pub fn aerodynamic_forces( &self, - ori: &Ori, rel_flow: &Vel, fluid_density: f32, - wings: Option<&RigidWings>, + character_state: Option<&CharacterState>, ) -> Vec3 { let v_sq = rel_flow.0.magnitude_squared(); if v_sq < 0.25 { @@ -104,39 +103,41 @@ impl Body { // All the coefficients come pre-multiplied by their reference area 0.5 * fluid_density * v_sq - * wings - .map(|wings| { - // Since we have wings, we proceed to calculate the lift and drag + * character_state + .and_then(|cs| match cs { + CharacterState::Glide(data) => { + Some((data.aspect_ratio, data.planform_area, data.ori)) + }, + _ => None, + }) + .map(|(ar, area, ori)| { + // We have an elliptical wing; proceed to calculate its lift and drag - let ar = wings.aspect_ratio(); // aoa will be positive when we're pitched up and negative otherwise - let aoa = angle_of_attack(ori, &rel_flow_dir); + let aoa = angle_of_attack(&ori, &rel_flow_dir); // c_l will be positive when aoa is positive (we have positive lift, // producing an upward force) and negative otherwise - let c_l = wings.lift_coefficient(aoa); + let c_l = lift_coefficient(ar, area, aoa); // lift dir will be orthogonal to the local relative flow vector. - // Local relative flow is the resulting vector of (relative) freestream flow - // + downwash (created by the vortices of the wing tips) + // Local relative flow is the resulting vector of (relative) freestream + // flow + downwash (created by the vortices + // of the wing tips) let lift_dir: Dir = { // induced angle of attack let aoa_i = c_l / (PI * ar); - // effective angle of attack; the aoa as seen by aerofoil after downwash + // effective angle of attack; the aoa as seen by aerofoil after + // downwash let aoa_eff = aoa - aoa_i; - /*println!( - "CL={:.1}, α={:.1}°, αᵢ={:.1}°, αₑ={:.1}°, AR={:.1}", - c_l, - aoa.to_degrees(), - aoa_i.to_degrees(), - aoa_eff.to_degrees(), - ar - );*/ - // Angle between chord line and local relative wind is aoa_eff radians. - // Direction of lift is perpendicular to local relative wind. - // At positive lift, local relative wind will be below our cord line at - // an angle of aoa_eff. Thus if we pitch down by aoa_eff radians then - // our chord line will be colinear with local relative wind vector and - // our up will be the direction of lift. + // Angle between chord line and local relative wind is aoa_eff + // radians. Direction of lift is + // perpendicular to local relative wind. + // At positive lift, local relative wind will be below our cord line + // at an angle of aoa_eff. Thus if + // we pitch down by aoa_eff radians then + // our chord line will be colinear with local relative wind vector + // and our up will be the direction + // of lift. ori.pitched_down(aoa_eff).up() }; @@ -146,20 +147,12 @@ impl Body { // (this definition should not be used for aspect ratios > 25) let e = 1.78 * (1.0 - 0.045 * ar.powf(0.68)) - 0.64; - wings.zero_lift_drag_coefficient() + zero_lift_drag_coefficient(area) + self.parasite_drag_coefficient() + c_l.powi(2) / (PI * e * ar) }; debug_assert!(c_d.is_sign_positive()); debug_assert!(c_l.is_sign_positive() || aoa.is_sign_negative()); - /*println!( - "L/D (at α={:.1}, AR={:.1}) = {:.1}/{:.1} = {:.1}", - aoa.to_degrees(), - ar, - 0.5 * fluid_density * v_sq * c_l, - 0.5 * fluid_density * v_sq * c_d, - c_l / c_d - );*/ c_l * *lift_dir + c_d * *rel_flow_dir }) @@ -192,24 +185,28 @@ impl Body { } else { 1.0 }; - cd * std::f32::consts::PI * dim.x * dim.z + cd * PI * dim.x * dim.z }, // Cross-section, zero-lift angle; exclude the wings (width * 0.2) Body::BirdMedium(_) | Body::BirdLarge(_) | Body::Dragon(_) => { let dim = self.dimensions().map(|a| a * 0.5); + // "Field Estimates of Body Drag Coefficient on the Basis of Dives in Passerine + // Birds", Anders Hedenström and Felix Liechti, 2001 let cd = match self { - Body::BirdMedium(_) => 0.2, - Body::BirdLarge(_) => 0.4, + Body::BirdLarge(_) | Body::BirdMedium(_) => 0.2, + // arbitrary _ => 0.7, }; - cd * std::f32::consts::PI * dim.x * 0.2 * dim.z + cd * PI * dim.x * 0.2 * dim.z }, // Cross-section, zero-lift angle; exclude the fins (width * 0.2) Body::FishMedium(_) | Body::FishSmall(_) => { let dim = self.dimensions().map(|a| a * 0.5); - 0.031 * std::f32::consts::PI * dim.x * 0.2 * dim.z + // "A Simple Method to Determine Drag Coefficients in Aquatic Animals", + // D. Bilo and W. Nachtigall, 1980 + 0.031 * PI * dim.x * 0.2 * dim.z }, Body::Object(object) => match object { @@ -225,7 +222,7 @@ impl Body { | object::Body::FireworkYellow | object::Body::MultiArrow => { let dim = self.dimensions().map(|a| a * 0.5); - 0.02 * std::f32::consts::PI * dim.x * dim.z + 0.02 * PI * dim.x * dim.z }, // spherical-ish objects @@ -243,12 +240,12 @@ impl Body { | object::Body::Pumpkin4 | object::Body::Pumpkin5 => { let dim = self.dimensions().map(|a| a * 0.5); - 0.5 * std::f32::consts::PI * dim.x * dim.z + 0.5 * PI * dim.x * dim.z }, _ => { let dim = self.dimensions(); - 2.0 * (std::f32::consts::PI / 6.0 * dim.x * dim.y * dim.z).powf(2.0 / 3.0) + 2.0 * (PI / 6.0 * dim.x * dim.y * dim.z).powf(2.0 / 3.0) }, }, @@ -256,7 +253,7 @@ impl Body { // Airships tend to use the square of the cube root of its volume for // reference area let dim = self.dimensions(); - (std::f32::consts::PI / 6.0 * dim.x * dim.y * dim.z).powf(2.0 / 3.0) + (PI / 6.0 * dim.x * dim.y * dim.z).powf(2.0 / 3.0) }, } } @@ -267,122 +264,73 @@ fn angle_of_attack(ori: &Ori, rel_flow_dir: &Dir) -> f32 { PI / 2.0 - ori.up().angle_between(rel_flow_dir.to_vec()) } -/// An elliptical fixed rigid wing. Plurally named simply because it's a shape -/// typically composed of two wings forming an elliptical lift distribution. -// -// Animal wings are technically flexible, not rigid, (difference being that the -// former's shape is affected by the flow) and usually has the ability to -// assume complex shapes with properties like curved camber line, span-wise -// twist, dihedral angle, sweep angle, and partitioned sections. However, we -// could make do with this model for fully extended animal wings, enabling them -// to glide. -#[derive(Copy, Clone, Debug, PartialEq, Serialize, Deserialize)] -pub struct RigidWings { - aspect_ratio: f32, - planform_area: f32, - // sweep_angle: Option, -} - -impl RigidWings { - /// Wings from total span (wing-tip to wing-tip) and - /// chord length (leading edge to trailing edge) - pub fn new(span_length: f32, chord_length: f32) -> Self { - let planform_area = std::f32::consts::PI * chord_length * span_length * 0.25; - Self { - aspect_ratio: span_length.powi(2) / planform_area, - planform_area, - } - } - - /// The aspect ratio is the ratio of the span squared to actual planform - /// area - pub fn aspect_ratio(&self) -> f32 { self.aspect_ratio } - - pub fn planform_area(&self) -> f32 { self.planform_area } -} - -impl RigidWings { - /// Total lift coefficient for a finite wing of symmetric aerofoil shape and - /// elliptical pressure distribution. - pub fn lift_coefficient(&self, aoa: f32) -> f32 { - let aoa_abs = aoa.abs(); - let stall_angle = PI * 0.1; - inline_tweak::tweak!(1.0) - * self.planform_area() - * if aoa_abs < stall_angle { - self.lift_slope(None) * aoa - } else if inline_tweak::tweak!(true) { - // This is when flow separation and turbulence starts to kick in. - // Going to just make something up (based on some data), as the alternative is - // to just throw your hands up and return 0 - let aoa_s = aoa.signum(); - let c_l_max = self.lift_slope(None) * stall_angle; - let deg_45 = PI / 4.0; - if aoa_abs < deg_45 { - // drop directly to 0.6 * max lift at stall angle - // then climb back to max at 45° - Lerp::lerp(0.6 * c_l_max, c_l_max, aoa_abs / deg_45) * aoa_s - } else { - // let's just say lift goes down linearly again until we're at 90° - Lerp::lerp(c_l_max, 0.0, (aoa_abs - deg_45) / deg_45) * aoa_s - } +/// Total lift coefficient for a finite wing of symmetric aerofoil shape and +/// elliptical pressure distribution. +pub fn lift_coefficient(aspect_ratio: f32, planform_area: f32, aoa: f32) -> f32 { + let aoa_abs = aoa.abs(); + let stall_angle = PI * 0.1; + inline_tweak::tweak!(1.0) + * planform_area + * if aoa_abs < stall_angle { + lift_slope(aspect_ratio, None) * aoa + } else if inline_tweak::tweak!(true) { + // This is when flow separation and turbulence starts to kick in. + // Going to just make something up (based on some data), as the alternative is + // to just throw your hands up and return 0 + let aoa_s = aoa.signum(); + let c_l_max = lift_slope(aspect_ratio, None) * stall_angle; + let deg_45 = PI / 4.0; + if aoa_abs < deg_45 { + // drop directly to 0.6 * max lift at stall angle + // then climb back to max at 45° + Lerp::lerp(0.6 * c_l_max, c_l_max, aoa_abs / deg_45) * aoa_s } else { - 0.0 + // let's just say lift goes down linearly again until we're at 90° + Lerp::lerp(c_l_max, 0.0, (aoa_abs - deg_45) / deg_45) * aoa_s } - } - - /// The zero-lift profile drag coefficient is the parasite drag on the wings - /// at the angle of attack which generates no lift - pub fn zero_lift_drag_coefficient(&self) -> f32 { - // avg value for Harris' hawk (Parabuteo unicinctus) [1] - self.planform_area() * 0.02 - } - - /// The change in lift over change in angle of attack¹. Multiplying by angle - /// of attack gives the lift coefficient (for a finite wing, not aerofoil). - /// - /// Aspect ratio is the ratio of total wing span squared over planform area. - /// - /// # Notes - /// - /// Only valid for symmetric, elliptical wings at small² angles of attack³. - /// Does not apply to twisted, cambered or delta wings. (It still gives a - /// reasonably accurate approximation if the wing shape is not truly - /// elliptical.) - /// - /// 1. geometric angle of attack, i.e. the pitch angle relative to - /// freestream flow - /// 2. up to around ~18°, at which point maximum lift has been achieved and - /// thereafter falls precipitously, causing a stall (this is the stall - /// angle) 3. effective aoa, i.e. geometric aoa - induced aoa; assumes - /// no sideslip - fn lift_slope(&self, sweep_angle: Option) -> f32 { - // lift slope for a thin aerofoil, given by Thin Aerofoil Theory - let ar = self.aspect_ratio(); - let a0 = 2.0 * PI; - if let Some(sweep) = sweep_angle { - // for swept wings we use Kuchemann's modification to Helmbold's - // equation - let a0_cos_sweep = a0 * sweep.cos(); - let x = a0_cos_sweep / (PI * ar); - a0_cos_sweep / ((1.0 + x.powi(2)).sqrt() + x) - } else if ar < 4.0 { - // for low aspect ratio wings (AR < 4) we use Helmbold's equation - let x = a0 / (PI * ar); - a0 / ((1.0 + x.powi(2)).sqrt() + x) } else { - // for high aspect ratio wings (AR > 4) we use the equation given by - // Prandtl's lifting-line theory - a0 / (1.0 + (a0 / (PI * ar))) + 0.0 } - } } -/* -## References: +/// The zero-lift profile drag coefficient is the parasite drag on the wings +/// at the angle of attack which generates no lift +pub fn zero_lift_drag_coefficient(planform_area: f32) -> f32 { + // avg value for Harris' hawk (Parabuteo unicinctus) [1] + planform_area * 0.02 +} -1. "Field Estimates of Body Drag Coefficient on the Basis of Dives in Passerine Birds", - Anders Hedenström and Felix Liechti, 2001 -2. "A Simple Method to Determine Drag Coefficients in Aquatic Animals", - D. Bilo and W. Nachtigall, 1980 -*/ +/// The change in lift over change in angle of attack¹. Multiplying by angle +/// of attack gives the lift coefficient (for a finite wing, not aerofoil). +/// Aspect ratio is the ratio of total wing span squared over planform area. +/// +/// # Notes +/// Only valid for symmetric, elliptical wings at small² angles of attack³. +/// Does not apply to twisted, cambered or delta wings. (It still gives a +/// reasonably accurate approximation if the wing shape is not truly +/// elliptical.) +/// 1. geometric angle of attack, i.e. the pitch angle relative to +/// freestream flow +/// 2. up to around ~18°, at which point maximum lift has been achieved and +/// thereafter falls precipitously, causing a stall (this is the stall +/// angle) 3. effective aoa, i.e. geometric aoa - induced aoa; assumes +/// no sideslip +fn lift_slope(aspect_ratio: f32, sweep_angle: Option) -> f32 { + // lift slope for a thin aerofoil, given by Thin Aerofoil Theory + let a0 = 2.0 * PI; + if let Some(sweep) = sweep_angle { + // for swept wings we use Kuchemann's modification to Helmbold's + // equation + let a0_cos_sweep = a0 * sweep.cos(); + let x = a0_cos_sweep / (PI * aspect_ratio); + a0_cos_sweep / ((1.0 + x.powi(2)).sqrt() + x) + } else if aspect_ratio < 4.0 { + // for low aspect ratio wings (AR < 4) we use Helmbold's equation + let x = a0 / (PI * aspect_ratio); + a0 / ((1.0 + x.powi(2)).sqrt() + x) + } else { + // for high aspect ratio wings (AR > 4) we use the equation given by + // Prandtl's lifting-line theory + a0 / (1.0 + (a0 / (PI * aspect_ratio))) + } +} diff --git a/common/src/comp/mod.rs b/common/src/comp/mod.rs index 53e86e65fb..5f563a3836 100644 --- a/common/src/comp/mod.rs +++ b/common/src/comp/mod.rs @@ -68,7 +68,7 @@ pub use self::{ InputKind, InventoryAction, InventoryEvent, InventoryManip, MountState, Mounting, }, energy::{Energy, EnergyChange, EnergySource}, - fluid_dynamics::{Fluid, RigidWings}, + fluid_dynamics::Fluid, group::Group, home_chunk::HomeChunk, inputs::CanBuild, diff --git a/common/src/states/glide.rs b/common/src/states/glide.rs index 10f030a623..48856dd6e3 100644 --- a/common/src/states/glide.rs +++ b/common/src/states/glide.rs @@ -1,6 +1,6 @@ use super::utils::handle_climb; use crate::{ - comp::{inventory::slot::EquipSlot, CharacterState, Ori, RigidWings, StateUpdate}, + comp::{inventory::slot::EquipSlot, CharacterState, Ori, StateUpdate}, states::behavior::{CharacterBehavior, JoinData}, util::Dir, }; @@ -9,14 +9,19 @@ use vek::*; #[derive(Copy, Clone, Debug, PartialEq, Serialize, Deserialize)] pub struct Data { - pub wings: RigidWings, + /// The aspect ratio is the ratio of the span squared to actual planform + /// area + pub aspect_ratio: f32, + pub planform_area: f32, pub ori: Ori, } impl Data { pub fn new(span_length: f32, chord_length: f32, ori: Ori) -> Self { + let planform_area = std::f32::consts::PI * chord_length * span_length * 0.25; Self { - wings: RigidWings::new(span_length, chord_length), + aspect_ratio: span_length.powi(2) / planform_area, + planform_area, ori, } } diff --git a/common/systems/src/phys.rs b/common/systems/src/phys.rs index cbd3f19c85..24f387218e 100644 --- a/common/systems/src/phys.rs +++ b/common/systems/src/phys.rs @@ -9,7 +9,6 @@ use common::{ event::{EventBus, ServerEvent}, outcome::Outcome, resources::DeltaTime, - states, terrain::{Block, TerrainGrid}, uid::Uid, util::{Projection, SpatialGrid}, @@ -43,7 +42,6 @@ fn fluid_density(height: f32, fluid: &Fluid) -> Density { fn integrate_forces( dt: &DeltaTime, mut vel: Vel, - ori: &Ori, body: &Body, density: &Density, mass: &Mass, @@ -61,17 +59,7 @@ fn integrate_forces( // Aerodynamic/hydrodynamic forces if !rel_flow.0.is_approx_zero() { debug_assert!(!rel_flow.0.map(|a| a.is_nan()).reduce_or()); - let glider: Option<&states::glide::Data> = character_state.and_then(|cs| match cs { - CharacterState::Glide(data) => Some(data), - _ => None, - }); - let impulse = dt.0 - * body.aerodynamic_forces( - glider.map(|g| &g.ori).unwrap_or(ori), - &rel_flow, - fluid_density.0, - glider.map(|g| g.wings).as_ref(), - ); + let impulse = dt.0 * body.aerodynamic_forces(&rel_flow, fluid_density.0, character_state); debug_assert!(!impulse.map(|a| a.is_nan()).reduce_or()); if !impulse.is_approx_zero() { let new_v = vel.0 + impulse / mass.0; @@ -576,7 +564,6 @@ impl<'a> PhysicsData<'a> { ( positions, velocities, - &write.orientations, read.stickies.maybe(), &read.bodies, read.character_states.maybe(), @@ -595,7 +582,6 @@ impl<'a> PhysicsData<'a> { ( pos, vel, - ori, sticky, body, character_state, @@ -627,7 +613,6 @@ impl<'a> PhysicsData<'a> { vel.0 = integrate_forces( &dt, *vel, - ori, body, density, mass,