From 85d4fb828a687f6880016466df9952a89f4f0995 Mon Sep 17 00:00:00 2001 From: grillkol Date: Fri, 23 Oct 2026 12:43:57 +0200 Subject: [PATCH] Add American V8 preset, carb fuel transport, and adaptive desktop audio --- CARB_PROGRESS.md | 63 +++++++++++++++ CMakeLists.txt | 1 + README.md | 106 ++++++++++++++++++++++--- include/FluidSim/Atmosphere.hpp | 1 + include/FluidSim/AudioOutputStream.hpp | 5 ++ include/FluidSim/Carburetor.hpp | 46 +++++++++++ include/FluidSim/Component.hpp | 2 + include/FluidSim/Crankcase.hpp | 9 ++- include/FluidSim/DesktopScenario.hpp | 5 +- include/FluidSim/Ecu.hpp | 5 +- include/FluidSim/Engine.hpp | 2 + include/FluidSim/Fluid.hpp | 83 +++++++++++++------ include/FluidSim/OutdoorAudio.hpp | 6 +- include/FluidSim/Scenario.hpp | 26 +++--- include/FluidSim/SoundEngine.hpp | 2 + include/FluidSim/TwoStrokeCylinder.hpp | 4 + include/FluidSim/TwoStrokeScenario.hpp | 3 + src/audio_output_stream.cpp | 29 ++++++- src/carburetor.cpp | 102 ++++++++++++++++++++++++ src/crankcase.cpp | 12 ++- src/desktop_main.cpp | 75 +++++++++++------ src/ecu.cpp | 8 +- src/engine.cpp | 40 +++++++--- src/main.cpp | 8 +- src/outdoor_audio.cpp | 8 +- src/scenario.cpp | 98 ++++++++++++++++++++--- src/turbocharger.cpp | 1 + src/two_stroke_cylinder.cpp | 23 ++++-- src/two_stroke_scenario.cpp | 32 +++++++- tests/fluid_validation.cpp | 16 +++- tests/two_stroke_physics.cpp | 22 +++++ tests/two_stroke_regression.cpp | 22 ++++- tools/audio_capture.cpp | 4 +- tools/scenario_benchmark.cpp | 32 ++++++-- tools/two_stroke_benchmark.cpp | 6 +- verification/YZ125_BASELINE.md | 42 ++++++---- 36 files changed, 805 insertions(+), 144 deletions(-) create mode 100644 CARB_PROGRESS.md create mode 100644 include/FluidSim/Carburetor.hpp create mode 100644 src/carburetor.cpp diff --git a/CARB_PROGRESS.md b/CARB_PROGRESS.md new file mode 100644 index 0000000..80b3425 --- /dev/null +++ b/CARB_PROGRESS.md @@ -0,0 +1,63 @@ +# Carburetor work in progress + +Resumed on 2026-10-23 to fix throttle-release crashes and remove scheduled pops. +All changes are uncommitted. The historical results below predate that revision. + +Current revision: bounded slide aperture to avoid negative roundoff near its +seat; removed the hardcoded firing intervals and forced late spark. Exhaust +combustion was subsequently removed at the user's request due to lag and limited +audible benefit. Fuel remains transported through pipes without reacting. +Cylinder combustion consumes air in proportion to burned fuel. +Warm-running startup now seeds fuel in the intake, crankcase and cylinder. +Pilot effective diameter is .40 mm with .08 base/.3 maximum idle slide opening. +The revised transient passes: peak 12954 RPM, idle 2583 RPM, peak audio .0985. + +## Implemented + +- Added `Carburetor`, with effective pilot, needle and main jet restrictions. +- Jet flow uses simulated air mass flow, air density, venturi dynamic pressure, + static pressure drop, fuel density, discharge coefficients and float-bowl + supply. It is not a cylinder-side target-AFR calculation. +- Added wall-film deposition and evaporation during throttle transients. +- Added fuel as a conserved transported species to pipe cells, boundaries, + 0-D volumes, crankcase, cylinder and turbo internal transfers. +- Carb fuel enters the intake plenum, then moves through the reed, crankcase, + transfer pipe and cylinder with the simulated gas flow. +- Cylinder heat release is limited by transported fuel, trapped oxygen and a + basic 7:1--24:1 flammability window. Rich excess remains unburned. +- Disabled electronic overrun fuel cut for the two-stroke. A deterministic + sparse late-fire cadence overlaps exhaust-port opening to make separated + coast-down bangs through the pressure/audio solver; no sound sample is mixed. +- Added fuel conservation and downstream transport tests, carb bounds/cadence + checks, transient firing-count checks, benchmark AFR columns and desktop carb + labeling. + +## Last verified state + +- Before the final throttle-slide geometry change, the entire 10-test CTest + suite passed. +- That transient reached about 15,954 RPM, recovered to about 3,151 RPM, and + produced 11 late firings in the first second after release with no clipping + or pipe positivity repairs. +- After replacing the old butterfly-area formula with a circular-segment slide + aperture and reducing the idle slide calibration, the project and desktop + rebuilt successfully and `fluidsim_two_stroke_physics` passed. + +## Current task / next step + +The `fluidsim_two_stroke_regression` run after the slide-area change was stopped +mid-run when work was paused. Run it first, then run full CTest. If idle recovery, +mixture, or overrun firing count moved, tune `idle_base_opening`, +`idle_max_opening`, and the effective jet dimensions in `Carburetor.hpp`. +Re-run `fluidsim_two_stroke_benchmark` afterward and update the YZ125 verification +numbers if the slide geometry changed them. + +The most relevant files are: + +- `include/FluidSim/Carburetor.hpp` +- `src/carburetor.cpp` +- `src/two_stroke_scenario.cpp` +- `src/two_stroke_cylinder.cpp` +- `include/FluidSim/Fluid.hpp` +- `tests/two_stroke_physics.cpp` +- `tests/two_stroke_regression.cpp` diff --git a/CMakeLists.txt b/CMakeLists.txt index e75df6c..e28e46b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -5,6 +5,7 @@ option(FLUIDSIM_BUILD_DESKTOP "Build the SFML windowed simulator with keyboard c set(FLUIDSIM_FONT_FILE "" CACHE FILEPATH "Optional font file copied beside the desktop executable") add_library(fluidsim_core + src/carburetor.cpp src/ecu.cpp src/turbocharger.cpp src/engine.cpp diff --git a/README.md b/README.md index df3eb70..8f326e5 100644 --- a/README.md +++ b/README.md @@ -144,6 +144,59 @@ physical cell diameter and length. Compare exhausts with `--exhaust racing` this option. See [YZ125 verification](verification/YZ125_BASELINE.md) for the geometry, assumptions, speed sweep and remaining unloaded-RPM limitation. +## 1980s American V8 + +Choose **1980s American 5.0 V8** in the desktop menu, or run: + +```sh +./build-desktop/fluidsim_desktop --scenario v8 --scene field +``` + +This is an estimated naturally aspirated 302-class cross-plane V8: eight +101.6 x 76.2 mm cylinders, 7:1 compression, mild pushrod-style valve timing, +1-5-4-2-6-3-7-8 firing order, 750 RPM idle target and 7000 RPM limiter +(6750 RPM recovery). Crank/flywheel inertia is 1.1 kg m², versus the original +.28. Eight unequal .60--.78 m headers feed two independent 1.5-litre bank +collectors: cylinders 1--4 and 5--8. Each collector has its own 50 mm outlet +pipe (1.8 / 1.95 m), and each outlet has a separate spatial audio emitter. A +nominal 3000 kg five-speed car uses the existing clutch/gear controls. +The V8 clutch has 1100 Nm capacity to shorten launch slip with the heavy +flywheel and car; engagement remains progressive rather than forcibly locking. +Final drive is approximately 3.624:1 (6.16 / 1.7), making gearing 1.7x taller +than the previous setup: 1.7x road speed at the same RPM and selected gear. +Throttle opening uses a 65 ms response and 2.3 travels/second limit; closing +uses 160 ms and 1.4 travels/second for a slower release. +The V8 uses the shared ideal EFI fuel dose with ECU throttle/idle control, +ignition, overrun fuel cut and rev limiting. Its earlier carb setup has been +removed. Independent exhaust banks retain their uneven cross-plane pulse +spacing rather than blending both banks in a shared collector. Numerical +diffusion in the pipe solver and audio filtering still limit high-frequency rasp. +Dimensions and ignition/cam calibration are design estimates; +the output is not a validated period engine dyno curve. + +The intake uses one 70 mm inlet duct, one 65 mm throttle body and a shared +4-litre plenum feeding eight 38 mm diameter, 12 cm runners (previously 30 cm). +The revised mesh has 19 pipes and 190 cells, down from 310. Runner cells are +6 cm long and exhaust cells approximately 5--6.5 cm long. This coarsening trades +some high-frequency resolution for more real-time headroom. + +The earlier lift-off slowdown came from the CFL timestep falling below the +22.7 microsecond audio interval: the whole fluid network advanced almost twice +per sample, doubling boundary solves and cell updates. The revised 16-second +rev/coast benchmark stayed at one substep per sample and took .51--.57 wall +seconds per simulated second here, including coast-down (previously up to +1.28). It returned to 759 RPM with finite audio and zero pipe positivity +repairs. The CFL limit remains enforced; other operating points may still +require subdivisions. The desktop retains the 4x playback gain. + +Benchmark or capture with: + +```sh +./build-desktop/fluidsim_scenario_benchmark 8 transient v8 +./build-desktop/fluidsim_scenario_benchmark 20 coast v8 +./build-desktop/fluidsim_audio_capture --scenario v8 --seconds 6 --throttle .3 --out /tmp/v8 +``` + ## Twin with a VZ21-size turbo Choose **600 cc twin + VZ21-size turbo** in the startup menu, or bypass it: @@ -228,26 +281,50 @@ thermodynamics, crankcase and atmosphere state, profiled-pipe geometry, inertive and closed boundaries, gearing and clutch momentum transfer, finite runs for every scenario, and finite audio/reverb output. -## ECU and throttle control +## Fuel and throttle control + +Desktop playback has automatic gain enabled by default for every scenario. +Press **G** to toggle it; the dashboard shows its state. It tracks stereo RMS +and peaks after DC removal, lifting quiet idle audio gradually (1.5-second +gain-rise time constant) and reducing gain faster for loud revs (30 ms). +Extra gain is capped at 12x above the existing fixed 4x playback gain, with +linked stereo peak protection. Silence does not request additional boost and +underruns freeze the detector. Disabling smoothly returns to fixed playback +gain. Simulation pressures and WAV capture levels are unaffected. Both desktop engines use `Ecu` (`include/FluidSim/Ecu.hpp`). The desktop and CLI throttle values are pedal requests; the actual throttle is advanced every simulation sample, independently of rendering. Opening and closing use an 80 ms first-order actuator response and a maximum travel rate of two full -travels per second. An ideal butterfly's projected opening determines orifice -area continuously; shaft obstruction and edge leakage are not modeled. +travels per second. The four-stroke uses an ideal butterfly projection; the +two-stroke uses the circular-segment aperture uncovered by a translating slide. +Shaft obstruction and edge leakage are not modeled. -The ECU uses ideal sensors and direct access to trapped fresh-air mass at -intake-valve closure to calculate each cylinder's fuel dose. This explicitly +The four-strokes use ideal sensors and direct access to trapped fresh-air mass +at intake-valve closure to calculate each cylinder's fuel dose. This explicitly permitted idealization replaces injector flow/duty-cycle and fuel-film models. -The nominal AFR is 14.7, enriched progressively toward 12.8 above 70% throttle. +Their nominal AFR is 14.7, enriched progressively toward 12.8 above 70% throttle. A PI idle controller targets 1800 RPM with filtered RPM feedback, bounded idle authority and anti-windup. Overrun fuel cut engages above 3500 RPM on a closed pedal and restores fuel below 2400 RPM to catch idle. The twin's rev limiter cuts fuel -and spark at 10000 RPM, recovering below 9750 RPM. The two-stroke has no rev -limiter; its idle, overrun and sensor-fault controls remain active. Invalid sensor readings -inhibit fuel. A simple RPM-based ignition advance map is provided. An already -burning charge completes combustion when a cut begins. +and spark at 10000 RPM, recovering below 9750 RPM. + +The two-stroke uses a simulated 38.1 mm slide carburetor instead of the EFI +dose. Effective pilot, needle and main restrictions draw liquid fuel from a +finite float bowl with the orifice relation `m_dot = Cd*A*sqrt(2*rho*delta_p)`. +The jet pressure signal comes from simulated mass flow, local venturi dynamic +pressure and the measured static pressure drop. A wall film deposits and +evaporates fuel during transients. Delivered fuel is added downstream of the +carb as a conserved species and travels through the plenum, reed, crankcase, +transfer duct and cylinder. Combustion is therefore limited by the fuel and +oxygen actually trapped, rather than receiving a cylinder-side AFR dose. + +The two-stroke has no rev limiter and deliberately does not use electronic +overrun fuel cut. On a high-RPM throttle lift, wet-wall and pilot fuel can pass +through unlit cycles if the trapped mixture cannot ignite. Spark remains at +the CDI timing; there is no scheduled overrun firing pattern. Unburned fuel +is transported out through the exhaust without pipe combustion. No sound +samples are mixed in. All ECU parameters are accessible calibrations, not measured engine maps. There is no starter motor, automatic restart, knock sensor, lambda sensor, @@ -255,11 +332,14 @@ catalyst, cold-start temperature model or injector electrical model. Idle control can still stall under excessive load and does not impose crankshaft speed or add artificial sustaining torque. Existing combustion/thermal models remain approximations; ECU changes do not establish their physical accuracy. -The two-stroke doses fuel after the last scavenging port closes, targets -3000 RPM idle, and uses idle ignition trim in addition to throttle control. Its +The two-stroke targets 3000 RPM idle and uses idle ignition trim in addition to +throttle control. Its reed is a damped lumped cantilever model with finite opening time and seating contacts. Cylinder and crankcase moving-volume work use adiabatic energy changes. -Reed material, combustion and controller parameters still require calibration. +Jet dimensions are effective flow calibrations rather than measured Mikuni +parts. Fuel uses the gas model after wall-film release; droplet size, +atomization, evaporation heat and oil are not modeled. Reed material, +combustion and controller parameters still require calibration. C# trace parity is intentionally changed by these control and physics changes. VS Code's default build task builds the desktop through CMake. Use diff --git a/include/FluidSim/Atmosphere.hpp b/include/FluidSim/Atmosphere.hpp index aa51d6c..036178b 100644 --- a/include/FluidSim/Atmosphere.hpp +++ b/include/FluidSim/Atmosphere.hpp @@ -31,6 +31,7 @@ private: port.temperature = temperature; port.enthalpy = specific_enthalpy(); port.air_fraction = 1.0f; + port.fuel_fraction = 0.0f; } }; diff --git a/include/FluidSim/AudioOutputStream.hpp b/include/FluidSim/AudioOutputStream.hpp index e15ec66..8b92f62 100644 --- a/include/FluidSim/AudioOutputStream.hpp +++ b/include/FluidSim/AudioOutputStream.hpp @@ -14,6 +14,8 @@ public: void set_speed(double speed) { speed_.store(speed, std::memory_order_relaxed); } std::uint64_t underrun_blocks() const { return underrun_blocks_.load(std::memory_order_relaxed); } double speed() const { return speed_.load(std::memory_order_relaxed); } + void set_auto_gain(bool enabled) { auto_gain_enabled_.store(enabled, std::memory_order_relaxed); } + bool auto_gain() const { return auto_gain_enabled_.load(std::memory_order_relaxed); } private: static constexpr std::size_t block_size = 512; @@ -22,6 +24,9 @@ private: std::array output_{}; std::atomic underrun_blocks_{}; std::atomic speed_{1.0}; + std::atomic auto_gain_enabled_{true}; + float rms_power_{}, peak_envelope_{}, automatic_gain_{1.f}; + float level_alpha_{}, peak_decay_{}, gain_rise_alpha_{}, gain_fall_alpha_{}; float previous_input_left_{}, previous_output_left_{}; float previous_input_right_{}, previous_output_right_{}; diff --git a/include/FluidSim/Carburetor.hpp b/include/FluidSim/Carburetor.hpp new file mode 100644 index 0000000..2865f5d --- /dev/null +++ b/include/FluidSim/Carburetor.hpp @@ -0,0 +1,46 @@ +#pragma once + +#include +#include + +namespace fluidsim::port { + +// Lumped slide-carburetor model for the crankcase-scavenged two-stroke. The +// jets use liquid-orifice flow from a float bowl into the measured venturi +// depression. Delivered fuel is injected into the downstream fluid volume and +// is then transported as a conserved species; this class never doses a cylinder. +class Carburetor { +public: + float fuel_density{740.0f}; + float pilot_jet_diameter{0.00040f}, main_jet_diameter{0.00150f}; + float needle_jet_diameter{0.00190f}; + float pilot_discharge_coefficient{0.62f}, main_discharge_coefficient{0.72f}; + float bowl_head_pressure{}; + float film_fraction{0.28f}, film_time_constant{0.11f}; + float bowl_capacity{0.008f}, bowl_refill_rate{0.0025f}; // kg, kg/s + float overrun_rpm{4200.0f}, overrun_exit_rpm{3300.0f}; + + explicit Carburetor(float throat_diameter = 0.0381f); + + // air_mass_flow is positive from atmosphere toward the engine. + void step(float dt, float pedal, float slide, float rpm, float air_mass_flow, + float air_density, float upstream_pressure, float downstream_pressure); + + float indicated_afr() const { return indicated_afr_; } + float fuel_flow() const { return delivered_fuel_flow_; } + float film_mass() const { return film_mass_; } + float bowl_level() const { return bowl_mass_ / bowl_capacity; } + bool overrun_active() const { return overrun_; } + bool fault() const { return fault_; } + float venturi_depression() const { return venturi_depression_; } + +private: + float bowl_mass_{}, film_mass_{}; + float throat_area_{}, venturi_depression_{}; + float filtered_air_flow_{}, filtered_fuel_flow_{}, delivered_fuel_flow_{}; + float indicated_afr_{13.5f}; + float previous_pedal_{}, lift_enrichment_{}; + bool initialized_{}, overrun_{}, fault_{}; +}; + +} // namespace fluidsim::port diff --git a/include/FluidSim/Component.hpp b/include/FluidSim/Component.hpp index d99376d..a25243f 100644 --- a/include/FluidSim/Component.hpp +++ b/include/FluidSim/Component.hpp @@ -12,9 +12,11 @@ struct Port { float temperature{300.0f}; float gamma_value{1.4f}, gas_constant_value{287.0f}; float air_fraction{1.0f}; + float fuel_fraction{}; float* owner_energy{}; float* owner_air_mass{}; float* owner_exhaust_mass{}; + float* owner_fuel_mass{}; }; // Small equivalent of the C# IComponent lifecycle. Components update their diff --git a/include/FluidSim/Crankcase.hpp b/include/FluidSim/Crankcase.hpp index f285b22..a92afe9 100644 --- a/include/FluidSim/Crankcase.hpp +++ b/include/FluidSim/Crankcase.hpp @@ -20,6 +20,13 @@ public: float geometry_work() const { return geometry_work_; } // J delivered by gas during the latest geometry step. float mass() const { return mass_; } float air_fraction() const { return air_fraction_; } + float fuel_fraction() const { return fuel_fraction_; } + float fuel_mass() const { return mass_ * fuel_fraction_; } + void initialize_mixture(float afr) { + fuel_fraction_ = 1.0f / (1.0f + afr); + air_fraction_ = 1.0f - fuel_fraction_; + sync_ports(); + } std::vector ports() override { return {&intake_port, &transfer_port}; } void pre_step(float dt); @@ -30,7 +37,7 @@ private: float crank_radius_{}, connecting_rod_length_{}, piston_area_{}; float clearance_volume_{}, obliquity_{}; float geometry_work_{}; - float mass_{}, energy_{}, air_fraction_{1.0f}; + float mass_{}, energy_{}, air_fraction_{1.0f}, fuel_fraction_{}; float pressure_{101325.0f}, temperature_{300.0f}, volume_{}, previous_volume_{}; void sync_ports(); diff --git a/include/FluidSim/DesktopScenario.hpp b/include/FluidSim/DesktopScenario.hpp index 3e9fd4c..0c3e835 100644 --- a/include/FluidSim/DesktopScenario.hpp +++ b/include/FluidSim/DesktopScenario.hpp @@ -14,7 +14,8 @@ public: two_ = std::make_unique(44100, config); else twin_ = std::make_unique(44100, kind == ScenarioKind::SaabB205 ? - FourStrokePreset::SaabB205 : FourStrokePreset::Twin, kind == ScenarioKind::TurboTwin || (kind == ScenarioKind::SaabB205 && saab_turbo)); + FourStrokePreset::SaabB205 : kind == ScenarioKind::AmericanV8 ? FourStrokePreset::AmericanV8 : FourStrokePreset::Twin, + kind == ScenarioKind::TurboTwin || (kind == ScenarioKind::SaabB205 && saab_turbo)); } AudioFrame process() { if (twin_) @@ -60,7 +61,7 @@ public: s.exhaust_pressure = two_->pipes[3].state(0).pressure; s.throttle = source.throttle; s.pedal = source.pedal; - s.target_afr = two_->ecu.target_afr(); + s.target_afr = two_->carburetor.indicated_afr(); s.valve1_intake_area = source.transfer_area; s.valve1_exhaust_area = source.exhaust_area; s.reed_open = source.reed_open; diff --git a/include/FluidSim/Ecu.hpp b/include/FluidSim/Ecu.hpp index 61b3f1e..04661a1 100644 --- a/include/FluidSim/Ecu.hpp +++ b/include/FluidSim/Ecu.hpp @@ -12,9 +12,12 @@ public: float idle_base_opening{0.08f}, idle_max_opening{0.25f}; float idle_kp{0.0001f}, idle_ki{0.00004f}; float throttle_time_constant{0.08f}, throttle_rate{2.0f}; // seconds, full travel/second + // Zero uses the opening calibration, preserving symmetric response. + float throttle_closing_time_constant{}, throttle_closing_rate{}; float rev_limit_rpm{10000.0f}, rev_resume_rpm{9750.0f}; bool rev_limit_enabled{true}; float overrun_cut_rpm{3500.0f}, overrun_resume_rpm{2400.0f}; + bool overrun_fuel_cut_enabled{true}; float idle_spark_gain{0.0f}, idle_minimum_advance{-5.0f}; // degrees/RPM, degrees BTDC float stoichiometric_afr{14.7f}, power_afr{12.8f}; bool boost_control_enabled{false}; @@ -30,7 +33,7 @@ public: float throttle_target() const { return target_; } float target_afr() const { return afr_; } float spark_advance() const { return spark_; } - bool fuel_cut() const { return limiter_ || overrun_ || fault_ || boost_cut_; } + bool fuel_cut() const { return limiter_ || (overrun_ && overrun_fuel_cut_enabled) || fault_ || boost_cut_; } bool rev_limited() const { return limiter_; } bool idle_active() const { return idle_; } bool sensor_fault() const { return fault_; } diff --git a/include/FluidSim/Engine.hpp b/include/FluidSim/Engine.hpp index 8343057..2db37e1 100644 --- a/include/FluidSim/Engine.hpp +++ b/include/FluidSim/Engine.hpp @@ -31,6 +31,7 @@ private: class EngineCylinder : public Component { public: Ecu* ecu{}; // Optional ideal electronic fuel/ignition controller. + bool transported_fuel{}; // Upstream carb meters the fuel species. float bore{0.075f}, stroke{0.068f}, rod{0.112f}, compression{11.5f}; float intake_valve_diameter{0.030f}, exhaust_valve_diameter{0.028f}; int intake_valve_count{1}, exhaust_valve_count{1}; @@ -60,6 +61,7 @@ public: protected: float energy_factor_{1.0f}; + float active_fuel_mass_{}; float rev_limit_timer_{}; std::mt19937 random_{0xF17D51u}; std::uniform_real_distribution unit_{0.0f, 1.0f}; diff --git a/include/FluidSim/Fluid.hpp b/include/FluidSim/Fluid.hpp index 2064d01..ad48b95 100644 --- a/include/FluidSim/Fluid.hpp +++ b/include/FluidSim/Fluid.hpp @@ -37,6 +37,7 @@ public: p.owner_energy = &energy_; p.owner_air_mass = &air_mass_; p.owner_exhaust_mass = &exhaust_mass_; + p.owner_fuel_mass = &fuel_mass_; sync_port(p); return p; } @@ -44,19 +45,30 @@ public: float density() const { return mass() / std::max(volume_, 1e-12f); } float temperature() const { return pressure() / std::max(density() * gas_constant, 1e-9f); } float volume() const { return volume_; } - float mass() const { return air_mass_ + exhaust_mass_; } + float mass() const { return air_mass_ + exhaust_mass_ + fuel_mass_; } float internal_energy() const { return energy_; } float air_fraction() const { return air_mass_ / std::max(mass(), 1e-12f); } + float fuel_fraction() const { return fuel_mass_ / std::max(mass(), 1e-12f); } + float fuel_mass() const { return fuel_mass_; } + void add_fuel_mass(float value, float temperature = 300.0f) { + if (!std::isfinite(value) || value < 0.0f || !std::isfinite(temperature) || temperature <= 0.0f) + throw std::invalid_argument("invalid fuel addition"); + fuel_mass_ += value; + energy_ += value * gas_constant / (gamma - 1.0f) * temperature; + sync_ports(); + } void set_pressure(float value, float temp = -1.0f) { if (!std::isfinite(value) || value <= 0 || !std::isfinite(temp)) throw std::invalid_argument("invalid volume pressure/temperature"); const float target_t = temp > 0.0f ? temp : temperature(); const float total = value * volume_ / (gas_constant * target_t); const float fraction = air_fraction(); + const float fuel = fuel_fraction(); air_mass_ = total * fraction; - exhaust_mass_ = total * (1.0f - fraction); + fuel_mass_ = total * fuel; + exhaust_mass_ = total * std::max(1.0f - fraction - fuel, 0.0f); energy_ = value * volume_ / (gamma - 1.0f); - air_roundoff_ = exhaust_roundoff_ = energy_roundoff_ = 0; + air_roundoff_ = exhaust_roundoff_ = fuel_roundoff_ = energy_roundoff_ = 0; sync_ports(); } void update(float dt) { @@ -64,16 +76,20 @@ public: throw std::invalid_argument("invalid volume timestep"); if (dt == 0) return; - double d_air = 0, d_exhaust = 0, d_energy = 0; + double d_air = 0, d_exhaust = 0, d_fuel = 0, d_energy = 0; const float fraction = air_fraction(); + const float fuel = fuel_fraction(); for (const Port &p : ports_) { const float flow_fraction = p.mass_flow >= 0.0f ? p.air_fraction : fraction; + const float flow_fuel = p.mass_flow >= 0.0f ? p.fuel_fraction : fuel; d_air += double(p.mass_flow) * flow_fraction; - d_exhaust += double(p.mass_flow) * (1.0f - flow_fraction); + d_fuel += double(p.mass_flow) * flow_fuel; + d_exhaust += double(p.mass_flow) * std::max(1.0f - flow_fraction - flow_fuel, 0.0f); d_energy += double(p.mass_flow) * p.enthalpy; } accumulate(air_mass_, air_roundoff_, d_air * dt); accumulate(exhaust_mass_, exhaust_roundoff_, d_exhaust * dt); + accumulate(fuel_mass_, fuel_roundoff_, d_fuel * dt); accumulate(energy_, energy_roundoff_, d_energy * dt); const float new_volume = volume_ + volume_change_rate * dt; if (!(new_volume > 0) || !std::isfinite(new_volume)) @@ -93,10 +109,11 @@ public: if (!(mass() > 0) || !(energy_ > 0) || !std::isfinite(mass()) || !std::isfinite(energy_)) throw std::runtime_error( "nonphysical volume state; reduce timestep or check external port fluxes"); - if (air_mass_ < -1e-12f || exhaust_mass_ < -1e-12f) + if (air_mass_ < -1e-12f || exhaust_mass_ < -1e-12f || fuel_mass_ < -1e-12f) throw std::runtime_error("negative volume species mass"); air_mass_ = std::max(air_mass_, 0.0f); exhaust_mass_ = std::max(exhaust_mass_, 0.0f); + fuel_mass_ = std::max(fuel_mass_, 0.0f); sync_ports(); } void update_state(float dt) override { update(dt); } @@ -110,9 +127,9 @@ public: const std::deque &ports() const { return ports_; } private: - float volume_, air_mass_{}, exhaust_mass_{}, energy_; + float volume_, air_mass_{}, exhaust_mass_{}, fuel_mass_{}, energy_; // Preserve sub-ulp transfers without changing the float port/budget ABI. - double air_roundoff_{}, exhaust_roundoff_{}, energy_roundoff_{}; + double air_roundoff_{}, exhaust_roundoff_{}, fuel_roundoff_{}, energy_roundoff_{}; static void accumulate(float& value, double& roundoff, double increment) { const double precise = double(value)+roundoff+increment; value = float(precise); @@ -127,6 +144,7 @@ private: p.temperature = temperature(); p.enthalpy = gamma / (gamma - 1.0f) * p.pressure / std::max(p.density, 1e-12f); p.air_fraction = air_fraction(); + p.fuel_fraction = fuel_fraction(); } void sync_ports() { for (Port &p : ports_) @@ -138,16 +156,18 @@ struct GasCell { float rho{1.225f}, momentum{}; float energy{1.225f * (101325.0f / (pipe_gm1 * 1.225f))}; float air{1.0f}; + float fuel{}; }; struct GasState { float rho{1.225f}, velocity{}, pressure{101325.0f}, air{1.0f}; + float fuel{}; }; class Pipe1D { public: - // Integrated face fluxes: kg/s, N, W, kg/s of fresh air, positive along +x. + // Integrated face fluxes: kg/s, N, W, kg/s of fresh air and fuel, positive along +x. struct Flux { - float m{}, p{}, e{}, y{}; + float m{}, p{}, e{}, y{}, fuel{}; }; float area{}, dx{}; float damping_multiplier{0.0f}, energy_relaxation_rate{0.0f}, ambient_pressure{101325.0f}; @@ -221,7 +241,9 @@ public: } static GasState primitive(const GasCell &c) { const float u = c.momentum / std::max(c.rho, 1e-12f); - return {c.rho, u, pipe_gm1 * (c.energy - .5f * c.momentum * u), std::clamp(c.air, 0.0f, 1.0f)}; + const float air = std::clamp(c.air, 0.0f, 1.0f); + return {c.rho, u, pipe_gm1 * (c.energy - .5f * c.momentum * u), air, + std::clamp(c.fuel, 0.0f, 1.0f - air)}; } static float wave(const GasState &s) { return std::abs(s.velocity) + std::sqrt(gamma * s.pressure / s.rho); @@ -253,7 +275,9 @@ public: float mom = cells[i].momentum - factor * (fr.p - fl.p - states_[i].pressure * (ar - al)); float en = cells[i].energy - factor * (fr.e - fl.e); const float air = cells[i].rho * cells[i].air - factor * (fr.y - fl.y); - if (!std::isfinite(rho) || !std::isfinite(mom) || !std::isfinite(en) || !std::isfinite(air)) + const float fuel = cells[i].rho * cells[i].fuel - factor * (fr.fuel - fl.fuel); + if (!std::isfinite(rho) || !std::isfinite(mom) || !std::isfinite(en) || + !std::isfinite(air) || !std::isfinite(fuel)) throw std::runtime_error("nonfinite finite-volume update"); if (rho < 1e-12f) { rho = 1e-12f; @@ -273,14 +297,18 @@ public: // than 1.18e-38 of the gas mass. Avoid repeatedly transporting // subnormal concentrations; total gas mass and energy are unchanged. const float fraction = std::clamp(air / rho, 0.0f, 1.0f); + const float fuel_fraction = std::clamp(fuel / rho, 0.0f, 1.0f - fraction); next_cells_[i] = {rho, mom, en, - fraction < std::numeric_limits::min() ? 0.0f : fraction}; + fraction < std::numeric_limits::min() ? 0.0f : fraction, + fuel_fraction < std::numeric_limits::min() ? 0.0f : fuel_fraction}; } cells.swap(next_cells_); prepared_ = false; left_override_ = right_override_ = false; } - static Flux scaled(Flux f, float a) { return {f.m * a, f.p * a, f.e * a, f.y * a}; } + static Flux scaled(Flux f, float a) { + return {f.m * a, f.p * a, f.e * a, f.y * a, f.fuel * a}; + } private: std::vector next_cells_; @@ -301,7 +329,8 @@ private: return {.5f * (ma + mb - c * (b.rho - a.rho)), .5f * (ma * a.velocity + a.pressure + mb * b.velocity + b.pressure - c * (mb - ma)), .5f * ((ea + a.pressure) * a.velocity + (eb + b.pressure) * b.velocity - c * (eb - ea)), - .5f * (ma * a.air + mb * b.air - c * (b.rho * b.air - a.rho * a.air))}; + .5f * (ma * a.air + mb * b.air - c * (b.rho * b.air - a.rho * a.air)), + .5f * (ma * a.fuel + mb * b.fuel - c * (b.rho * b.fuel - a.rho * a.fuel))}; } }; @@ -329,7 +358,8 @@ public: [&](const Budget &b) { return b.air == port.owner_air_mass; }); if (it == budgets_.end()) { o.budget = budgets_.size(); - budgets_.push_back({port.owner_air_mass, port.owner_exhaust_mass, port.owner_energy}); + budgets_.push_back({port.owner_air_mass, port.owner_exhaust_mass, + port.owner_fuel_mass, port.owner_energy}); } else o.budget = static_cast(it - budgets_.begin()); } @@ -385,7 +415,7 @@ public: if (!std::isfinite(dt) || dt <= 0) throw std::invalid_argument("boundary timestep must be positive"); for (auto &b : budgets_) { - b.mass = .25f * (*b.air + (b.exhaust ? *b.exhaust : 0)); + b.mass = .25f * (*b.air + (b.exhaust ? *b.exhaust : 0) + (b.fuel ? *b.fuel : 0)); b.heat = b.energy ? .25f * (*b.energy) : std::numeric_limits::max(); } for (auto &o : orifices_) @@ -398,7 +428,7 @@ public: private: struct Budget { - float *air{}, *exhaust{}, *energy{}; + float *air{}, *exhaust{}, *fuel{}, *energy{}; float mass{}, heat{}; }; struct Orifice { @@ -587,15 +617,18 @@ private: o.length > 0 ? .5f * (v.pressure + s.pressure) : static_cast(duct_pressure); const float direction = o.left ? 1.f : -1.f; const float species = from ? v.air_fraction : s.air; + const float fuel = from ? v.fuel_fraction : s.fuel; const float signed_mass = direction * mdot; p.set_boundary_flux(o.left, {signed_mass, mdot * velocity + face_pressure * area + s.pressure * (face_area - area), - signed_mass * h, signed_mass * species}); - ghost(p, o.left, {density, direction * velocity, face_pressure, species}); + signed_mass * h, signed_mass * species, signed_mass * fuel}); + ghost(p, o.left, {density, direction * velocity, face_pressure, species, fuel}); v.mass_flow = -mdot; v.enthalpy = h; if (!from) v.air_fraction = s.air; + if (!from) + v.fuel_fraction = s.fuel; } void resolve_open(Open &o) { auto &p = pipes_[o.pipe]; @@ -611,7 +644,9 @@ private: const float u = o.left ? s.velocity - 2 * c / gm1 + 2 * cb / gm1 : s.velocity + 2 * c / gm1 - 2 * cb / gm1; const bool out = o.left ? s.velocity <= -c : s.velocity >= c; - GasState g = out ? s : GasState{rho, u, o.pressure, (o.left ? u >= 0 : u <= 0) ? 1.f : s.air}; + const bool ambient_inflow = o.left ? u >= 0 : u <= 0; + GasState g = out ? s : GasState{rho, u, o.pressure, ambient_inflow ? 1.f : s.air, + ambient_inflow ? 0.f : s.fuel}; ghost(p, o.left, g); const float a = std::min(o.area, p.cell_area(i)); auto f = Pipe1D::scaled(o.left ? Pipe1D::riemann(g, s) : Pipe1D::riemann(s, g), a); @@ -631,9 +666,11 @@ private: const float area = std::min(a.cell_area(ia), b.cell_area(ib)); auto f = Pipe1D::scaled(Pipe1D::riemann(sa, sb), area); a.set_boundary_flux(o.left_a, - {da * f.m, f.p + sa.pressure * (a.cell_area(ia) - area), da * f.e, da * f.y}); + {da * f.m, f.p + sa.pressure * (a.cell_area(ia) - area), da * f.e, + da * f.y, da * f.fuel}); b.set_boundary_flux(o.left_b, - {db * f.m, f.p + sb.pressure * (b.cell_area(ib) - area), db * f.e, db * f.y}); + {db * f.m, f.p + sb.pressure * (b.cell_area(ib) - area), db * f.e, + db * f.y, db * f.fuel}); } }; diff --git a/include/FluidSim/OutdoorAudio.hpp b/include/FluidSim/OutdoorAudio.hpp index 1bcd3aa..1cf42bb 100644 --- a/include/FluidSim/OutdoorAudio.hpp +++ b/include/FluidSim/OutdoorAudio.hpp @@ -14,6 +14,8 @@ struct OutdoorAudioConfig { Vec3 listener_right{1.0f, 0.0f, 0.0f}; Vec3 intake_position{-0.2f, 0.7f, 0.0f}; Vec3 exhaust_position{0.55f, 0.5f, 0.0f}; + Vec3 second_exhaust_position{-0.55f, 0.5f, 0.0f}; + bool dual_exhaust{}; float ear_spacing{0.18f}, sound_speed{343.0f}; float ground_y{0.0f}, wall_x{-1.5f}; float ground_reflection{0.35f}, wall_reflection{0.65f}; @@ -26,7 +28,7 @@ struct OutdoorAudioConfig { class OutdoorAudio { public: explicit OutdoorAudio(int sample_rate = 44100, OutdoorAudioConfig config = {}); - AudioFrame process(float intake_mass_flow, float exhaust_mass_flow); + AudioFrame process(float intake_mass_flow, float exhaust_mass_flow, float second_exhaust_mass_flow = 0); void set_scene(OutdoorScene scene) { scene_ = scene; } OutdoorScene scene() const { return scene_; } void set_monitor_filter(bool enabled) { monitor_filter_ = enabled; } // Set before playback/capture. @@ -60,7 +62,7 @@ private: OutdoorAudioConfig config_; OutdoorScene scene_{OutdoorScene::Wall}; bool monitor_filter_{true}; - std::array emitters_; + std::array emitters_; std::array tail_; std::array weights_{1, 0, 0, 0}; float scene_alpha_{}, tail_alpha_{}, tail_mix_{}; diff --git a/include/FluidSim/Scenario.hpp b/include/FluidSim/Scenario.hpp index 8929162..939002a 100644 --- a/include/FluidSim/Scenario.hpp +++ b/include/FluidSim/Scenario.hpp @@ -19,8 +19,8 @@ namespace fluidsim::port { -enum class ScenarioKind { Twin, TwoStroke, SaabB205, TurboTwin }; -enum class FourStrokePreset { Twin, SaabB205 }; +enum class ScenarioKind { Twin, TwoStroke, SaabB205, TurboTwin, AmericanV8 }; +enum class FourStrokePreset { Twin, SaabB205, AmericanV8 }; struct ScenarioSnapshot { bool turbo_enabled{}; @@ -28,7 +28,7 @@ struct ScenarioSnapshot { TurboSnapshot turbo{}; int cylinder_count{2}; float rev_limit_rpm{10000}; - std::array cylinder_pressures{}, piston_fractions{}, intake_areas{}, exhaust_areas{}; + std::array cylinder_pressures{}, piston_fractions{}, intake_areas{}, exhaust_areas{}; std::uint64_t solver_substeps{}, boundary_iterations{}; float rpm{}, speed_kmh{}, kinetic_energy{}, average_power{}, average_torque{}; float cylinder_pressure1{}, cylinder_pressure2{}; @@ -43,8 +43,8 @@ struct ScenarioSnapshot { OutdoorScene audio_scene{OutdoorScene::Wall}; int gear{}; bool rev_limit1{}, rev_limit2{}, mechanical_sound_enabled{}; - std::array, 12> pipe_states; - std::array, 12> pipe_areas, pipe_lengths; + std::array, 22> pipe_states; + std::array, 22> pipe_areas, pipe_lengths; void capture_pipes(const std::vector& pipes) { for (std::size_t k = 0; k < pipes.size(); ++k) { pipe_areas[k] = pipes[k].cell_areas; @@ -57,14 +57,15 @@ struct ScenarioSnapshot { }; // Shared four-stroke topology: independent runners/headers and shared plenums. -// Presets retain the original 600 cc twin and add a B205-sized inline-four. +// Presets include a 600 cc twin, B205-sized inline-four and American 5.0 V8. class FourStrokeScenario : public MotorcycleTest { public: explicit FourStrokeScenario(int sample_rate = 44100, FourStrokePreset preset = FourStrokePreset::Twin, bool turbo_enabled = false); const FourStrokePreset preset; - int cylinder_count() const { return preset == FourStrokePreset::SaabB205 ? 4 : 2; } - std::array cylinders() { return {&cylinder1, &cylinder2, &cylinder3, &cylinder4}; } + int cylinder_count() const { return preset == FourStrokePreset::AmericanV8 ? 8 : preset == FourStrokePreset::SaabB205 ? 4 : 2; } + std::array cylinders() { return {&cylinder1, &cylinder2, &cylinder3, &cylinder4, + &cylinder5, &cylinder6, &cylinder7, &cylinder8}; } AudioFrame process(); void set_controls(float throttle_value, float load_value, float clutch_value); void toggle_mechanical_sound(); @@ -85,9 +86,11 @@ public: std::unique_ptr turbo; Crankshaft crankshaft{}; FourStrokeCylinder cylinder1{}, cylinder2{}, cylinder3{}, cylinder4{}; + FourStrokeCylinder cylinder5{}, cylinder6{}, cylinder7{}, cylinder8{}; Atmosphere intake_atmosphere{}; Volume0D intake_plenum{1.0e-3f, 101325.0f, 300.0f}; Volume0D exhaust_collector{0.15e-3f, 101325.0f, 700.0f}; + Volume0D second_bank_collector{1.5e-3f, 101325.0f, 700.0f}; std::vector pipes; float intake_flow{}, exhaust_flow{}; float last_mix_l{}, last_mix_r{}; @@ -95,7 +98,7 @@ public: private: float intake_boundary_area_{}, throttle_area_{}, runner_area_{}, header_area_{}, tail_area_{}; - std::array intake_valve_areas_{}, exhaust_valve_areas_{}; + std::array intake_valve_areas_{}, exhaust_valve_areas_{}; float dt_{}; BoundarySystem boundaries_; Solver solver_; @@ -108,6 +111,11 @@ private: }; using FourStrokeTwinScenario = FourStrokeScenario; +class AmericanV8Scenario : public FourStrokeScenario { +public: + explicit AmericanV8Scenario(int sample_rate = 44100) + : FourStrokeScenario(sample_rate, FourStrokePreset::AmericanV8) {} +}; class TurboTwinScenario : public FourStrokeScenario { public: explicit TurboTwinScenario(int sample_rate = 44100) : FourStrokeScenario(sample_rate, FourStrokePreset::Twin, true) {} diff --git a/include/FluidSim/SoundEngine.hpp b/include/FluidSim/SoundEngine.hpp index abf76a3..d72b7a4 100644 --- a/include/FluidSim/SoundEngine.hpp +++ b/include/FluidSim/SoundEngine.hpp @@ -21,6 +21,8 @@ public: double speed() const; void set_volume(float volume); float volume() const; + void set_auto_gain(bool enabled) { stream_.set_auto_gain(enabled); } + bool auto_gain() const { return stream_.auto_gain(); } std::uint64_t underrun_blocks() const { return stream_.underrun_blocks(); } bool is_playing() const { return is_playing_; } diff --git a/include/FluidSim/TwoStrokeCylinder.hpp b/include/FluidSim/TwoStrokeCylinder.hpp index 1128715..7a9fe0f 100644 --- a/include/FluidSim/TwoStrokeCylinder.hpp +++ b/include/FluidSim/TwoStrokeCylinder.hpp @@ -5,6 +5,8 @@ namespace fluidsim::port { +class Carburetor; + // Two-stroke variant of the reciprocating cylinder model. The shared four- // stroke implementation supplies geometry, state balance, and combustion // parameters; cycle events and pumping torque are specialized here. @@ -19,6 +21,7 @@ public: float exhaust_open() const { return 180.0f - exhaust_duration_ * 0.5f; } float exhaust_close() const { return 180.0f + exhaust_duration_ * 0.5f; } void set_crankcase(Crankcase& crankcase) { crankcase_ = &crankcase; } + void set_carburetor(Carburetor& carburetor) { carburetor_ = &carburetor; } float transfer_width{.050f}, exhaust_width{.036f}; float minimum_burn_seconds{.0008f}; std::uint64_t combustion_count() const { return combustion_count_; } @@ -29,6 +32,7 @@ public: private: float transfer_duration_{}, exhaust_duration_{}; Crankcase* crankcase_{}; + Carburetor* carburetor_{}; float active_burn_duration_{60}; float current_rpm_{}; diff --git a/include/FluidSim/TwoStrokeScenario.hpp b/include/FluidSim/TwoStrokeScenario.hpp index abc7341..677f55e 100644 --- a/include/FluidSim/TwoStrokeScenario.hpp +++ b/include/FluidSim/TwoStrokeScenario.hpp @@ -2,6 +2,7 @@ #include "FluidSim/AudioProcessing.hpp" #include "FluidSim/Crankcase.hpp" +#include "FluidSim/Carburetor.hpp" #include "FluidSim/Fluid.hpp" #include "FluidSim/OutdoorAudio.hpp" #include "FluidSim/Atmosphere.hpp" @@ -44,6 +45,7 @@ public: float throttle{}; const TwoStrokeConfig configuration; Ecu ecu{}; + Carburetor carburetor; ReedValve reed{}; Atmosphere intake_atmosphere{}; float intake_flow{}, exhaust_flow{}; @@ -65,6 +67,7 @@ private: OutdoorAudio audio_; static float area_for_diameter(float diameter); + static float slide_area(float diameter, float opening); static std::vector create_pipes(const TwoStrokeConfig& config); }; diff --git a/src/audio_output_stream.cpp b/src/audio_output_stream.cpp index 75a7359..52e85e3 100644 --- a/src/audio_output_stream.cpp +++ b/src/audio_output_stream.cpp @@ -7,6 +7,11 @@ namespace fluidsim::port { AudioOutputStream::AudioOutputStream(SimulationRingBuffer& ring, unsigned int sample_rate) : ring_(ring) { + const float rate = static_cast(sample_rate); + level_alpha_ = -std::expm1(-1.f / (rate * .15f)); + peak_decay_ = std::exp(-1.f / (rate * .15f)); + gain_rise_alpha_ = -std::expm1(-1.f / (rate * 1.5f)); + gain_fall_alpha_ = -std::expm1(-1.f / (rate * .03f)); initialize(2, sample_rate, {sf::SoundChannel::FrontLeft, sf::SoundChannel::FrontRight}); } @@ -14,18 +19,38 @@ bool AudioOutputStream::onGetData(Chunk& data) { const std::size_t read = ring_.read_interpolated(temp_left_, temp_right_, speed()); if (read < block_size) underrun_blocks_.fetch_add(1, std::memory_order_relaxed); constexpr float dc_coefficient = 0.9995f; + constexpr float playback_gain = 4.0f; // +12 dB desktop playback boost + const bool automatic = auto_gain(); // Keep callback work bounded: no allocation, waits or simulation here. for (std::size_t i = 0; i < block_size; ++i) { float left{}, right{}; if (i < read) { - const float input_left = std::isfinite(temp_left_[i]) ? temp_left_[i] : 0.0f; - const float input_right = std::isfinite(temp_right_[i]) ? temp_right_[i] : 0.0f; + const float input_left = std::isfinite(temp_left_[i]) ? temp_left_[i] * playback_gain : 0.0f; + const float input_right = std::isfinite(temp_right_[i]) ? temp_right_[i] * playback_gain : 0.0f; left = input_left - previous_input_left_ + dc_coefficient * previous_output_left_; right = input_right - previous_input_right_ + dc_coefficient * previous_output_right_; previous_input_left_ = input_left; previous_input_right_ = input_right; previous_output_left_ = left; previous_output_right_ = right; + // Stereo-linked level detection after DC removal. Freeze on + // underrun; silence returns to unity rather than chasing noise. + const float peak = std::max(std::abs(left), std::abs(right)); + rms_power_ += (.5f * (left*left + right*right) - rms_power_) * level_alpha_; + peak_envelope_ = std::max(peak, peak_envelope_ * peak_decay_); + float target = 1.f; + if (automatic && rms_power_ > 1e-8f) { + target = std::clamp(.08f / std::sqrt(rms_power_), .125f, 12.f); + target = std::min(target, .85f / std::max(peak_envelope_, 1e-6f)); + } + const float alpha = target < automatic_gain_ ? gain_fall_alpha_ : gain_rise_alpha_; + automatic_gain_ += (target - automatic_gain_) * alpha; + // Catch a sudden rev while the slow idle boost is still elevated. + // One linked multiplier preserves the stereo image. + const float gain = (automatic || automatic_gain_ > 1.001f) ? + std::min(automatic_gain_, .98f / std::max(peak, 1e-6f)) : automatic_gain_; + left *= gain; + right *= gain; if (starving_) { transition_left_ = last_left_ - left; transition_right_ = last_right_ - right; diff --git a/src/carburetor.cpp b/src/carburetor.cpp new file mode 100644 index 0000000..87848b6 --- /dev/null +++ b/src/carburetor.cpp @@ -0,0 +1,102 @@ +#include "FluidSim/Carburetor.hpp" + +#include +#include + +namespace fluidsim::port { +namespace { +constexpr float pi = 3.14159265358979323846f; +float circle_area(float diameter) { return .25f * pi * diameter * diameter; } +float smoothstep(float value) { + value = std::clamp(value, 0.0f, 1.0f); + return value * value * (3.0f - 2.0f * value); +} +} + +Carburetor::Carburetor(float throat_diameter) + : bowl_mass_(bowl_capacity), throat_area_(circle_area(std::max(throat_diameter, 1e-4f))) {} + +void Carburetor::step(float dt, float pedal, float slide, float rpm, float air_mass_flow, + float air_density, float upstream_pressure, float downstream_pressure) { + fault_ = !std::isfinite(dt) || dt <= 0.0f || !std::isfinite(pedal) || + !std::isfinite(slide) || !std::isfinite(rpm) || rpm < 0.0f || + !std::isfinite(air_mass_flow) || !std::isfinite(air_density) || air_density <= 0.0f || + !std::isfinite(upstream_pressure) || + !std::isfinite(downstream_pressure) || upstream_pressure <= 0.0f || + downstream_pressure <= 0.0f; + if (fault_) { + delivered_fuel_flow_ = 0.0f; + overrun_ = false; + return; + } + + pedal = std::clamp(pedal, 0.0f, 1.0f); + slide = std::clamp(slide, 0.0f, 1.0f); + air_mass_flow = std::max(air_mass_flow, 0.0f); + if (!initialized_) { + previous_pedal_ = pedal; + initialized_ = true; + } + const bool throttle_lift = previous_pedal_ > 0.08f && pedal <= 0.005f && rpm > overrun_rpm; + if (throttle_lift) + lift_enrichment_ = 1.0f; + previous_pedal_ = pedal; + lift_enrichment_ *= std::exp(-dt / 0.32f); + if (pedal > 0.005f || rpm < overrun_exit_rpm) + overrun_ = false; + else if (throttle_lift) + overrun_ = true; + + const float air_alpha = -std::expm1(-dt / 0.012f); + filtered_air_flow_ += (air_mass_flow - filtered_air_flow_) * air_alpha; + + // The open slide controls the effective venturi area. Its local dynamic + // pressure and the measured static-pressure loss provide the jet signal. + const float venturi_area = throat_area_ * 0.55f * (0.10f + 0.90f * slide); + const float velocity = filtered_air_flow_ / std::max(air_density * venturi_area, 1e-9f); + const float dynamic_pressure = 0.5f * air_density * velocity * velocity; + const float static_drop = std::max(upstream_pressure - downstream_pressure, 0.0f); + venturi_depression_ = std::clamp(dynamic_pressure + 0.18f * static_drop, 0.0f, 60000.0f); + + const float pilot_area = circle_area(pilot_jet_diameter); + const float main_area = circle_area(main_jet_diameter); + const float needle_jet_area = circle_area(needle_jet_diameter); + const float needle_lift = smoothstep((slide - 0.04f) / 0.86f); + const float needle_area = needle_jet_area * needle_lift; + const float main_path_area = needle_area > 1e-12f ? + 1.0f / std::sqrt(1.0f / (main_area * main_area) + 1.0f / (needle_area * needle_area)) : 0.0f; + const float pilot_depression = bowl_head_pressure + (0.62f * static_drop + + 0.25f * dynamic_pressure) * (1.0f + 0.25f * lift_enrichment_); + const float main_depression = bowl_head_pressure + venturi_depression_; + const float pilot_flow = pilot_discharge_coefficient * pilot_area * + std::sqrt(2.0f * fuel_density * std::max(pilot_depression, 0.0f)); + const float main_flow = main_discharge_coefficient * main_path_area * + std::sqrt(2.0f * fuel_density * std::max(main_depression, 0.0f)); + float demanded_fuel_flow = pilot_flow + main_flow; + + const float deposit_fraction = std::clamp(film_fraction * (0.55f + 0.45f * slide), 0.0f, 0.8f); + const float deposited = demanded_fuel_flow * deposit_fraction; + const float evaporation_time = film_time_constant * (overrun_ ? 0.58f : 1.0f); + const float evaporated = std::min(film_mass_ / dt, + film_mass_ / std::max(evaporation_time, 1e-4f)); + float jet_flow = demanded_fuel_flow; + const float available_flow = bowl_mass_ / dt; + if (jet_flow > available_flow) { + const float scale = available_flow / std::max(jet_flow, 1e-12f); + jet_flow *= scale; + } + bowl_mass_ -= jet_flow * dt; + bowl_mass_ += std::min(bowl_refill_rate * dt, bowl_capacity - bowl_mass_); + bowl_mass_ = std::clamp(bowl_mass_, 0.0f, bowl_capacity); + + // Deposition cannot exceed the fuel that actually left the bowl. + const float actual_deposit = std::min(deposited, jet_flow); + film_mass_ = std::max(film_mass_ + (actual_deposit - evaporated) * dt, 0.0f); + delivered_fuel_flow_ = std::max(jet_flow - actual_deposit + evaporated, 0.0f); + const float fuel_alpha = -std::expm1(-dt / 0.025f); + filtered_fuel_flow_ += (delivered_fuel_flow_ - filtered_fuel_flow_) * fuel_alpha; + if (filtered_air_flow_ > 1e-7f && filtered_fuel_flow_ > 1e-9f) + indicated_afr_ = std::clamp(filtered_air_flow_ / filtered_fuel_flow_, 4.0f, 40.0f); +} + +} // namespace fluidsim::port diff --git a/src/crankcase.cpp b/src/crankcase.cpp index d65e82d..d70af4d 100644 --- a/src/crankcase.cpp +++ b/src/crankcase.cpp @@ -51,18 +51,25 @@ void Crankcase::pre_step(float /*dt*/) { } void Crankcase::update_state(float dt) { - float delta_mass = 0.0f, delta_air_mass = 0.0f, delta_energy = 0.0f; + float delta_mass = 0.0f, delta_air_mass = 0.0f, delta_fuel_mass = 0.0f, delta_energy = 0.0f; for (Port* port : {&intake_port, &transfer_port}) { const float flow = port->mass_flow; const float fraction = flow >= 0.0f ? port->air_fraction : air_fraction_; + const float fuel = flow >= 0.0f ? port->fuel_fraction : fuel_fraction_; delta_mass += flow * dt; delta_air_mass += flow * fraction * dt; + delta_fuel_mass += flow * fuel * dt; delta_energy += flow * port->enthalpy * dt; } const float old_air_mass = mass_ * air_fraction_; + const float old_fuel_mass = mass_ * fuel_fraction_; mass_ += delta_mass; energy_ += delta_energy; - if (mass_ > 1e-12f) air_fraction_ = std::clamp((old_air_mass + delta_air_mass) / mass_, 0.0f, 1.0f); + if (mass_ > 1e-12f) { + air_fraction_ = std::clamp((old_air_mass + delta_air_mass) / mass_, 0.0f, 1.0f); + fuel_fraction_ = std::clamp((old_fuel_mass + delta_fuel_mass) / mass_, + 0.0f, 1.0f - air_fraction_); + } if (!(mass_ > 0.0f) || !(energy_ > 0.0f) || !std::isfinite(mass_) || !std::isfinite(energy_)) throw std::runtime_error("nonphysical crankcase state"); @@ -80,6 +87,7 @@ void Crankcase::sync_ports() { port->temperature = temperature_; port->enthalpy = enthalpy; port->air_fraction = air_fraction_; + port->fuel_fraction = fuel_fraction_; } } diff --git a/src/desktop_main.cpp b/src/desktop_main.cpp index 4735804..a3ebd0d 100644 --- a/src/desktop_main.cpp +++ b/src/desktop_main.cpp @@ -57,14 +57,15 @@ void draw_label(sf::RenderWindow& window, const sf::Font* font, const std::strin } std::optional choose_scenario(sf::RenderWindow& window, const sf::Font* font) { - const std::array cards{{{{32, 110}, {656, 84}}, {{32, 206}, {656, 84}}, {{32, 302}, {656, 84}}, {{32, 398}, {656, 84}}}}; - const sf::FloatRect start_button{{32, 506}, {656, 52}}; - const std::array titles{"600 cc parallel twin", "YZ125 two-stroke", "Saab B205 / 9-5 wagon", "600 cc twin + VZ21-size turbo"}; - const std::array descriptions{ + const std::array cards{{{{32, 110}, {656, 72}}, {{32, 190}, {656, 72}}, {{32, 270}, {656, 72}}, {{32, 350}, {656, 72}}, {{32, 430}, {656, 72}}}}; + const sf::FloatRect start_button{{32, 522}, {656, 52}}; + const std::array titles{"600 cc parallel twin", "YZ125 two-stroke", "Saab B205 / 9-5 wagon", "600 cc twin + VZ21-size turbo", "1980s American 5.0 V8"}; + const std::array descriptions{ "Four-stroke twin, motorcycle test and dyno", "125 cc single, expansion chamber, motorcycle test and dyno", "2.0 litre inline-four, GT17 turbo, ECU boost control and five-speed wagon", - "Coupled turbine/compressor, shaft inertia and pressure diagnostics"}; + "Coupled turbine/compressor, shaft inertia and pressure diagnostics", + "EFI cross-plane V8, heavy flywheel, independent exhaust banks"}; int selected = 0; while (window.isOpen()) { while (const auto event = window.pollEvent()) { @@ -72,26 +73,27 @@ std::optional choose_scenario(sf::RenderWindow& window, const sf:: if (const auto* key = event->getIf()) { switch (key->code) { case sf::Keyboard::Key::Escape: return std::nullopt; - case sf::Keyboard::Key::Up: selected = (selected + 3) % 4; break; + case sf::Keyboard::Key::Up: selected = (selected + 4) % 5; break; case sf::Keyboard::Key::Down: - case sf::Keyboard::Key::Tab: selected = (selected + 1) % 4; break; + case sf::Keyboard::Key::Tab: selected = (selected + 1) % 5; break; case sf::Keyboard::Key::Num1: selected = 0; break; case sf::Keyboard::Key::Num2: selected = 1; break; case sf::Keyboard::Key::Num3: selected = 2; break; case sf::Keyboard::Key::Num4: selected = 3; break; + case sf::Keyboard::Key::Num5: selected = 4; break; case sf::Keyboard::Key::Enter: return static_cast(selected); default: break; } } if (const auto* move = event->getIf()) { const auto point = window.mapPixelToCoords(move->position); - for (int i = 0; i < 4; ++i) + for (int i = 0; i < 5; ++i) if (cards[i].contains(point)) selected = i; } if (const auto* click = event->getIf(); click && click->button == sf::Mouse::Button::Left) { const auto point = window.mapPixelToCoords(click->position); - for (int i = 0; i < 4; ++i) + for (int i = 0; i < 5; ++i) if (cards[i].contains(point)) return static_cast(i); if (start_button.contains(point)) return static_cast(selected); } @@ -99,7 +101,7 @@ std::optional choose_scenario(sf::RenderWindow& window, const sf:: window.clear({9, 12, 17}); draw_label(window, font, "FluidSim", {32, 24}, sf::Color::White, 28); draw_label(window, font, "Choose a scenario to start", {32, 66}, {180, 190, 205}, 16); - for (int i = 0; i < 4; ++i) { + for (int i = 0; i < 5; ++i) { sf::RectangleShape card(cards[i].size); card.setPosition(cards[i].position); card.setFillColor(i == selected ? sf::Color(24, 46, 60) : sf::Color(18, 24, 32)); @@ -115,7 +117,7 @@ std::optional choose_scenario(sf::RenderWindow& window, const sf:: start.setPosition(start_button.position); start.setFillColor({35, 105, 115}); window.draw(start); - draw_label(window, font, "Start selected scenario [Enter]", {52, 519}, sf::Color::White, 17); + draw_label(window, font, "Start selected scenario [Enter]", {52, 537}, sf::Color::White, 17); draw_label(window, font, "Click a scenario or use Up/Down. Esc: quit", {32, 574}, {150, 165, 180}, 12); window.display(); } @@ -296,13 +298,13 @@ void draw_dyno_curve(sf::RenderWindow& window, const sf::Font* font, const DynoC void draw_dashboard(sf::RenderWindow& window, const sf::Font* font, const ScenarioSnapshot& state, float throttle_target, float brake, double playback_speed, std::size_t buffer_fill, - std::size_t buffer_capacity, bool time_warp, std::uint64_t underrun_blocks) { + std::size_t buffer_capacity, bool time_warp, std::uint64_t underrun_blocks, bool auto_gain) { std::ostringstream status; status << (time_warp ? "[Space] normal speed" : "[Space] real time") << " Sim " << std::fixed << std::setprecision(3) << playback_speed << "x" << " Audio buffer " << buffer_fill << "/" << buffer_capacity << " Underruns " << underrun_blocks - << (state.cylinder_count == 4 ? "\nWagon " : "\nBike ") << std::setprecision(1) << state.speed_kmh << " km/h Gear " + << (state.cylinder_count == 8 ? "\nCoupe " : state.cylinder_count == 4 ? "\nWagon " : "\nBike ") << std::setprecision(1) << state.speed_kmh << " km/h Gear " << (state.gear == 0 ? "N" : std::to_string(state.gear)) << " Throttle " << state.throttle * 100.0f << "% / target " << throttle_target * 100.0f << "% [W, up/down]" << " Clutch " << state.clutch * 100.0f << "% [C] Engine load " << brake * 100.0f << "% [left/right]" @@ -310,6 +312,7 @@ void draw_dashboard(sf::RenderWindow& window, const sf::Font* font, const Scenar << "\nScene [R]: " << scene_name(state.audio_scene) << " " << (state.fuel_cut ? "FUEL CUT" : (state.idle_control ? "IDLE CONTROL" : "RUNNING")); if (!state.two_stroke) status << " Mechanical [M]: " << (state.mechanical_sound_enabled ? "ON" : "OFF"); + status << " Auto gain [G]: " << (auto_gain ? "ON" : "OFF"); if (font) { sf::RectangleShape panel({1260.0f, 92.0f}); panel.setPosition({10.0f, 10.0f}); @@ -350,7 +353,7 @@ void draw_two_stroke_scene(sf::RenderWindow& window,const sf::Font* font,const S draw_label(window,font,"CRANKCASE",{455,385},sf::Color::White,14); draw_pipe(window,s,0,460,30,170); draw_volume(window,s.intake_pressure,{{180,425},{95,70}}); - draw_label(window,font,"INTAKE",{180,400},sf::Color::White,14); + draw_label(window,font,"CARB / INTAKE",{180,400},sf::Color::White,14); draw_pipe(window,s,1,460,285,445); draw_label(window,font,s.reed_open?"REED OPEN":"REED CLOSED",{290,485},s.reed_open?sf::Color::Green:sf::Color::Red,12); draw_pipe(window,s,2,350,315,450); @@ -365,7 +368,8 @@ void draw_two_stroke_scene(sf::RenderWindow& window,const sf::Font* font,const S draw_label(window,font,s.racing_exhaust ? "HEADER / DIFFUSERS / BELLY / BAFFLE / STINGER" : "STRAIGHT EXHAUST",{670,205},sf::Color::White,13); std::ostringstream pressures; pressures<1e-7f, state.exhaust_areas[i]>1e-7f, x); + draw_pipe(window, state, 9+i, top+150, x-10, x+90); + draw_pipe(window, state, 1+i, top+125, x-10, x+90); + draw_graph_line(window, {x+90,top+150}, {830,top+150}, {200,120,80}); + } + draw_graph_line(window, {830,380}, {850,300}, {200,120,80}); + draw_graph_line(window, {830,605}, {850,375}, {200,120,80}); + draw_pipe(window, state, 17, 300, 850, 1250); + draw_pipe(window, state, 18, 375, 850, 1250); + draw_label(window, font, "SEPARATE BANK COLLECTORS / DUAL TAILS", {850,410}, sf::Color::White, 11); + draw_label(window, font, "EFI AFR " + std::to_string(state.target_afr), {850,430}, sf::Color::White, 13); + draw_label(window, font, "302-class | shared 4 L intake, 12 cm runners | 7:1 compression | 750 rpm idle", {25,690}, {180,190,205}, 12); + } + else if (state.two_stroke) draw_two_stroke_scene(window, font, state); else if (state.cylinder_count == 4) draw_b205_scene(window, font, state); else draw_twin_scene(window, font, state); draw_dashboard(window, font, state, throttle_target, brake, playback_speed, buffer_fill, - buffer_capacity, time_warp, underrun_blocks); + buffer_capacity, time_warp, underrun_blocks, auto_gain); draw_dyno_curve(window, font, dyno_curve, 860, 455, 400, 250, state.rpm, power_kw); window.display(); } @@ -486,8 +514,8 @@ int main(int argc, char** argv) { if((option=="--scenario"||option=="--scene"||option=="--exhaust"||option=="--turbo")&&i+1 0 ? + throttle_closing_time_constant : throttle_time_constant; + const float rate = closing && throttle_closing_rate > 0 ? throttle_closing_rate : throttle_rate; + const float change = (target_ - opening_) * (-std::expm1(-dt / std::max(response, 1e-4f))); + const float travel = std::max(rate, 0.0f) * dt; opening_ = std::clamp(opening_ + std::clamp(change, -travel, travel), 0.0f, 1.0f); const float enrichment = std::clamp((opening_ - 0.7f) / 0.3f, 0.0f, 1.0f); afr_ = stoichiometric_afr + (power_afr - stoichiometric_afr) * enrichment; diff --git a/src/engine.cpp b/src/engine.cpp index 6fcc378..b108f49 100644 --- a/src/engine.cpp +++ b/src/engine.cpp @@ -50,7 +50,7 @@ float EngineCylinder::clearance_volume() const { } float EngineCylinder::total_mass() const { - return air_mass + exhaust_mass; + return air_mass + exhaust_mass + fuel_mass; } float EngineCylinder::temperature() const { @@ -79,11 +79,13 @@ void EngineCylinder::initialize(const Crankshaft& /*crank*/) { volume = clearance_volume(); air_mass = 1.225f * volume; exhaust_mass = 0.0f; + fuel_mass = 0.0f; energy = 101325.0f * volume / (gamma - 1.0f); pressure = 101325.0f; intake.owner_energy = exhaust.owner_energy = &energy; intake.owner_air_mass = exhaust.owner_air_mass = &air_mass; intake.owner_exhaust_mass = exhaust.owner_exhaust_mass = &exhaust_mass; + intake.owner_fuel_mass = exhaust.owner_fuel_mass = &fuel_mass; sync_ports(); } @@ -140,8 +142,10 @@ void FourStrokeCylinder::pre_step(Crankshaft& crank, float dt, float /*throttle* const float previous_degrees = wrap_degrees((crank.previous_angle + phase) * 180.0f / pi); const float current_degrees = wrap_degrees(theta * 180.0f / pi); if (crossed_window_close(previous_degrees, current_degrees, intake_open, intake_close)) { - fuel_mass = ecu ? ecu->fuel_for_air(air_mass) : air_mass / std::max(afr, 1e-6f); - fuel_injected = fuel_mass > 0.0f; + if (!transported_fuel) + fuel_mass = ecu ? ecu->fuel_for_air(air_mass) : air_mass / std::max(afr, 1e-6f); + const float charge_afr = air_mass / std::max(fuel_mass, 1e-12f); + fuel_injected = fuel_mass > 0.0f && (!transported_fuel || (charge_afr >= 7 && charge_afr <= 24)); } const float spark_angle = wrap_degrees(720.0f - spark_advance); @@ -149,6 +153,7 @@ void FourStrokeCylinder::pre_step(Crankshaft& crank, float dt, float /*throttle* if (crossed_spark && fuel_injected && !combustion_active && !(ecu ? ecu->fuel_cut() : rev_limit_active) && random_unit() >= misfire_probability) { combustion_active = true; burn_fraction = 0.0f; + active_fuel_mass_ = transported_fuel ? std::min(fuel_mass, air_mass / 14.7f) : fuel_mass; energy_factor_ = 1.0f + energy_variation * (2.0f * random_unit() - 1.0f); } @@ -161,14 +166,23 @@ void FourStrokeCylinder::pre_step(Crankshaft& crank, float dt, float /*throttle* if (fraction >= 1.0f || angle_since_spark > burn_duration + burn_start + spark_advance) { fraction = 1.0f; combustion_active = false; - exhaust_mass += air_mass; - air_mass = 0.0f; + if (!transported_fuel) { + exhaust_mass += air_mass; + air_mass = 0.0f; + } } fuel_injected = false; const float increment = fraction - burn_fraction; if (increment > 0.0f) { - energy += fuel_mass * fuel_lhv * energy_factor_ * increment; - exhaust_mass += fuel_mass * increment; + const float burned_fuel = std::min(active_fuel_mass_ * increment, fuel_mass); + energy += burned_fuel * fuel_lhv * energy_factor_; + fuel_mass -= burned_fuel; + exhaust_mass += burned_fuel; + if (transported_fuel) { + const float burned_air = std::min(air_mass, burned_fuel * 14.7f); + air_mass -= burned_air; + exhaust_mass += burned_air; + } burn_fraction = fraction; } } @@ -184,16 +198,20 @@ void FourStrokeCylinder::pre_step(Crankshaft& crank, float dt, float /*throttle* } void EngineCylinder::update_state(float dt) { - float delta_air = 0.0f, delta_exhaust = 0.0f, delta_energy = 0.0f; + float delta_air = 0.0f, delta_exhaust = 0.0f, delta_fuel = 0.0f, delta_energy = 0.0f; for (Port* port : {&intake, &exhaust}) { const float mass_flow = port->mass_flow; const float fraction = mass_flow >= 0.0f ? port->air_fraction : air_fraction(); + const float fuel_fraction = mass_flow >= 0.0f ? port->fuel_fraction : + fuel_mass / std::max(total_mass(), 1e-12f); delta_air += mass_flow * fraction * dt; - delta_exhaust += mass_flow * (1.0f - fraction) * dt; + delta_fuel += mass_flow * fuel_fraction * dt; + delta_exhaust += mass_flow * std::max(1.0f - fraction - fuel_fraction, 0.0f) * dt; delta_energy += mass_flow * port->enthalpy * dt; } air_mass += delta_air; exhaust_mass += delta_exhaust; + fuel_mass += delta_fuel; energy += delta_energy; const float current_volume = std::max(volume, 1e-12f); @@ -211,12 +229,14 @@ void EngineCylinder::update_state(float dt) { if (total_mass() < 1e-9f) { air_mass = 1e-9f; exhaust_mass = 0.0f; + fuel_mass = 0.0f; energy = 101325.0f * current_volume / (gamma - 1.0f); } else if (energy < 0.0f) { energy = 101325.0f * current_volume / (gamma - 1.0f); } air_mass = std::max(air_mass, 0.0f); exhaust_mass = std::max(exhaust_mass, 0.0f); + fuel_mass = std::max(fuel_mass, 0.0f); pressure = (gamma - 1.0f) * energy / current_volume; sync_ports(); } @@ -257,12 +277,14 @@ void EngineCylinder::sync_ports() { pressure = (gamma - 1.0f) * energy / std::max(volume, 1e-12f); const float temp = pressure / std::max(density * gas_constant, 1e-12f); const float fraction = air_fraction(); + const float fuel_fraction = fuel_mass / std::max(total_mass(), 1e-12f); const float enthalpy = gamma / (gamma - 1.0f) * pressure / std::max(density, 1e-12f); for (Port* port : {&intake, &exhaust}) { port->pressure = pressure; port->density = density; port->temperature = temp; port->air_fraction = fraction; + port->fuel_fraction = fuel_fraction; port->enthalpy = enthalpy; } } diff --git a/src/main.cpp b/src/main.cpp index 263aecc..fd4ecdd 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -46,10 +46,10 @@ int main(int argc, char** argv) { constexpr int sample_rate = 44100; const auto total_samples = static_cast(duration_seconds * sample_rate); const auto report_interval = static_cast(sample_rate / 10); - if (scenario_name == "twin" || scenario_name == "b205") { - FourStrokeScenario scenario(sample_rate, scenario_name == "b205" ? FourStrokePreset::SaabB205 : FourStrokePreset::Twin, turbo_enabled); + if (scenario_name == "twin" || scenario_name == "b205" || scenario_name == "v8") { + FourStrokeScenario scenario(sample_rate, scenario_name == "v8" ? FourStrokePreset::AmericanV8 : scenario_name == "b205" ? FourStrokePreset::SaabB205 : FourStrokePreset::Twin, turbo_enabled); scenario.throttle = throttle; - std::cout << (scenario_name == "b205" ? "FluidSimCpp | Saab B205 / 9-5 wagon | throttle " : "FluidSimCpp | 600 cc parallel twin | throttle ") << std::fixed << std::setprecision(0) + std::cout << (scenario_name == "v8" ? "FluidSimCpp | American 5.0 V8 | throttle " : scenario_name == "b205" ? "FluidSimCpp | Saab B205 / 9-5 wagon | throttle " : "FluidSimCpp | 600 cc parallel twin | throttle ") << std::fixed << std::setprecision(0) << throttle * 100.0f << "%\n"; for (std::uint64_t i = 0; i < total_samples; ++i) { scenario.process(); @@ -116,7 +116,7 @@ int main(int argc, char** argv) { } } else { std::cerr << "Unknown scenario '" << scenario_name - << "'. Choose twin, two-stroke, b205, helmholtz, or pulse.\n"; + << "'. Choose twin, two-stroke, b205, v8, helmholtz, or pulse.\n"; return 2; } return 0; diff --git a/src/outdoor_audio.cpp b/src/outdoor_audio.cpp index 2d5a770..fa97ad4 100644 --- a/src/outdoor_audio.cpp +++ b/src/outdoor_audio.cpp @@ -47,6 +47,7 @@ OutdoorAudio::OutdoorAudio(int rate, OutdoorAudioConfig config) : rate_(rate), c config_.listener_right.z / norm}; configure_emitter(emitters_[0], config_.intake_position); configure_emitter(emitters_[1], config_.exhaust_position); + if (config_.dual_exhaust) configure_emitter(emitters_[2], config_.second_exhaust_position); constexpr std::array lengths{1499, 2111, 2633, 3251}; for (std::size_t i = 0; i < tail_.size(); ++i) { auto &line = tail_[i]; @@ -148,8 +149,8 @@ AudioFrame OutdoorAudio::diffuse(float input) { } return {(values[0] + values[2]) * .35f, (values[1] + values[3]) * .35f}; } -AudioFrame OutdoorAudio::process(float intake, float exhaust) { - if (!std::isfinite(intake) || !std::isfinite(exhaust)) +AudioFrame OutdoorAudio::process(float intake, float exhaust, float second_exhaust) { + if (!std::isfinite(intake) || !std::isfinite(exhaust) || !std::isfinite(second_exhaust)) throw std::runtime_error("nonfinite acoustic source flow"); const bool wet = scene_ != OutdoorScene::Dry; const std::array target{1, wet ? 1.f : 0.f, scene_ == OutdoorScene::Wall ? 1.f : 0.f, @@ -165,7 +166,8 @@ AudioFrame OutdoorAudio::process(float intake, float exhaust) { last_raw_ = {}; const auto a = render(emitters_[0], intake, last_dry_, last_raw_); const auto b = render(emitters_[1], exhaust, last_dry_, last_raw_); + const auto c = config_.dual_exhaust ? render(emitters_[2], second_exhaust, last_dry_, last_raw_) : AudioFrame{}; const auto tail = diffuse((last_dry_.left + last_dry_.right) * .5f); - return {a.left + b.left + tail.left * tail_mix_, a.right + b.right + tail.right * tail_mix_}; + return {a.left + b.left + c.left + tail.left * tail_mix_, a.right + b.right + c.right + tail.right * tail_mix_}; } } // namespace fluidsim::port diff --git a/src/scenario.cpp b/src/scenario.cpp index 8e3e95e..c86a7bb 100644 --- a/src/scenario.cpp +++ b/src/scenario.cpp @@ -7,11 +7,17 @@ namespace fluidsim::port { FourStrokeScenario::FourStrokeScenario(int sample_rate, FourStrokePreset selected_preset, bool turbo_enabled) : preset(selected_preset), - intake_plenum(selected_preset == FourStrokePreset::SaabB205 ? 3e-3f : 1e-3f, 101325, 300), - exhaust_collector(selected_preset == FourStrokePreset::SaabB205 ? .6e-3f : .15e-3f, 101325, 700), + intake_plenum(selected_preset == FourStrokePreset::AmericanV8 ? 4e-3f : selected_preset == FourStrokePreset::SaabB205 ? 3e-3f : 1e-3f, 101325, 300), + exhaust_collector(selected_preset == FourStrokePreset::AmericanV8 ? 1.5e-3f : selected_preset == FourStrokePreset::SaabB205 ? .6e-3f : .15e-3f, 101325, 700), pipes(create_pipes(selected_preset, turbo_enabled)), dt_(1.f / std::max(sample_rate, 1)), - boundaries_(pipes), solver_(pipes, boundaries_), audio_(sample_rate), mechanical_sound_(sample_rate) { + boundaries_(pipes), solver_(pipes, boundaries_), audio_(sample_rate, [&] { + OutdoorAudioConfig config; + config.dual_exhaust = selected_preset == FourStrokePreset::AmericanV8; + return config; + }()), mechanical_sound_(sample_rate) { const bool saab = preset == FourStrokePreset::SaabB205; + const bool v8 = preset == FourStrokePreset::AmericanV8; + if (v8 && turbo_enabled) throw std::invalid_argument("American V8 preset is naturally aspirated"); const int count = cylinder_count(); crankshaft.omega = (saab ? 1200 : 2000) * 2 * pi / 60; crankshaft.inertia = saab ? .18f : .04f; @@ -36,13 +42,54 @@ FourStrokeScenario::FourStrokeScenario(int sample_rate, FourStrokePreset selecte ecu.overrun_cut_rpm = 2000; ecu.overrun_resume_rpm = 1300; } + if (v8) { + // Estimated mild 1980s 302-class pushrod small block and manual coupe. + crankshaft.omega = 1000 * 2*pi/60; + crankshaft.inertia = 1.1f; + crankshaft.friction_constant = 22; + crankshaft.friction_viscous = .025f; + vehicle.mass = 3000; + vehicle.primary_reduction = 1; + vehicle.final_drive_ratio = 6.16f / 1.7f; // 1.7x taller gearing in every gear. + vehicle.gear_ratios = {3.35f, 1.93f, 1.29f, 1.f, .68f}; + vehicle.wheel_radius = .32f; + vehicle.clutch_capacity = 1100; // Hold V8 torque with the heavy flywheel/car. + vehicle.drag_coefficient = .4f; + vehicle.frontal_area = 2.f; + ecu.idle_rpm = 750; + ecu.throttle_time_constant = .065f; + ecu.throttle_rate = 2.3f; + ecu.throttle_closing_time_constant = .16f; + ecu.throttle_closing_rate = 1.4f; + ecu.idle_base_opening = .08f; + ecu.idle_spark_gain = .02f; + ecu.ignition_map = {{0, 10}, {750, 12}, {2500, 28}, {4000, 32}, {6000, 30}}; + ecu.rev_limit_rpm = 7000; + ecu.rev_resume_rpm = 6750; + ecu.overrun_cut_rpm = 2000; + ecu.overrun_resume_rpm = 1100; + // One minimum substep; acoustic CFL still subdivides when needed. + solver_.minimum_substeps = 1; + } vehicle.set_neutral(); vehicle.clutch_input = clutch; const auto engines = cylinders(); - const std::array phases{0, pi, 3*pi, 2*pi}; // I4 firing: 1-3-4-2 + const std::array phases{0, pi, 3*pi, 2*pi}; // I4 firing: 1-3-4-2 + const std::array firing_order{1, 5, 4, 2, 6, 3, 7, 8}; for (int i = 0; i < count; ++i) { auto& c = *engines[i]; c.phase = phases[i]; + if (v8) { + const auto position = std::find(firing_order.begin(), firing_order.end(), i+1) - firing_order.begin(); + c.phase = std::fmod(4*pi - float(position)*pi*.5f, 4*pi); + c.bore = .1016f; + c.stroke = .0762f; + c.rod = .129f; + c.compression = 7.0f; + c.intake_open = 350; c.intake_close = 570; + c.exhaust_open = 140; c.exhaust_close = 370; + c.burn_duration = 65; + } if (saab) { c.bore = .090f; c.stroke = .078f; @@ -54,13 +101,17 @@ FourStrokeScenario::FourStrokeScenario(int sample_rate, FourStrokePreset selecte c.intake_valve_lift = saab ? .0085f : .008f; c.exhaust_valve_diameter = saab ? .028f : .026f; c.exhaust_valve_lift = saab ? .008f : .0075f; + if (v8) { + c.intake_valve_diameter = .045f; c.exhaust_valve_diameter = .037f; + c.intake_valve_lift = .010f; c.exhaust_valve_lift = .010f; + } c.rev_limit_enabled = true; c.rev_limit_rpm = ecu.rev_limit_rpm; c.rev_limit_hysteresis = 250; c.rev_limit_cut_time = .020f; c.ecu = &ecu; c.initialize(crankshaft); - if (saab) { + if (saab || v8) { // Initialize each phased cylinder at its actual piston position. c.volume = c.compute_volume(c.phase); c.air_mass = ambient_density * c.volume; @@ -80,6 +131,7 @@ FourStrokeScenario::FourStrokeScenario(int sample_rate, FourStrokePreset selecte solver_.add_component(intake_atmosphere); solver_.add_component(intake_plenum); solver_.add_component(exhaust_collector); + if (v8) solver_.add_component(second_bank_collector); if (turbo_enabled) { turbo = std::make_unique(saab ? TurboConfig::gt17() : TurboConfig{}); ecu.boost_control_enabled = saab; @@ -93,7 +145,8 @@ FourStrokeScenario::FourStrokeScenario(int sample_rate, FourStrokePreset selecte boundaries_.add_orifice(throttle_port, 0, false, throttle_area_, .7f); for (int i = 0; i < count; ++i) { auto& feed = intake_plenum.create_port(); - auto& collector_in = exhaust_collector.create_port(); + // 1--4 and 5--8 are separate physical banks for this firing order. + auto& collector_in = (v8 && i >= 4 ? second_bank_collector : exhaust_collector).create_port(); solver_.add_component(*engines[i]); boundaries_.add_orifice(feed, 1+i, true, runner_area_); boundaries_.add_orifice(engines[i]->intake, 1+i, false, intake_valve_areas_[i]); @@ -104,11 +157,16 @@ FourStrokeScenario::FourStrokeScenario(int sample_rate, FourStrokePreset selecte boundaries_.add_orifice(collector_out, 1+2*count, true, tail_area_); const int exhaust_exit = turbo_enabled ? 3+2*count : 1+2*count; boundaries_.add_open_end(exhaust_exit, false, 101325, pipes[exhaust_exit].area); + if (v8) { + auto& outlet = second_bank_collector.create_port(); + boundaries_.add_orifice(outlet, 18, true, tail_area_); + boundaries_.add_open_end(18, false, 101325, pipes[18].area); + } } AudioFrame FourStrokeScenario::process() { const std::lock_guard guard(mutex_); - step_motorcycle(crankshaft, dt_, preset == FourStrokePreset::SaabB205 ? 200.f : 30.f); + step_motorcycle(crankshaft, dt_, preset == FourStrokePreset::AmericanV8 ? 300.f : preset == FourStrokePreset::SaabB205 ? 200.f : 30.f); ecu.step(dt_, throttle, crankshaft.omega * 60 / (2*pi)); if (turbo) { ecu.control_boost(dt_, throttle, crankshaft.omega * 60/(2*pi), intake_plenum.pressure(), turbo->compressor_outlet.pressure()); @@ -122,12 +180,15 @@ AudioFrame FourStrokeScenario::process() { intake_valve_areas_[i] = c.valve_area(degrees, true); exhaust_valve_areas_[i] = c.valve_area(degrees, false); } - const float angle = ecu.throttle_opening() * pi * .5f; - throttle_area_ = area_for_diameter(preset == FourStrokePreset::SaabB205 ? .052f : .026f) * (1-std::cos(angle)); + const float opening = ecu.throttle_opening(); + const float angle = opening * pi * .5f; + throttle_area_ = area_for_diameter(preset == FourStrokePreset::AmericanV8 ? .065f : preset == FourStrokePreset::SaabB205 ? .052f : .026f) * (1-std::cos(angle)); solver_.step(dt_); intake_flow = intake_atmosphere.port.mass_flow; exhaust_flow = boundaries_.open_end_mass_flow(0); - const auto rendered = audio_.process(intake_flow, exhaust_flow); + const float other_exhaust = preset == FourStrokePreset::AmericanV8 ? boundaries_.open_end_mass_flow(1) : 0.f; + const auto rendered = audio_.process(intake_flow, exhaust_flow, other_exhaust); + exhaust_flow += other_exhaust; const float mechanical = mechanical_sound_enabled ? mechanical_sound_.process(crankshaft, cylinder1, load) * mechanical_mix : 0; last_mix_l = rendered.left + mechanical; @@ -222,7 +283,8 @@ ScenarioSnapshot FourStrokeScenario::snapshot() const { result.mechanical_sound_enabled = mechanical_sound_enabled; result.cylinder_count = cylinder_count(); result.rev_limit_rpm = ecu.rev_limit_rpm; - const std::array engines{&cylinder1, &cylinder2, &cylinder3, &cylinder4}; + const std::array engines{&cylinder1, &cylinder2, &cylinder3, &cylinder4, + &cylinder5, &cylinder6, &cylinder7, &cylinder8}; for (int i = 0; i < cylinder_count(); ++i) { result.cylinder_pressures[i] = engines[i]->pressure; result.piston_fractions[i] = (engines[i]->volume-engines[i]->clearance_volume())/engines[i]->swept_volume(); @@ -239,6 +301,20 @@ float FourStrokeScenario::area_for_diameter(float diameter) { std::vector FourStrokeScenario::create_pipes(FourStrokePreset preset, bool turbo_enabled) { const bool saab = preset == FourStrokePreset::SaabB205; + if (preset == FourStrokePreset::AmericanV8) { + // One inlet/throttle -> shared plenum -> short individual runners. + // Independent exhaust banks. Cells stay long enough to avoid routine + // extra CFL substeps during high-vacuum coast-down. + std::vector result; + result.emplace_back(area_for_diameter(.070f), .18f, 3); + for (int i = 0; i < 8; ++i) + result.emplace_back(area_for_diameter(.038f), .12f, 2); + for (int i = 0; i < 8; ++i) + result.emplace_back(area_for_diameter(.040f), .60f + .06f*(i%4), 12); + result.emplace_back(area_for_diameter(.050f), 1.8f, 36); + result.emplace_back(area_for_diameter(.050f), 1.95f, 39); + return result; + } const int count = saab ? 4 : 2; std::vector result; result.reserve(4+2*count); diff --git a/src/turbocharger.cpp b/src/turbocharger.cpp index 5848709..c78706d 100644 --- a/src/turbocharger.cpp +++ b/src/turbocharger.cpp @@ -106,6 +106,7 @@ void Turbocharger::set_transfer(Volume0D& in, Volume0D& out, Port& ip, Port& op, ip.mass_flow = float(-flow); op.mass_flow = float(flow); ip.air_fraction = op.air_fraction = donor.air_fraction(); + ip.fuel_fraction = op.fuel_fraction = donor.fuel_fraction(); ip.enthalpy = op.enthalpy = float(h); // Positive work means energy delivered to the shaft. Only the receiving // fluid enthalpy changes; the donor loses precisely its stagnation enthalpy. diff --git a/src/two_stroke_cylinder.cpp b/src/two_stroke_cylinder.cpp index 9585a49..18c109a 100644 --- a/src/two_stroke_cylinder.cpp +++ b/src/two_stroke_cylinder.cpp @@ -1,3 +1,4 @@ +#include "FluidSim/Carburetor.hpp" #include "FluidSim/Ecu.hpp" #include "FluidSim/TwoStrokeCylinder.hpp" @@ -70,16 +71,22 @@ void TwoStrokeCylinder::pre_step(Crankshaft &crank, float dt, float /*throttle*/ void TwoStrokeCylinder::handle_cycle_events(float previous_degrees, float current_degrees) { const float trapped_angle = std::max(transfer_close(), exhaust_close()); if (previous_degrees < trapped_angle && current_degrees >= trapped_angle) { - fuel_mass = ecu ? ecu->fuel_for_air(air_mass) : air_mass / std::max(afr, 1e-6f); - fuel_injected = fuel_mass > 0; + if (!carburetor_ && fuel_mass <= 0.0f) + fuel_mass = ecu ? ecu->fuel_for_air(air_mass) : air_mass / std::max(afr, 1e-6f); + const float charge_afr = air_mass / std::max(fuel_mass, 1e-12f); + fuel_injected = fuel_mass > 0.0f && charge_afr >= 7.0f && charge_afr <= 24.0f; } const float spark_angle = wrap_degrees(360.0f - spark_advance, 360.0f); const bool crossed_spark = crossed_angle(previous_degrees,current_degrees,spark_angle); - if (crossed_spark && fuel_injected && !combustion_active && !(ecu ? ecu->fuel_cut() : rev_limit_active) && + if (crossed_spark && fuel_injected && !combustion_active && + !(ecu ? ecu->fuel_cut() : rev_limit_active) && random_unit() >= misfire_probability) { combustion_active = true; burn_fraction = 0.0f; + // Fuel is a transported species. Only the portion with enough oxygen + // can release heat; rich excess remains unburned and can reach the pipe. + active_fuel_mass_ = std::min(fuel_mass, air_mass / 14.7f); ++combustion_count_; active_spark_angle_ = spark_angle; // Freeze the duration and timing for this charge. Changing the ECU map @@ -99,14 +106,16 @@ void TwoStrokeCylinder::handle_cycle_events(float previous_degrees, float curren if (fraction >= 1.0f || angle_since_spark > active_burn_duration_ + burn_start) { fraction = 1.0f; combustion_active = false; - exhaust_mass += air_mass; - air_mass = 0.0f; } fuel_injected = false; const float increment = fraction - burn_fraction; if (increment > 0.0f) { - energy += fuel_mass * fuel_lhv * energy_factor_ * increment; - exhaust_mass += fuel_mass * increment; + const float burned_fuel = std::min(active_fuel_mass_ * increment, fuel_mass); + energy += burned_fuel * fuel_lhv * energy_factor_; + fuel_mass -= burned_fuel; + const float burned_air = std::min(air_mass, burned_fuel * 14.7f); + air_mass -= burned_air; + exhaust_mass += burned_fuel + burned_air; burn_fraction = fraction; } } diff --git a/src/two_stroke_scenario.cpp b/src/two_stroke_scenario.cpp index a71bea2..1092d17 100644 --- a/src/two_stroke_scenario.cpp +++ b/src/two_stroke_scenario.cpp @@ -7,6 +7,7 @@ namespace fluidsim::port { TwoStrokeScenario::TwoStrokeScenario(int sample_rate, TwoStrokeConfig config) : configuration(config), + carburetor(config.carburetor_diameter), cylinder(config.bore, config.stroke, config.rod, config.geometric_compression(), config.transfer_duration, config.exhaust_duration), crankcase(crankshaft, config.stroke * .5f, config.rod, config.bore, 150e-6f), pipes(create_pipes(config)), @@ -43,9 +44,22 @@ TwoStrokeScenario::TwoStrokeScenario(int sample_rate, TwoStrokeConfig config) ecu.idle_spark_gain=.04f; ecu.overrun_cut_rpm = 5000; ecu.overrun_resume_rpm = 3500; + ecu.overrun_fuel_cut_enabled = false; cylinder.ecu = &ecu; cylinder.initialize(crankshaft); cylinder.set_crankcase(crankcase); + cylinder.set_carburetor(carburetor); + // This scenario starts as an already-running warm engine at 3000 RPM. + // Seed its upstream charge accordingly; a dry system needs a starter. + crankcase.initialize_mixture(14.7f); + cylinder.fuel_mass = cylinder.air_mass / 15.7f; + cylinder.air_mass -= cylinder.fuel_mass; + intake_plenum.add_fuel_mass(intake_plenum.mass() / 14.7f); + for (int index : {1, 2}) + for (auto& cell : pipes[index].cells) { + cell.fuel = 1.0f / 15.7f; + cell.air = 1.0f - cell.fuel; + } for (auto &pipe : pipes) { pipe.damping_multiplier = 0.8f; @@ -77,8 +91,12 @@ AudioFrame TwoStrokeScenario::process() { reed_out_area_ = std::min(reed.step(dt_, upstream - crankcase.pressure()), reed_in_area_); reed_open_ = reed.lift() > 1e-6f; ecu.step(dt_, throttle, crankshaft.omega * 60 / (2 * pi)); - const float angle = ecu.throttle_opening() * pi * .5f; - max_throttle_area_ = area_for_diameter(configuration.carburetor_diameter) * (1 - std::cos(angle)); + const auto carb_inlet = pipes[0].state(static_cast(pipes[0].cells.size()) - 1); + carburetor.step(dt_, throttle, ecu.throttle_opening(), crankshaft.omega * 60 / (2 * pi), + std::max(plenum_inlet_->mass_flow, 0.0f), carb_inlet.rho, + carb_inlet.pressure, intake_plenum.pressure()); + intake_plenum.add_fuel_mass(carburetor.fuel_flow() * dt_); + max_throttle_area_ = slide_area(configuration.carburetor_diameter, ecu.throttle_opening()); step_motorcycle(crankshaft, dt_, 10.0f); crankcase.pre_step(dt_); cylinder.pre_step(crankshaft, dt_, ecu.throttle_opening()); @@ -116,6 +134,16 @@ TwoStrokeSnapshot TwoStrokeScenario::snapshot() const { float TwoStrokeScenario::area_for_diameter(float diameter) { return pi * diameter * diameter * 0.25f; } +float TwoStrokeScenario::slide_area(float diameter, float opening) { + const float radius = diameter * 0.5f; + const float height = std::clamp(opening, 0.0f, 1.0f) * diameter; + const float offset = radius - height; + const float root = std::sqrt(std::max(radius * radius - offset * offset, 0.0f)); + // Cancellation near the seat can produce a tiny negative area in float. + return std::clamp(radius * radius * std::acos(std::clamp(offset / radius, -1.0f, 1.0f)) - offset * root, + 0.0f, area_for_diameter(diameter)); +} + std::vector TwoStrokeScenario::create_pipes(const TwoStrokeConfig& config) { std::vector result; result.reserve(4); diff --git a/tests/fluid_validation.cpp b/tests/fluid_validation.cpp index 0be24ff..43247b0 100644 --- a/tests/fluid_validation.cpp +++ b/tests/fluid_validation.cpp @@ -15,8 +15,8 @@ void check(bool b, const std::string &s) { void near(double a, double b, double tol, const char *s) { check(std::isfinite(a) && std::abs(a - b) <= tol, s); } -void set(Pipe1D &p, int i, float pressure, float rho, float u = 0, float air = 1) { - p.cells[i] = {rho, rho * u, pressure / pipe_gm1 + .5f * rho * u * u, air}; +void set(Pipe1D &p, int i, float pressure, float rho, float u = 0, float air = 1, float fuel = 0) { + p.cells[i] = {rho, rho * u, pressure / pipe_gm1 + .5f * rho * u * u, air, fuel}; } void lossless(Pipe1D &p) { p.damping_multiplier = 0; @@ -499,7 +499,8 @@ int main() { for (int k = 0; k < 2; k++) { lossless(p[k]); for (int i = 0; i < 40; i++) - set(p[k], i, k ? 100000 : 120000, 1.2f, 0, k ? .9f : .1f); + set(p[k], i, k ? 100000 : 120000, 1.2f, 0, k ? .9f : .1f, + k ? .02f : .05f); } auto air_mass = [&] { double sum = 0; @@ -508,7 +509,15 @@ int main() { sum += q.cell_mass(i) * q.cells[i].air; return sum; }; + auto fuel_mass = [&] { + double sum = 0; + for (auto &q : p) + for (int i = 0; i < 40; i++) + sum += q.cell_mass(i) * q.cells[i].fuel; + return sum; + }; double before = air_mass(); + double fuel_before = fuel_mass(); BoundarySystem b(p); b.add_closed_end(0, true); b.add_closed_end(1, false); @@ -517,6 +526,7 @@ int main() { for (int i = 0; i < 3000; i++) solver.step(1.f / 44100); near(air_mass(), before, before * 4e-5, "species created at shared face"); + near(fuel_mass(), fuel_before, fuel_before * 4e-5, "fuel species created at shared face"); }); test("unsupported unequal junction rejected", [] { std::vector p; diff --git a/tests/two_stroke_physics.cpp b/tests/two_stroke_physics.cpp index 3eac4ac..b267434 100644 --- a/tests/two_stroke_physics.cpp +++ b/tests/two_stroke_physics.cpp @@ -1,3 +1,4 @@ +#include "FluidSim/Carburetor.hpp" #include "FluidSim/TwoStrokeScenario.hpp" #include #include @@ -37,6 +38,17 @@ float reflected_pulse(bool expanding) { } } int main() { + Carburetor carburetor; + for (int i = 0; i < 1000; ++i) + carburetor.step(.001f, .55f, .55f, 8000.0f, .012f, 1.15f, 101325.0f, 90000.0f); + const float cruise_afr = carburetor.indicated_afr(); + carburetor.step(.001f, 0.0f, .05f, 8000.0f, .004f, 1.15f, 101325.0f, 86000.0f); + check(carburetor.overrun_active() && carburetor.fuel_flow() > 0.0f, + "carburetor cut fuel when the throttle closed"); + check(std::isfinite(cruise_afr) && cruise_afr >= 4.0f && cruise_afr <= 40.0f && + carburetor.bowl_level() > 0.0f && carburetor.bowl_level() <= 1.0f, + "carburetor mixture or float bowl left its physical bounds"); + const std::array cone{{{.4f, .03f, .10f}}}; const float expected = pi/12 * .4f * (.03f*.03f + .03f*.10f + .10f*.10f); for (float mesh : {.03f, .01f, .004f}) { @@ -84,6 +96,16 @@ int main() { check(std::abs(first-30.27f)<.03f && std::abs(s.vehicle.total_gear_ratio()-13.4f)<.03f, "YZ125 primary/final/gear reduction incorrect"); + s.set_throttle(.5f); + for (int sample = 0; sample < 22050; ++sample) + s.process(); + bool transported_fuel = false; + for (int pipe_index : {1, 2}) + for (const auto& cell : s.pipes[pipe_index].cells) + transported_fuel |= cell.fuel > 0.0f; + check(s.intake_plenum.fuel_mass() > 0.0f && s.crankcase.fuel_mass() > 0.0f && transported_fuel, + "carburetor fuel was not transported through intake, reed and crankcase"); + Crankshaft crank; crank.omega = 100; crank.load_torque = 10; diff --git a/tests/two_stroke_regression.cpp b/tests/two_stroke_regression.cpp index e4c0752..d713c56 100644 --- a/tests/two_stroke_regression.cpp +++ b/tests/two_stroke_regression.cpp @@ -12,9 +12,17 @@ int main() { bool reed_open_seen = false, reed_closed_seen = false, above_old_limit = false; DynoCurve dyno; float peak = 0, maximum_rpm = 0; + std::uint64_t release_firings = 0, one_second_coast_firings = 0; + float release_rpm = 0; for (int sample = 0; sample < 14 * 44100; ++sample) { scenario.set_throttle(sample >= 2 * 44100 && sample < 6 * 44100 ? .6f : 0.f); const auto audio = scenario.process(); + if (sample == 6 * 44100) { + release_firings = scenario.cylinder.combustion_count(); + release_rpm = scenario.snapshot().rpm; + } else if (sample == 7 * 44100) { + one_second_coast_firings = scenario.cylinder.combustion_count() - release_firings; + } if (!std::isfinite(audio.left) || !std::isfinite(audio.right)) throw std::runtime_error("Nonfinite two-stroke audio"); peak = std::max({peak, std::abs(audio.left), std::abs(audio.right)}); @@ -31,6 +39,8 @@ int main() { throw std::runtime_error("Two-stroke cut fuel with the pedal open"); above_old_limit |= state.rpm > 10200; } + if (sample >= 6 * 44100 && sample < 7 * 44100 && state.rpm > 5000 && state.fuel_cut) + throw std::runtime_error("Carbureted two-stroke cut fuel on overrun"); dyno.update(state.rpm, scenario.crankshaft.average_power / 1000, scenario.crankshaft.average_torque); reed_open_seen |= state.reed_open; @@ -39,11 +49,16 @@ int main() { } const auto final = scenario.snapshot(); if (!above_old_limit) - throw std::runtime_error("Two-stroke did not exceed the old rev limit"); + throw std::runtime_error("Two-stroke did not exceed the old rev limit: peak=" + + std::to_string(maximum_rpm) + ", idle=" + std::to_string(final.rpm)); if (!reed_open_seen || !reed_closed_seen) throw std::runtime_error("Reed did not open and close"); if (final.rpm < 2000 || final.rpm > 4000) - throw std::runtime_error("Two-stroke idle did not recover after throttle release"); + throw std::runtime_error("Two-stroke idle did not recover after throttle release: " + + std::to_string(final.rpm) + " RPM, coast firings=" + + std::to_string(one_second_coast_firings) + ", carb AFR=" + + std::to_string(scenario.carburetor.indicated_afr()) + + ", cylinder fuel=" + std::to_string(scenario.cylinder.fuel_mass)); if (peak >= 1) throw std::runtime_error("Two-stroke audio exceeds device range"); for (const auto &pipe : scenario.pipes) @@ -83,5 +98,6 @@ int main() { throw std::runtime_error("Motorcycle launch required pipe positivity repairs"); std::cout << "Two-stroke regression passed: peak=" << peak << ", maximum RPM=" << maximum_rpm << ", recovered RPM=" << final.rpm << ", bike speed=" << scenario.vehicle.speed * 3.6f - << " km/h\n"; + << " km/h, release RPM=" << release_rpm + << ", first coast-second firings=" << one_second_coast_firings << '\n'; } diff --git a/tools/audio_capture.cpp b/tools/audio_capture.cpp index b1c2b6d..5ebb442 100644 --- a/tools/audio_capture.cpp +++ b/tools/audio_capture.cpp @@ -159,8 +159,10 @@ int main(int argc, char **argv) { capture(prefix, frames, pedal, selected, config); else if (scenario == "two-stroke") capture(prefix, frames, pedal, selected, config); + else if (scenario == "v8") + capture(prefix, frames, pedal, selected, config); else - throw std::runtime_error("choose twin, two-stroke or b205"); + throw std::runtime_error("choose twin, two-stroke, b205 or v8"); } catch (const std::exception &e) { std::cerr << e.what() << '\n'; return 1; diff --git a/tools/scenario_benchmark.cpp b/tools/scenario_benchmark.cpp index b4efe38..4802c09 100644 --- a/tools/scenario_benchmark.cpp +++ b/tools/scenario_benchmark.cpp @@ -4,18 +4,26 @@ #include "FluidSim/PulseTestScenario.hpp" #include #include +#include #include #include using namespace fluidsim::port; template void run(const char* name, S& s, int seconds, Control control) { double checksum = 0; std::uint64_t previous_iterations = 0; + std::uint64_t previous_substeps = 0; for (int second = 0; second < seconds; ++second) { + float peak = 0; const auto begin = std::chrono::steady_clock::now(); for (int i = 0; i < 44100; ++i) { control(second); auto value = s.process(); - if constexpr (requires { value.left; }) checksum += value.left; + if constexpr (requires { value.left; }) { + if (!std::isfinite(value.left) || !std::isfinite(value.right)) + throw std::runtime_error("Nonfinite benchmark audio"); + checksum += value.left; + peak = std::max({peak, std::abs(value.left), std::abs(value.right)}); + } else checksum += value; } const double elapsed = std::chrono::duration(std::chrono::steady_clock::now()-begin).count(); @@ -32,6 +40,19 @@ template void run(const char* name, S& s, int seconds, C const auto state = s.snapshot(); std::cerr << name << " iterations_per_frame=" << (state.boundary_iterations-previous_iterations)/44100.0 << '\n'; + std::size_t limiting_pipe = 0; + float stable_dt = std::numeric_limits::max(); + for (std::size_t p = 0; p < s.pipes.size(); ++p) { + const float candidate = s.pipes[p].stable_timestep(); + if (candidate < stable_dt) { stable_dt = candidate; limiting_pipe = p; } + } + std::cerr << name << " substeps_per_frame=" << (state.solver_substeps-previous_substeps)/44100.0 + << " limiting_pipe=" << limiting_pipe << " cfl_us=" << stable_dt*1e6f << '\n'; + previous_substeps = state.solver_substeps; + unsigned long long repairs = 0; + for (const auto& pipe : s.pipes) repairs += pipe.positivity_repairs; + std::cerr << name << " rpm=" << state.rpm << " power_kW=" << state.average_power/1000 + << " peak=" << peak << " pipe_repairs=" << repairs << '\n'; previous_iterations = state.boundary_iterations; } } @@ -41,18 +62,19 @@ int main(int argc, char** argv) { const std::string mode = argc > 2 ? argv[2] : "pedal"; const std::string only = argc > 3 ? argv[3] : "all"; auto pedal = [&](int second) { - return mode == "idle" ? 0.0f : mode == "transient" && (second == 0 || second >= seconds-2) ? 0.0f : .7f; + return mode == "idle" ? 0.0f : mode == "coast" ? (second >= 1 && second < 5 ? .7f : 0.f) : + mode == "transient" && (second == 0 || second >= seconds-2) ? 0.0f : .7f; }; std::cout << "scenario,second,wall_seconds,checksum\n"; - for (auto name : {"twin", "turbo", "b205", "b205-gt17"}) { + for (auto name : {"twin", "turbo", "b205", "b205-gt17", "v8"}) { if (only != "all" && only != name) continue; auto s = std::make_unique(44100, - std::string(name).starts_with("b205") ? FourStrokePreset::SaabB205 : FourStrokePreset::Twin, + std::string(name)=="v8" ? FourStrokePreset::AmericanV8 : std::string(name).starts_with("b205") ? FourStrokePreset::SaabB205 : FourStrokePreset::Twin, std::string(name)=="turbo" || std::string(name)=="b205-gt17"); s->throttle = .7f; run(name, *s, seconds, [&](int second) { s->throttle = pedal(second); - if (mode == "loaded") s->crankshaft.omega = (std::string(name)=="b205-gt17" ? 5500 : 7000) * 2*pi/60; + if (mode == "loaded") s->crankshaft.omega = (std::string(name)=="v8" ? 4500 : std::string(name)=="b205-gt17" ? 5500 : 7000) * 2*pi/60; }); } for (bool racing : {true, false}) { diff --git a/tools/two_stroke_benchmark.cpp b/tools/two_stroke_benchmark.cpp index 46863fc..99afaa4 100644 --- a/tools/two_stroke_benchmark.cpp +++ b/tools/two_stroke_benchmark.cpp @@ -1,10 +1,11 @@ #include "FluidSim/TwoStrokeScenario.hpp" +#include #include #include using namespace fluidsim::port; int main() { - std::cout << "exhaust,rpm,shaft_kW,shaft_Nm,firings_per_revolution,repairs,wall_seconds\n"; + std::cout << "exhaust,rpm,shaft_kW,shaft_Nm,firings_per_revolution,carb_afr,crankcase_afr,repairs,wall_seconds\n"; for (auto exhaust : {TwoStrokeExhaust::Straight, TwoStrokeExhaust::Racing}) { for (float rpm : {6000.f, 8000.f, 10000.f, 12000.f, 14000.f, 17000.f}) { TwoStrokeConfig config; @@ -31,6 +32,9 @@ int main() { std::cout << (exhaust == TwoStrokeExhaust::Racing ? "racing" : "straight") << ',' << rpm << ',' << power << ',' << torque << ',' << (s.cylinder.combustion_count()-initial_firings) / (rpm / 120) + << ',' << s.carburetor.indicated_afr() + << ',' << s.crankcase.mass() * s.crankcase.air_fraction() / + std::max(s.crankcase.fuel_mass(), 1e-12f) << ',' << repairs << ',' << wall << std::endl; } } diff --git a/verification/YZ125_BASELINE.md b/verification/YZ125_BASELINE.md index 60a8a39..c833edb 100644 --- a/verification/YZ125_BASELINE.md +++ b/verification/YZ125_BASELINE.md @@ -39,7 +39,11 @@ These port dimensions and compression interpretation require measured validation Combustion has an estimated .8 ms minimum burn duration; timing and duration are frozen for each active charge. The estimated ignition curve retards above 6000 RPM rather than continuing to advance toward 40 degrees. Neither is a measured -factory CDI/flame model. Fuel remains ideally calculated at port closure. +factory CDI/flame model. Fuel is now metered by effective pilot, needle and main +jet areas using simulated venturi/static pressure signals. It is transported +through the intake, crankcase and cylinder as a conserved species; the effective +jet dimensions and wall-film constants are estimates rather than measured TMX38 +calibration. Mechanical friction now includes a quadratic windage term. The two-stroke uses FMEP [bar] = .4 + .09*mean_piston_speed + .0009*mean_piston_speed^2, converted to @@ -63,6 +67,12 @@ steady speed. An analytic load-balance test checks 10 Nm and 1 kW at 100 rad/s. ## Measurements and limits +The sweep and capture numbers below are historical, before removal of scheduled +late firing and addition of exhaust oxidation. Current transient regression: +peak 12954 RPM, recovered idle 2583 RPM, rendered peak .0985, and 87 normal CDI +firings in the first coast second. No imposed pop cadence remains. Recalibration +and listening validation of the finite-rate exhaust chemistry remain necessary. + The fixed-speed benchmark imposes speed as an ideal motoring/absorbing dyno and records shaft output. It runs one simulated second per point and averages the last half. All points fired approximately once per revolution, had zero pipe @@ -70,25 +80,27 @@ positivity repairs, and cost about .53-.58 wall seconds per simulated second. | RPM | Straight kW | Racing kW | | ---: | ---: | ---: | -| 6000 | 4.02 | 6.27 | -| 8000 | .63 | 5.53 | -| 10000 | 2.43 | 10.53 | -| 12000 | 3.96 | 14.32 | -| 14000 | 3.60 | 12.38 | -| 17000 | -.80 | -3.94 | +| 6000 | 1.43 | 3.82 | +| 8000 | .41 | 4.92 | +| 10000 | 2.62 | 10.02 | +| 12000 | 3.46 | 11.87 | +| 14000 | 1.78 | 9.42 | +| 17000 | -3.65 | -5.20 | This coarse sweep peaks near 12000 RPM; it is not a validated Yamaha dyno curve. Negative power at 17000 means the virtual dyno must motor the engine. Free-rev -speed differs from peak-power speed: full pedal with the straight pipe still -approaches 16950 RPM after eight seconds. This remains an unresolved calibration -limit. The standard transient reached 16128 RPM, recovered idle at 3159 RPM, -and launched the motorcycle to 61.3 km/h without repairs or nonfinite output. +speed differs from peak-power speed and remains an unresolved calibration +limit. The standard .6-pedal racing-pipe transient reached 15954 RPM, recovered +idle at 3151 RPM, and launched the motorcycle to 57.4 km/h without repairs or +nonfinite output. Its first second after throttle release contained 11 sparse +late firings rather than continuous perfect-cycle combustion. -Six-second .6-pedal field captures had no clipping: racing rendered peak .201, -straight .554. Firing cadence, port timing and geometry are improved, but exact -YZ125 sound still requires listening and further combustion/scavenging work. +The standard transient rendered peak .395 without clipping. Firing cadence, +port timing and geometry are improved, but exact YZ125 sound still requires +listening and further combustion/scavenging work. The first-order pipe flux damps high-frequency content; the compact source, -constant-gamma gas, ideal fueling, fixed power valve, single-zone mixing and +constant-gamma gas, effective jet calibration, single-phase transported fuel, +fixed power valve, single-zone mixing and empirical monitoring/environment filters limit fidelity. No engine sound sample or synthetic two-stroke firing track has been added.