From 74f9945ab3fe00a601402b202641a186c9de680d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludvig=20B=C3=B6klin?= Date: Fri, 21 May 2021 11:13:47 +0200 Subject: [PATCH] Refactor and abstract fluid dynamics --- common/src/comp/body.rs | 126 ++++++++- common/src/comp/body/bird_large.rs | 73 ++++- common/src/comp/body/bird_medium.rs | 73 ++++- common/src/comp/body/dragon.rs | 73 ++++- common/src/comp/fluid_dynamics.rs | 414 ++++++++++++---------------- common/src/consts.rs | 2 +- common/src/lib.rs | 1 + common/src/states/glide.rs | 23 +- common/systems/src/phys.rs | 48 ++-- 9 files changed, 560 insertions(+), 273 deletions(-) diff --git a/common/src/comp/body.rs b/common/src/comp/body.rs index 5d51fc3b36..eb21faecdf 100644 --- a/common/src/comp/body.rs +++ b/common/src/comp/body.rs @@ -16,6 +16,7 @@ pub mod theropod; use crate::{ assets::{self, Asset}, + comp::{fluid_dynamics::*, Ori, Vel}, consts::{HUMAN_DENSITY, WATER_DENSITY}, make_case_elim, npc::NpcKind, @@ -23,6 +24,7 @@ use crate::{ use serde::{Deserialize, Serialize}; use specs::{Component, DerefFlaggedStorage}; use specs_idvs::IdvStorage; +use std::f32::consts::PI; use vek::*; use super::{BuffKind, Density, Mass}; @@ -316,9 +318,9 @@ impl Body { _ => Vec3::new(1.0, 0.75, 1.4), }, - Body::BirdMedium(_) => Vec3::new(2.0, 1.0, 1.5), - Body::BirdLarge(_) => Vec3::new(2.0, 6.0, 3.5), - Body::Dragon(_) => Vec3::new(16.0, 10.0, 16.0), + Body::BirdMedium(body) => body.dimensions(), + Body::BirdLarge(body) => body.dimensions(), + Body::Dragon(body) => body.dimensions(), Body::FishMedium(_) => Vec3::new(0.5, 2.0, 0.8), Body::FishSmall(_) => Vec3::new(0.3, 1.2, 0.6), Body::Golem(_) => Vec3::new(5.0, 5.0, 7.5), @@ -732,3 +734,121 @@ impl Body { impl Component for Body { type Storage = DerefFlaggedStorage>; } + +impl Body { + pub fn aerodynamic_forces( + &self, + ori: Option<&Ori>, + rel_flow: &Vel, + fluid_density: f32, + ) -> Vec3 { + match (self, ori) { + (Body::BirdMedium(bird), Some(ori)) => bird_medium::FlyingBirdMedium::from((bird, ori)) + .aerodynamic_forces(rel_flow, fluid_density), + (Body::BirdLarge(bird), Some(ori)) => bird_large::FlyingBirdLarge::from((bird, ori)) + .aerodynamic_forces(rel_flow, fluid_density), + (Body::Dragon(not_bird), Some(ori)) => dragon::FlyingDragon::from((not_bird, ori)) + .aerodynamic_forces(rel_flow, fluid_density), + _ => self.drag(rel_flow, fluid_density), + } + } +} + +impl Drag for Body { + /// Parasite drag is the sum of pressure drag and skin friction. + /// Skin friction is the drag arising from the shear forces between a fluid + /// and a surface, while pressure drag is due to flow separation. Both are + /// viscous effects. + /// + /// This coefficient includes the reference area. + fn parasite_drag_coefficient(&self) -> f32 { + // Reference area and drag coefficient assumes best-case scenario of the + // orientation producing least amount of drag + match self { + // Cross-section, head/feet first + Body::BipedLarge(_) | Body::BipedSmall(_) | Body::Golem(_) | Body::Humanoid(_) => { + const SCALE: f32 = 0.7; + let radius = self.dimensions().xy().map(|a| SCALE * a * 0.5); + const CD: f32 = 0.7; + CD * PI * radius.x * radius.y + }, + + // Cross-section, nose/tail first + Body::Theropod(_) + | Body::QuadrupedMedium(_) + | Body::QuadrupedSmall(_) + | Body::QuadrupedLow(_) => { + let radius = self.dimensions().map(|a| a * 0.5); + let cd: f32 = if matches!(self, Body::QuadrupedLow(_)) { + 0.7 + } else { + 1.0 + }; + cd * PI * radius.x * radius.z + }, + + Body::BirdMedium(bird) => bird.parasite_drag_coefficient(), + Body::BirdLarge(bird) => bird.parasite_drag_coefficient(), + Body::Dragon(not_bird) => not_bird.parasite_drag_coefficient(), + + // Cross-section, zero-lift angle; exclude the fins + Body::FishMedium(_) | Body::FishSmall(_) => { + let radius = self.dimensions().map(|a| a * 0.5); + // "A Simple Method to Determine Drag Coefficients in Aquatic Animals", + // D. Bilo and W. Nachtigall, 1980 + const CD: f32 = 0.031; + CD * PI * radius.x * radius.z + }, + + Body::Object(object) => match object { + // very streamlined objects + object::Body::Arrow + | object::Body::ArrowSnake + | object::Body::ArrowTurret + | object::Body::FireworkBlue + | object::Body::FireworkGreen + | object::Body::FireworkPurple + | object::Body::FireworkRed + | object::Body::FireworkWhite + | object::Body::FireworkYellow + | object::Body::MultiArrow => { + let radius = self.dimensions().map(|a| a * 0.5); + const CD: f32 = 0.02; + CD * PI * radius.x * radius.z + }, + + // spherical-ish objects + object::Body::BoltFire + | object::Body::BoltFireBig + | object::Body::BoltNature + | object::Body::Bomb + | object::Body::PotionBlue + | object::Body::PotionGreen + | object::Body::PotionRed + | object::Body::Pouch + | object::Body::Pumpkin + | object::Body::Pumpkin2 + | object::Body::Pumpkin3 + | object::Body::Pumpkin4 + | object::Body::Pumpkin5 => { + let radius = self.dimensions().xy().map(|a| a * 0.5); + const CD: f32 = 0.5; + CD * PI * radius.product() + }, + + _ => { + let dim = self.dimensions(); + const CD: f32 = 2.0; + CD * (PI / 6.0 * dim.x * dim.y * dim.z).powf(2.0 / 3.0) + }, + }, + + Body::Ship(_) => { + // Airships tend to use the square of the cube root of its volume for + // reference area + let dim = self.dimensions(); + (PI / 6.0 * dim.x * dim.y * dim.z).powf(2.0 / 3.0) + }, + } + } +} diff --git a/common/src/comp/body/bird_large.rs b/common/src/comp/body/bird_large.rs index 8d282c0480..ea332d0db7 100644 --- a/common/src/comp/body/bird_large.rs +++ b/common/src/comp/body/bird_large.rs @@ -1,6 +1,13 @@ -use crate::{make_case_elim, make_proj_elim}; +use crate::{ + comp::{ + fluid_dynamics::{Drag, WingShape, WingState, Glide}, + Ori, + }, + make_case_elim, make_proj_elim, +}; use rand::{seq::SliceRandom, thread_rng}; use serde::{Deserialize, Serialize}; +use vek::*; make_proj_elim!( body, @@ -23,6 +30,18 @@ impl Body { let body_type = *(&ALL_BODY_TYPES).choose(rng).unwrap(); Self { species, body_type } } + + /// Dimensions of the body (wings folded) + pub const fn dimensions(&self) -> Vec3 { Vec3::new(1.0, 5.0, 2.4) } + + /// Distance from wing tip to wing tip and leading edge to trailing edge + /// respectively + // TODO: Check + pub const fn wing_dimensions(&self) -> Vec2 { Vec2::new(7.0, 1.0) } + + pub fn flying<'a>(&'a self, ori: &'a Ori) -> FlyingBirdLarge<'a> { + FlyingBirdLarge::from((self, ori)) + } } impl From for super::Body { @@ -83,3 +102,55 @@ make_case_elim!( ); pub const ALL_BODY_TYPES: [BodyType; 2] = [BodyType::Female, BodyType::Male]; + +#[derive(Copy, Clone)] +pub struct FlyingBirdLarge<'a> { + wing_shape: WingShape, + wing_state: WingState, + planform_area: f32, + body: &'a Body, + ori: &'a Ori, +} + +impl<'a> From<(&'a Body, &'a Ori)> for FlyingBirdLarge<'a> { + fn from((body, ori): (&'a Body, &'a Ori)) -> Self { + let Vec2 { + x: span_length, + y: chord_length, + } = body.wing_dimensions(); + let planform_area = WingShape::elliptical_planform_area(span_length, chord_length); + FlyingBirdLarge { + wing_shape: WingShape::elliptical(span_length, chord_length), + wing_state: WingState::Flapping, + planform_area, + body, + ori, + } + } +} + +impl Drag for Body { + fn parasite_drag_coefficient(&self) -> f32 { + let radius = self.dimensions().map(|a| a * 0.5); + // "Field Estimates of body::Body Drag Coefficient on the Basis of + // Dives in Passerine Birds", Anders Hedenström and Felix Liechti, 2001 + const CD: f32 = 0.2; + CD * std::f32::consts::PI * radius.x * radius.z + } +} + +impl Drag for FlyingBirdLarge<'_> { + fn parasite_drag_coefficient(&self) -> f32 { + self.body.parasite_drag_coefficient() + self.planform_area * 0.004 + } +} + +impl Glide for FlyingBirdLarge<'_> { + fn wing_shape(&self) -> &WingShape { &self.wing_shape } + + fn is_gliding(&self) -> bool { matches!(self.wing_state, WingState::Fixed) } + + fn planform_area(&self) -> f32 { self.planform_area } + + fn ori(&self) -> &Ori { self.ori } +} diff --git a/common/src/comp/body/bird_medium.rs b/common/src/comp/body/bird_medium.rs index 5eb7dbee65..06dd0cc0d8 100644 --- a/common/src/comp/body/bird_medium.rs +++ b/common/src/comp/body/bird_medium.rs @@ -1,6 +1,13 @@ -use crate::{make_case_elim, make_proj_elim}; +use crate::{ + comp::{ + fluid_dynamics::{Drag, WingShape, WingState, Glide}, + Ori, + }, + make_case_elim, make_proj_elim, +}; use rand::{seq::SliceRandom, thread_rng}; use serde::{Deserialize, Serialize}; +use vek::*; make_proj_elim!( body, @@ -23,6 +30,18 @@ impl Body { let body_type = *(&ALL_BODY_TYPES).choose(rng).unwrap(); Self { species, body_type } } + + /// Dimensions of the body (wings folded) + pub const fn dimensions(&self) -> Vec3 { Vec3::new(0.5, 1.0, 1.1) } + + /// Distance from wing tip to wing tip and leading edge to trailing edge + /// respectively + // TODO: Check + pub const fn wing_dimensions(&self) -> Vec2 { Vec2::new(2.0, 0.4) } + + pub fn flying<'a>(&'a self, ori: &'a Ori) -> FlyingBirdMedium<'a> { + FlyingBirdMedium::from((self, ori)) + } } impl From for super::Body { @@ -102,3 +121,55 @@ make_case_elim!( } ); pub const ALL_BODY_TYPES: [BodyType; 2] = [BodyType::Female, BodyType::Male]; + +#[derive(Copy, Clone)] +pub struct FlyingBirdMedium<'a> { + wing_shape: WingShape, + wing_state: WingState, + planform_area: f32, + body: &'a Body, + ori: &'a Ori, +} + +impl<'a> From<(&'a Body, &'a Ori)> for FlyingBirdMedium<'a> { + fn from((body, ori): (&'a Body, &'a Ori)) -> Self { + let Vec2 { + x: span_length, + y: chord_length, + } = body.wing_dimensions(); + let planform_area = WingShape::elliptical_planform_area(span_length, chord_length); + FlyingBirdMedium { + wing_shape: WingShape::elliptical(span_length, chord_length), + wing_state: WingState::Flapping, + planform_area, + body, + ori, + } + } +} + +impl Drag for Body { + fn parasite_drag_coefficient(&self) -> f32 { + let radius = self.dimensions().map(|a| a * 0.5); + // "Field Estimates of body::Body Drag Coefficient on the Basis of + // Dives in Passerine Birds", Anders Hedenström and Felix Liechti, 2001 + const CD: f32 = 0.2; + CD * std::f32::consts::PI * radius.x * radius.z + } +} + +impl Drag for FlyingBirdMedium<'_> { + fn parasite_drag_coefficient(&self) -> f32 { + self.body.parasite_drag_coefficient() + self.planform_area * 0.004 + } +} + +impl Glide for FlyingBirdMedium<'_> { + fn wing_shape(&self) -> &WingShape { &self.wing_shape } + + fn is_gliding(&self) -> bool { matches!(self.wing_state, WingState::Fixed) } + + fn planform_area(&self) -> f32 { self.planform_area } + + fn ori(&self) -> &Ori { self.ori } +} diff --git a/common/src/comp/body/dragon.rs b/common/src/comp/body/dragon.rs index 07c3991a44..e86f60711b 100644 --- a/common/src/comp/body/dragon.rs +++ b/common/src/comp/body/dragon.rs @@ -1,6 +1,13 @@ -use crate::{make_case_elim, make_proj_elim}; +use crate::{ + comp::{ + fluid_dynamics::{Drag, Glide, WingShape, WingState}, + Ori, + }, + make_case_elim, make_proj_elim, +}; use rand::{seq::SliceRandom, thread_rng}; use serde::{Deserialize, Serialize}; +use vek::*; make_proj_elim!( body, @@ -23,6 +30,18 @@ impl Body { let body_type = *(&ALL_BODY_TYPES).choose(rng).unwrap(); Self { species, body_type } } + + /// Dimensions of the body (wings folded) + pub const fn dimensions(&self) -> Vec3 { Vec3::new(5.0, 10.0, 16.0) } + + /// Distance from wing tip to wing tip and leading edge to trailing edge + /// respectively + // TODO: Check + pub const fn wing_dimensions(&self) -> Vec2 { Vec2::new(16.0, 5.0) } + + pub fn flying<'a>(&'a self, ori: &'a Ori) -> FlyingDragon<'a> { + FlyingDragon::from((self, ori)) + } } impl From for super::Body { @@ -77,3 +96,55 @@ make_case_elim!( ); pub const ALL_BODY_TYPES: [BodyType; 2] = [BodyType::Female, BodyType::Male]; + +#[derive(Copy, Clone)] +pub struct FlyingDragon<'a> { + wing_shape: WingShape, + wing_state: WingState, + planform_area: f32, + body: &'a Body, + ori: &'a Ori, +} + +impl<'a> From<(&'a Body, &'a Ori)> for FlyingDragon<'a> { + fn from((body, ori): (&'a Body, &'a Ori)) -> Self { + let Vec2 { + x: span_length, + y: chord_length, + } = body.wing_dimensions(); + let planform_area = WingShape::elliptical_planform_area(span_length, chord_length); + FlyingDragon { + wing_shape: WingShape::elliptical(span_length, chord_length), + wing_state: WingState::Flapping, + planform_area, + body, + ori, + } + } +} + +impl Drag for Body { + fn parasite_drag_coefficient(&self) -> f32 { + let radius = self.dimensions().map(|a| a * 0.5); + // "Field Estimates of body::Body Drag Coefficient on the Basis of + // Dives in Passerine Birds", Anders Hedenström and Felix Liechti, 2001 + const CD: f32 = 0.2; + CD * std::f32::consts::PI * radius.x * radius.z + } +} + +impl Drag for FlyingDragon<'_> { + fn parasite_drag_coefficient(&self) -> f32 { + self.body.parasite_drag_coefficient() + self.planform_area * 0.004 + } +} + +impl Glide for FlyingDragon<'_> { + fn wing_shape(&self) -> &WingShape { &self.wing_shape } + + fn is_gliding(&self) -> bool { matches!(self.wing_state, WingState::Fixed) } + + fn planform_area(&self) -> f32 { self.planform_area } + + fn ori(&self) -> &Ori { self.ori } +} diff --git a/common/src/comp/fluid_dynamics.rs b/common/src/comp/fluid_dynamics.rs index 70632ea387..2f2335a4ed 100644 --- a/common/src/comp/fluid_dynamics.rs +++ b/common/src/comp/fluid_dynamics.rs @@ -1,5 +1,3 @@ -#[cfg(not(target_arch = "wasm32"))] -use super::body::{object, Body}; use super::{Density, Ori, Vel}; use crate::{ consts::{AIR_DENSITY, LAVA_DENSITY, WATER_DENSITY}, @@ -120,20 +118,12 @@ impl Default for Fluid { } } -pub struct Wings { - pub aspect_ratio: f32, - pub planform_area: f32, - pub ori: Ori, -} +pub trait Drag: Clone { + /// Drag coefficient for skin friction and flow separation. + /// (Multiplied by reference area.) + fn parasite_drag_coefficient(&self) -> f32; -#[cfg(not(target_arch = "wasm32"))] -impl Body { - pub fn aerodynamic_forces( - &self, - rel_flow: &Vel, - fluid_density: f32, - wings: Option<&Wings>, - ) -> Vec3 { + fn drag(&self, rel_flow: &Vel, fluid_density: f32) -> Vec3 { let v_sq = rel_flow.0.magnitude_squared(); if v_sq < 0.25 { // don't bother with miniscule forces @@ -141,181 +131,100 @@ impl Body { } else { let rel_flow_dir = Dir::new(rel_flow.0 / v_sq.sqrt()); // All the coefficients come pre-multiplied by their reference area - 0.5 * fluid_density - * v_sq - * match wings { - Some(&Wings { - aspect_ratio, - planform_area, - ori, - }) => { - if aspect_ratio > 25.0 { - tracing::warn!( - "Calculating lift for wings with an aspect ratio of {}. The \ - formulas are only valid for aspect ratios below 25.", - aspect_ratio - ) - }; - let ar = aspect_ratio.min(24.0); - // We have an elliptical wing; proceed to calculate its lift and drag - - // aoa will be positive when we're pitched up and negative otherwise - 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 = lift_coefficient(ar, planform_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) - 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 - let aoa_eff = aoa - aoa_i; - // 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() - }; - - // drag coefficient - let c_d = { - // Oswald's efficiency factor (empirically derived--very magical) - // (this definition should not be used for aspect ratios > 25) - let e = 1.78 * (1.0 - 0.045 * ar.powf(0.68)) - 0.64; - // induced drag coefficient (drag due to lift) - let cdi = c_l.powi(2) / (PI * e * ar); - - zero_lift_drag_coefficient(planform_area) - + self.parasite_drag_coefficient() - + cdi - }; - debug_assert!(c_d.is_sign_positive()); - debug_assert!(c_l.is_sign_positive() || aoa.is_sign_negative()); - - c_l * *lift_dir + c_d * *rel_flow_dir - }, - - _ => self.parasite_drag_coefficient() * *rel_flow_dir, - } + 0.5 * fluid_density * v_sq * self.parasite_drag_coefficient() * *rel_flow_dir } } +} - /// Parasite drag is the sum of pressure drag and skin friction. - /// Skin friction is the drag arising from the shear forces between a fluid - /// and a surface, while pressure drag is due to flow separation. Both are - /// viscous effects. - /// - /// This coefficient includes the reference area. - fn parasite_drag_coefficient(&self) -> f32 { - // Reference area and drag coefficient assumes best-case scenario of the - // orientation producing least amount of drag - match self { - // Cross-section, head/feet first - Body::BipedLarge(_) | Body::BipedSmall(_) | Body::Golem(_) | Body::Humanoid(_) => { - let dim = self.dimensions().xy().map(|a| 0.7 * a * 0.5); - const CD: f32 = 0.7; - CD * PI * dim.x * dim.y - }, +pub trait Glide: Drag { + const STALL_ANGLE: f32 = PI * 0.1; - // Cross-section, nose/tail first - Body::Theropod(_) - | Body::QuadrupedMedium(_) - | Body::QuadrupedSmall(_) - | Body::QuadrupedLow(_) => { - let dim = self.dimensions().map(|a| a * 0.5); - let cd: f32 = if matches!(self, Body::QuadrupedLow(_)) { - 0.7 + fn wing_shape(&self) -> &WingShape; + + fn planform_area(&self) -> f32; + + fn ori(&self) -> &Ori; + + fn is_gliding(&self) -> bool; + + /// Total lift coefficient for a finite wing of symmetric aerofoil shape and + /// elliptical pressure distribution. (Multiplied by reference area.) + fn lift_coefficient(&self, angle_of_attack: f32) -> f32 { + let aoa_abs = angle_of_attack.abs(); + self.planform_area() + * if self.is_gliding() { 1.0 } else { 0.2 } + * if aoa_abs <= Self::STALL_ANGLE { + self.wing_shape().lift_slope() * angle_of_attack + } else { + // 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 = angle_of_attack.signum(); + let c_l_max = self.wing_shape().lift_slope() * Self::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 { - 1.0 - }; - cd * PI * dim.x * dim.z - }, + // 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 + } + } + } - // 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); - let cd: f32 = match self { - // "Field Estimates of Body Drag Coefficient on the Basis of Dives in Passerine - // Birds", Anders Hedenström and Felix Liechti, 2001 - Body::BirdLarge(_) | Body::BirdMedium(_) => 0.2, - // arbitrary - _ => 0.7, - }; - cd * PI * dim.x * 0.2 * dim.z - }, + /// Total drag coefficient (multiplied by reference area) + fn drag_coefficient(&self, aspect_ratio: f32, lift_coefficient: f32) -> f32 { + self.parasite_drag_coefficient() + induced_drag_coefficient(aspect_ratio, lift_coefficient) + } - // Cross-section, zero-lift angle; exclude the fins (width * 0.2) - Body::FishMedium(_) | Body::FishSmall(_) => { - let dim = self.dimensions().map(|a| a * 0.5); - // "A Simple Method to Determine Drag Coefficients in Aquatic Animals", - // D. Bilo and W. Nachtigall, 1980 - const CD: f32 = 0.031; - CD * PI * dim.x * 0.2 * dim.z - }, + fn aerodynamic_forces(&self, rel_flow: &Vel, fluid_density: f32) -> Vec3 { + let v_sq = rel_flow.0.magnitude_squared(); + if v_sq < 0.1 { + // don't bother with miniscule forces + Vec3::zero() + } else { + let q = 0.5 * fluid_density * v_sq; + let rel_flow_dir = Dir::new(rel_flow.0 / v_sq.sqrt()); + if self.is_gliding() { + let ori = self.ori(); + let ar = self.wing_shape().aspect_ratio(); + // aoa will be positive when we're pitched up and negative otherwise + 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 = self.lift_coefficient(aoa); - Body::Object(object) => match object { - // very streamlined objects - object::Body::Arrow - | object::Body::ArrowSnake - | object::Body::ArrowTurret - | object::Body::FireworkBlue - | object::Body::FireworkGreen - | object::Body::FireworkPurple - | object::Body::FireworkRed - | object::Body::FireworkWhite - | object::Body::FireworkYellow - | object::Body::MultiArrow => { - let dim = self.dimensions().map(|a| a * 0.5); - const CD: f32 = 0.02; - CD * PI * dim.x * dim.z - }, + let lift = q * c_l * *lift_dir(ori, ar, c_l, aoa); + let drag = q * self.drag_coefficient(ar, c_l) * *rel_flow_dir; - // spherical-ish objects - object::Body::BoltFire - | object::Body::BoltFireBig - | object::Body::BoltNature - | object::Body::Bomb - | object::Body::PotionBlue - | object::Body::PotionGreen - | object::Body::PotionRed - | object::Body::Pouch - | object::Body::Pumpkin - | object::Body::Pumpkin2 - | object::Body::Pumpkin3 - | object::Body::Pumpkin4 - | object::Body::Pumpkin5 => { - let dim = self.dimensions().map(|a| a * 0.5); - const CD: f32 = 0.5; - CD * PI * dim.x * dim.z - }, - - _ => { - let dim = self.dimensions(); - const CD: f32 = 2.0; - CD * (PI / 6.0 * dim.x * dim.y * dim.z).powf(2.0 / 3.0) - }, - }, - - Body::Ship(_) => { - // Airships tend to use the square of the cube root of its volume for - // reference area - let dim = self.dimensions(); - (PI / 6.0 * dim.x * dim.y * dim.z).powf(2.0 / 3.0) - }, + lift + drag + } else { + q * self.parasite_drag_coefficient() * *rel_flow_dir + } } } } + +pub fn lift_dir(ori: &Ori, aspect_ratio: f32, lift_coefficient: f32, angle_of_attack: f32) -> Dir { + // 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) + + // induced angle of attack + let aoa_i = lift_coefficient / (PI * aspect_ratio); + // effective angle of attack; the aoa as seen by aerofoil after downwash + let aoa_eff = angle_of_attack - aoa_i; + + // 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() +} + /// Geometric angle of attack /// /// # Note @@ -329,68 +238,99 @@ pub fn angle_of_attack(ori: &Ori, rel_flow_dir: &Dir) -> f32 { .unwrap_or(0.0) } -/// 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; - planform_area - * if aoa_abs < stall_angle { - lift_slope(aspect_ratio, None) * aoa - } else { - // 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 { - // 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 - } - } +#[derive(Copy, Clone)] +pub enum WingState { + Fixed, + Flapping, + Retracted, } -/// 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 { planform_area * 0.004 } +#[derive(Copy, Clone, Debug, PartialEq, Serialize, Deserialize)] +pub enum WingShape { + Elliptical { aspect_ratio: f32 }, + // Tapered { aspect_ratio: f32, e: f32 }, + Swept { aspect_ratio: f32, angle: f32 }, + // Delta, +} -/// 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 -// TODO: Look into handling tapered wings -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))) +impl WingShape { + pub fn aspect_ratio(&self) -> f32 { + match self { + Self::Elliptical { aspect_ratio } => *aspect_ratio, + Self::Swept { aspect_ratio, .. } => *aspect_ratio, + } + } + + pub fn elliptical_planform_area(span_length: f32, chord_length: f32) -> f32 { + std::f32::consts::PI * span_length * chord_length * 0.25 + } + + pub fn elliptical(span_length: f32, chord_length: f32) -> Self { + Self::Elliptical { + aspect_ratio: span_length.powi(2) + / Self::elliptical_planform_area(span_length, chord_length), + } + } + + /// 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 + pub fn lift_slope(&self) -> f32 { + // lift slope for a thin aerofoil, given by Thin Aerofoil Theory + let a0 = 2.0 * PI; + match self { + WingShape::Elliptical { aspect_ratio } => { + 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))) + } + }, + // WingShape::Tapered { aspect_ratio, e } => todo!(), + WingShape::Swept { + aspect_ratio, + angle, + } => { + // for swept wings we use Kuchemann's modification to Helmbold's + // equation + let a0_cos_sweep = a0 * angle.cos(); + let x = a0_cos_sweep / (PI * aspect_ratio); + a0_cos_sweep / ((1.0 + x.powi(2)).sqrt() + x) + }, + } } } + +/// Induced drag coefficient (drag due to lift) +pub fn induced_drag_coefficient(aspect_ratio: f32, lift_coefficient: f32) -> f32 { + let ar = aspect_ratio; + if ar > 25.0 { + tracing::warn!( + "Calculating induced drag for wings with a given aspect ratio of {}. The formulas are \ + only valid for aspect ratios below 25, so it will be substituted.", + ar + ) + }; + let ar = ar.min(24.0); + // Oswald's efficiency factor (empirically derived--very magical) + // (this definition should not be used for aspect ratios > 25) + let e = 1.78 * (1.0 - 0.045 * ar.powf(0.68)) - 0.64; + // induced drag coefficient (drag due to lift) + lift_coefficient.powi(2) / (PI * e * ar) +} diff --git a/common/src/consts.rs b/common/src/consts.rs index 59b7dbd8bf..4cbaa9598e 100644 --- a/common/src/consts.rs +++ b/common/src/consts.rs @@ -2,7 +2,7 @@ pub const MAX_PICKUP_RANGE: f32 = 8.0; pub const MAX_MOUNT_RANGE: f32 = 14.0; -pub const GRAVITY: f32 = 9.8; +pub const GRAVITY: f32 = 25.0; pub const FRIC_GROUND: f32 = 0.15; // Values for air taken from http://www-mdp.eng.cam.ac.uk/web/library/enginfo/aerothermal_dvd_only/aero/atmos/atmos.html diff --git a/common/src/lib.rs b/common/src/lib.rs index ec08f93737..e554b84456 100644 --- a/common/src/lib.rs +++ b/common/src/lib.rs @@ -6,6 +6,7 @@ #![feature( arbitrary_enum_discriminant, associated_type_defaults, + bindings_after_at, bool_to_option, const_generics, fundamental, diff --git a/common/src/states/glide.rs b/common/src/states/glide.rs index 72aed78390..0ba59a47f6 100644 --- a/common/src/states/glide.rs +++ b/common/src/states/glide.rs @@ -1,8 +1,9 @@ use super::utils::handle_climb; use crate::{ comp::{ - fluid_dynamics::angle_of_attack, inventory::slot::EquipSlot, CharacterState, Ori, - StateUpdate, Vel, + fluid_dynamics::{angle_of_attack, Drag, Glide, WingShape}, + inventory::slot::EquipSlot, + CharacterState, Ori, StateUpdate, Vel, }, states::behavior::{CharacterBehavior, JoinData}, util::{Dir, Plane, Projection}, @@ -17,7 +18,7 @@ const PITCH_SLOW_TIME: f32 = 0.5; pub struct Data { /// The aspect ratio is the ratio of the span squared to actual planform /// area - pub aspect_ratio: f32, + pub wing_shape: WingShape, pub planform_area: f32, pub ori: Ori, last_vel: Vel, @@ -36,7 +37,7 @@ impl Data { let planform_area = PI * chord_length * span_length * 0.25; let aspect_ratio = span_length.powi(2) / planform_area; Self { - aspect_ratio, + wing_shape: WingShape::Elliptical { aspect_ratio }, planform_area, ori, last_vel: Vel::zero(), @@ -69,6 +70,20 @@ impl Data { } } +impl Drag for Data { + fn parasite_drag_coefficient(&self) -> f32 { self.planform_area * 0.004 } +} + +impl Glide for Data { + fn wing_shape(&self) -> &WingShape { &self.wing_shape } + + fn is_gliding(&self) -> bool { true } + + fn planform_area(&self) -> f32 { self.planform_area } + + fn ori(&self) -> &Ori { &self.ori } +} + impl CharacterBehavior for Data { fn behavior(&self, data: &JoinData) -> StateUpdate { let mut update = StateUpdate::from(data); diff --git a/common/systems/src/phys.rs b/common/systems/src/phys.rs index 81d1cec1c5..eda32efbb6 100644 --- a/common/systems/src/phys.rs +++ b/common/systems/src/phys.rs @@ -1,7 +1,7 @@ use common::{ comp::{ body::ship::figuredata::{VoxelCollider, VOXEL_COLLIDER_MANIFEST}, - fluid_dynamics::{Fluid, LiquidKind, Wings}, + fluid_dynamics::{Fluid, LiquidKind, Glide}, BeamSegment, Body, CharacterState, Collider, Density, Mass, Mounting, Ori, PhysicsState, Pos, PosVelDefer, PreviousPhysCache, Projectile, Scale, Shockwave, Stats, Sticky, Vel, }, @@ -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,9 @@ fn fluid_density(height: f32, fluid: &Fluid) -> Density { fn integrate_forces( dt: &DeltaTime, mut vel: Vel, - (body, wings): (&Body, Option<&Wings>), + ori: Option<&Ori>, + body: &Body, + character_state: Option<&CharacterState>, density: &Density, mass: &Mass, fluid: &Fluid, @@ -59,10 +60,18 @@ 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 impulse = dt.0 * body.aerodynamic_forces(&rel_flow, fluid_density.0, wings); - debug_assert!(!impulse.map(|a| a.is_nan()).reduce_or()); - if !impulse.is_approx_zero() { - let new_v = vel.0 + impulse / mass.0; + let aerodynamic_forces = body.aerodynamic_forces(ori, &rel_flow, fluid_density.0) + + match character_state { + Some(CharacterState::Glide(glider)) => { + glider.aerodynamic_forces(&rel_flow, fluid_density.0) + }, + _ => Vec3::zero(), + }; + // let impulse = dt.0 * body.aerodynamic_forces(ori, &rel_flow, + // fluid_density.0); + debug_assert!(!aerodynamic_forces.map(|a| a.is_nan()).reduce_or()); + if !aerodynamic_forces.is_approx_zero() { + let new_v = vel.0 + dt.0 * aerodynamic_forces / mass.0; // If the new velocity is in the opposite direction, it's because the forces // involved are too high for the current tick to handle. We deal with this by // removing the component of our velocity vector along the direction of force. @@ -71,7 +80,7 @@ fn integrate_forces( if new_v.dot(vel.0) < 0.0 { // Multiply by a factor to prevent full stop, as this can cause things to get // stuck in high-density medium - vel.0 -= vel.0.projected(&impulse) * 0.9; + vel.0 -= vel.0.projected(&aerodynamic_forces) * 0.9; } else { vel.0 = new_v; } @@ -81,8 +90,7 @@ fn integrate_forces( // Hydrostatic/aerostatic forces // modify gravity to account for the effective density as a result of buoyancy - let down_force = dt.0 * gravity * (density.0 - fluid_density.0) / density.0; - vel.0.z -= down_force; + vel.0.z -= dt.0 * gravity * (density.0 - fluid_density.0) / density.0; vel } @@ -576,6 +584,7 @@ impl<'a> PhysicsData<'a> { read.character_states.maybe(), &write.physics_states, &read.masses, + write.orientations.maybe(), &read.densities, !&read.mountings, ) @@ -594,6 +603,7 @@ impl<'a> PhysicsData<'a> { character_state, physics_state, mass, + ori, density, _, )| { @@ -617,24 +627,12 @@ impl<'a> PhysicsData<'a> { vel.0.z -= dt.0 * GRAVITY; }, Some(fluid) => { - let wings = match character_state { - Some(&CharacterState::Glide(states::glide::Data { - aspect_ratio, - planform_area, - ori, - .. - })) => Some(Wings { - aspect_ratio, - planform_area, - ori, - }), - - _ => None, - }; vel.0 = integrate_forces( &dt, *vel, - (body, wings.as_ref()), + ori, + body, + character_state, density, mass, &fluid,