From 7b4bb2de99b1764f72286e3310c9293e614f2818 Mon Sep 17 00:00:00 2001 From: Joshua Barretto Date: Thu, 18 May 2023 17:17:22 +0100 Subject: [PATCH] Use non-physical mode for liquid drag --- common/src/comp/fluid_dynamics.rs | 77 +++++++++++++++++++++++++++++++ common/src/states/utils.rs | 35 +++++++------- common/systems/src/phys.rs | 54 ++++++++++++++-------- 3 files changed, 131 insertions(+), 35 deletions(-) diff --git a/common/src/comp/fluid_dynamics.rs b/common/src/comp/fluid_dynamics.rs index c7b3558497..c4edad4d87 100644 --- a/common/src/comp/fluid_dynamics.rs +++ b/common/src/comp/fluid_dynamics.rs @@ -211,6 +211,83 @@ impl Body { } } + /// Physically incorrect (but relatively dt-independent) way to calculate drag coefficients for liquids. + // TODO: Remove this in favour of `aerodynamic_forces` (see: `integrate_forces`) + pub fn drag_coefficient_liquid( + &self, + rel_flow: &Vel, + fluid_density: f32, + wings: Option<&Wings>, + scale: f32, + ) -> Vec3 { + let v_sq = rel_flow.0.magnitude_squared(); + let rel_flow_dir = Dir::new(rel_flow.0 / v_sq.sqrt()); + 0.5 * fluid_density + * 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, 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() + }; + + // induced drag coefficient (drag due to lift) + let cdi = { + // 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; + c_l.powi(2) / (PI * e * ar) + }; + + // drag coefficient + let c_d = zero_lift_drag_coefficient() + cdi; + debug_assert!(c_d.is_sign_positive()); + debug_assert!(c_l.is_sign_positive() || aoa.is_sign_negative()); + + planform_area * scale.powf(2.0) * (c_l * *lift_dir + c_d * *rel_flow_dir) + + self.parasite_drag(scale).powf(0.75) * *rel_flow_dir + }, + + _ => self.parasite_drag(scale).powf(0.75) * *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 diff --git a/common/src/states/utils.rs b/common/src/states/utils.rs index 5cd49ec014..eaefc3e033 100644 --- a/common/src/states/utils.rs +++ b/common/src/states/utils.rs @@ -220,7 +220,7 @@ impl Body { quadruped_low::Species::Mossdrake => 1.7, _ => 2.0, }, - Body::Ship(ship) if ship.has_water_thrust() => 0.2, + Body::Ship(ship) if ship.has_water_thrust() => 0.35, Body::Ship(_) => 0.12, Body::Arthropod(_) => 3.5, } @@ -236,16 +236,16 @@ impl Body { match self { Body::Object(_) => return None, Body::ItemDrop(_) => return None, - Body::Ship(ship) if ship.has_water_thrust() => 0.2 * self.mass().0, + Body::Ship(ship) if ship.has_water_thrust() => 350.0 * self.mass().0, Body::Ship(_) => return None, Body::BipedLarge(_) => 120.0 * self.mass().0, - Body::Golem(_) => 0.5 * self.mass().0, + Body::Golem(_) => 100.0 * self.mass().0, Body::BipedSmall(_) => 1000.0 * self.mass().0, - Body::BirdMedium(_) => 1500.0 * self.mass().0, - Body::BirdLarge(_) => 35.0 * self.mass().0, - Body::FishMedium(_) => 75.0 * self.mass().0, - Body::FishSmall(_) => 120.0 * self.mass().0, - Body::Dragon(_) => 0.5 * self.mass().0, + Body::BirdMedium(_) => 400.0 * self.mass().0, + Body::BirdLarge(_) => 400.0 * self.mass().0, + Body::FishMedium(_) => 900.0 * self.mass().0, + Body::FishSmall(_) => 1000.0 * self.mass().0, + Body::Dragon(_) => 50.0 * self.mass().0, // Humanoids are a bit different: we try to give them thrusts that result in similar // speeds for gameplay reasons Body::Humanoid(_) => 4_000_000.0 / self.mass().0, @@ -256,16 +256,16 @@ impl Body { | theropod::Species::Woodraptor | theropod::Species::Dodarock | theropod::Species::Axebeak - | theropod::Species::Yale => 700.0 * self.mass().0, - _ => 3.0 * self.mass().0, + | theropod::Species::Yale => 500.0 * self.mass().0, + _ => 150.0 * self.mass().0, }, - Body::QuadrupedLow(_) => 14.0 * self.mass().0, + Body::QuadrupedLow(_) => 1200.0 * self.mass().0, Body::QuadrupedMedium(body) => match body.species { - quadruped_medium::Species::Mammoth => 75.0 * self.mass().0, - _ => 500.0 * self.mass().0, + quadruped_medium::Species::Mammoth => 150.0 * self.mass().0, + _ => 1000.0 * self.mass().0, }, Body::QuadrupedSmall(_) => 1500.0 * self.mass().0, - Body::Arthropod(_) => 300.0 * self.mass().0, + Body::Arthropod(_) => 500.0 * self.mass().0, } * front_profile, ) } @@ -671,14 +671,15 @@ fn swim_move( // Assume that feet/flippers get less efficient as we become less submerged let move_z = move_z.min((submersion * 1.5 - 0.5).clamp(0.0, 1.0).powi(2)); - update.vel.0 += Vec3::broadcast(data.dt.0) - * Vec3::new(dir.x, dir.y, move_z) + update.vel.0 += Vec3::new(dir.x, dir.y, move_z) // TODO: Should probably be normalised, but creates odd discrepancies when treading water // .try_normalized() // .unwrap_or_default() * water_accel // Gives a good balance between submerged and surface speed - * (submersion * 3.0).clamp(0.0, 1.0); + * submersion.clamp(0.0, 1.0).sqrt() + // Good approximate compensation for dt-dependent effects + * data.dt.0 * 0.04; true } else { diff --git a/common/systems/src/phys.rs b/common/systems/src/phys.rs index d3ded17a6f..570329d612 100644 --- a/common/systems/src/phys.rs +++ b/common/systems/src/phys.rs @@ -62,29 +62,45 @@ 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( + // HACK: We should really use the latter logic (i.e: `aerodynamic_forces`) for liquids, but it results in + // pretty serious dt-dependent problems that are extremely difficult to resolve. This is a compromise: for + // liquids only, we calculate drag using an incorrect (but still visually plausible) model that is much more + // resistant to differences in dt. Because it's physically incorrect anyway, there are magic coefficients that + // exist simply to get us closer to what water 'should' feel like. + if fluid.is_liquid() { + let fric = body.drag_coefficient_liquid( &rel_flow, fluid_density.0, wings, scale.map_or(1.0, |s| s.0), - ); - debug_assert!(!impulse.map(|a| a.is_nan()).reduce_or()); - if !impulse.is_approx_zero() { - let new_v = vel.0 + impulse / 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. - // This way we can only ever lose velocity and will never experience a reverse - // in direction from events such as falling into water at high velocities. - 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; - } else { - vel.0 = new_v; + ).magnitude() * 0.02; + + vel.0 *= (1.0 / (1.0 + fric)).powf(dt.0 * 10.0); + } else { + let impulse = dt.0 + * body.aerodynamic_forces( + &rel_flow, + fluid_density.0, + wings, + scale.map_or(1.0, |s| s.0), + ); + debug_assert!(!impulse.map(|a| a.is_nan()).reduce_or()); + if !impulse.is_approx_zero() { + let new_v = vel.0 + impulse / 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. + // This way we can only ever lose velocity and will never experience a reverse + // in direction from events such as falling into water at high velocities. + 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; + } else { + vel.0 = new_v; + } } - }; + } debug_assert!(!vel.0.map(|a| a.is_nan()).reduce_or()); }; @@ -1732,6 +1748,8 @@ fn box_voxel_collision + ReadVol>( (old_depth + old_pos.z - pos.0.z).max(depth) }); + physics_state.ground_vel = ground_vel; + Fluid::Liquid { kind, depth: new_depth,