48 lines
2.0 KiB
C++
48 lines
2.0 KiB
C++
#pragma once
|
|
#include <algorithm>
|
|
#include <cmath>
|
|
#include <stdexcept>
|
|
|
|
namespace fluidsim::port {
|
|
// Lumped cantilever-petal model; geometry/material values require calibration.
|
|
class ReedValve {
|
|
public:
|
|
float length{0.030f}, width{0.016f}, thickness{0.00035f};
|
|
float youngs_modulus{70e9f}, material_density{1600.0f};
|
|
float damping_ratio{0.30f}, maximum_lift{0.004f};
|
|
int petals{4};
|
|
float lift() const { return lift_; }
|
|
float step(float dt, float pressure_difference) {
|
|
if (!std::isfinite(dt) || dt <= 0 || !std::isfinite(pressure_difference))
|
|
throw std::invalid_argument("invalid reed timestep or pressure");
|
|
if (!(length > 0 && width > 0 && thickness > 0 && youngs_modulus > 0 && material_density > 0 &&
|
|
maximum_lift > 0 && petals > 0 && damping_ratio >= 0))
|
|
throw std::invalid_argument("invalid reed geometry/material");
|
|
const float mass = .236f * material_density * length * width * thickness * petals;
|
|
const float stiffness = petals * youngs_modulus * width * thickness * thickness * thickness /
|
|
(4 * length * length * length);
|
|
const float damping = 2 * damping_ratio * std::sqrt(mass * stiffness);
|
|
const float force = pressure_difference * (3.0f / 8.0f) * length * width * petals;
|
|
const int steps = std::max(1, int(std::ceil(dt * std::sqrt(stiffness / mass) / .15f)));
|
|
const float h = dt / steps;
|
|
for (int i = 0; i < steps; i++) {
|
|
velocity_ += (force - stiffness * lift_ - damping * velocity_) / mass * h;
|
|
lift_ += velocity_ * h;
|
|
if (lift_ < 0) {
|
|
lift_ = 0;
|
|
velocity_ = 0;
|
|
}
|
|
if (lift_ > maximum_lift) {
|
|
lift_ = maximum_lift;
|
|
velocity_ = 0;
|
|
}
|
|
}
|
|
// Front-edge plus triangular side-gap areas for all petals.
|
|
return petals * (width + length) * lift_;
|
|
}
|
|
|
|
private:
|
|
float lift_{}, velocity_{};
|
|
};
|
|
} // namespace fluidsim::port
|