From 6cb999f57a0bf9f1c554e7f01b3c141530e309da Mon Sep 17 00:00:00 2001 From: Tony Martinez Date: Wed, 15 Jul 2026 17:47:01 -0600 Subject: [PATCH 01/16] First prototype of the drone interface --- src/wind_energy/actuator/ActuatorModel.H | 3 + src/wind_energy/actuator/CMakeLists.txt | 1 + src/wind_energy/actuator/actuator_types.H | 10 + src/wind_energy/actuator/drone/CMakeLists.txt | 6 + src/wind_energy/actuator/drone/Drone.H | 65 ++++ src/wind_energy/actuator/drone/Drone.cpp | 49 +++ src/wind_energy/actuator/drone/drone_ops.H | 71 ++++ src/wind_energy/actuator/drone/drone_ops.cpp | 293 ++++++++++++++++ .../actuator/sector/ActuatorSector.H | 9 + .../actuator/sector/actuator_sector_ops.H | 7 + .../actuator/sector/actuator_sector_ops.cpp | 48 ++- test/CMakeLists.txt | 1 + .../act_drone_quad/act_drone_quad.inp | 83 +++++ .../wind_energy/actuator/CMakeLists.txt | 1 + .../wind_energy/actuator/test_drone.cpp | 327 ++++++++++++++++++ 15 files changed, 973 insertions(+), 1 deletion(-) create mode 100644 src/wind_energy/actuator/drone/CMakeLists.txt create mode 100644 src/wind_energy/actuator/drone/Drone.H create mode 100644 src/wind_energy/actuator/drone/Drone.cpp create mode 100644 src/wind_energy/actuator/drone/drone_ops.H create mode 100644 src/wind_energy/actuator/drone/drone_ops.cpp create mode 100644 test/test_files/act_drone_quad/act_drone_quad.inp create mode 100644 unit_tests/wind_energy/actuator/test_drone.cpp diff --git a/src/wind_energy/actuator/ActuatorModel.H b/src/wind_energy/actuator/ActuatorModel.H index 33c0ce82d3..b67436876f 100644 --- a/src/wind_energy/actuator/ActuatorModel.H +++ b/src/wind_energy/actuator/ActuatorModel.H @@ -125,6 +125,9 @@ public: //! Return the meta info object for this actuator instance [[nodiscard]] const auto& meta() const { return m_data.meta(); } + //! Return the actuator grid for inspection and diagnostics + [[nodiscard]] const auto& grid() const { return m_data.grid(); } + void read_inputs(const utils::ActParser& pp) override { ops::ReadInputsOp()(m_data, pp); diff --git a/src/wind_energy/actuator/CMakeLists.txt b/src/wind_energy/actuator/CMakeLists.txt index e4307f0337..2b115aa85a 100644 --- a/src/wind_energy/actuator/CMakeLists.txt +++ b/src/wind_energy/actuator/CMakeLists.txt @@ -10,5 +10,6 @@ target_sources(${kynema_sgf_lib_name} add_subdirectory(aero) add_subdirectory(wing) add_subdirectory(sector) +add_subdirectory(drone) add_subdirectory(turbine) add_subdirectory(disk) diff --git a/src/wind_energy/actuator/actuator_types.H b/src/wind_energy/actuator/actuator_types.H index fdc13c5d9c..0b360aebde 100644 --- a/src/wind_energy/actuator/actuator_types.H +++ b/src/wind_energy/actuator/actuator_types.H @@ -78,6 +78,16 @@ struct ActSrcSector : ActSrcType static constexpr bool is_sector = true; }; +/** Composite source representation used by a drone containing sector rotors. */ +struct ActSrcDrone : ActSrcType +{ + static std::string identifier() { return ""; } + + static constexpr bool is_line = false; + static constexpr bool is_disk = false; + static constexpr bool is_sector = false; +}; + using RealList = amrex::Vector; using RealSlice = ::kynema_sgf::utils::Slice; using VecList = amrex::Vector; diff --git a/src/wind_energy/actuator/drone/CMakeLists.txt b/src/wind_energy/actuator/drone/CMakeLists.txt new file mode 100644 index 0000000000..612fab6c00 --- /dev/null +++ b/src/wind_energy/actuator/drone/CMakeLists.txt @@ -0,0 +1,6 @@ +target_sources(${kynema_sgf_lib_name} + PRIVATE + + Drone.cpp + drone_ops.cpp + ) diff --git a/src/wind_energy/actuator/drone/Drone.H b/src/wind_energy/actuator/drone/Drone.H new file mode 100644 index 0000000000..329f240ebd --- /dev/null +++ b/src/wind_energy/actuator/drone/Drone.H @@ -0,0 +1,65 @@ +#ifndef DRONE_H +#define DRONE_H + +#include "src/wind_energy/actuator/actuator_types.H" +#include "src/wind_energy/actuator/sector/ActuatorSector.H" +#include "src/wind_energy/actuator/sector/actuator_sector_ops.H" + +#include +#include + +namespace kynema_sgf::actuator { + +struct DroneRotor +{ + using SectorData = ActuatorSector::DataType; + + DroneRotor(CFDSim& sim, const std::string& label, int id); + + SectorData data; + ops::ActSrcOp source; + ops::ProcessOutputsOp output; + vs::Vector body_offset{vs::Vector::zero()}; +}; + +struct DroneData +{ + int num_rotors{0}; + amrex::Real arm_length{0.0_rt}; + RealList arm_lengths; + RealList arm_angles_degrees; + amrex::Real arm_phase_degrees{0.0_rt}; + vs::Vector center{vs::Vector::zero()}; + vs::Vector translation_velocity{vs::Vector::zero()}; + vs::Vector body_orientation_degrees{vs::Vector::zero()}; + vs::Tensor body_orientation{vs::Tensor::identity()}; + RealList rotor_omegas; + amrex::Vector mirror_blades; + vs::Vector total_force{vs::Vector::zero()}; + vs::Vector total_moment{vs::Vector::zero()}; + amrex::Vector> rotors; +}; + +struct Drone : public ActuatorType +{ + using InfoType = ActInfo; + using GridType = ActGrid; + using MetaType = DroneData; + using DataType = ActDataHolder; + + static std::string identifier() { return "Drone"; } +}; + +namespace drone { + +vs::Tensor body_rotation(const vs::Vector& roll_pitch_yaw_degrees); + +RealList uniform_arm_angles(int num_rotors, amrex::Real phase_degrees); + +VecList +rotor_body_offsets(const RealList& arm_lengths, const RealList& angles_degrees); + +} // namespace drone +} // namespace kynema_sgf::actuator + +#endif diff --git a/src/wind_energy/actuator/drone/Drone.cpp b/src/wind_energy/actuator/drone/Drone.cpp new file mode 100644 index 0000000000..09f8d935c1 --- /dev/null +++ b/src/wind_energy/actuator/drone/Drone.cpp @@ -0,0 +1,49 @@ +#include "src/wind_energy/actuator/drone/Drone.H" +#include "src/wind_energy/actuator/drone/drone_ops.H" +#include "src/wind_energy/actuator/ActuatorModel.H" + +#include +#include + +namespace kynema_sgf::actuator { + +DroneRotor::DroneRotor(CFDSim& sim, const std::string& label, const int id) + : data(sim, label, id), source(data), output(data) +{} + +namespace drone { + +vs::Tensor body_rotation(const vs::Vector& angles) +{ + return vs::zrot(angles.z()) & vs::yrot(angles.y()) & vs::xrot(angles.x()); +} + +RealList uniform_arm_angles(const int num_rotors, const amrex::Real phase) +{ + RealList angles(num_rotors); + for (int i = 0; i < num_rotors; ++i) { + angles[i] = phase + 360.0_rt * static_cast(i) / + static_cast(num_rotors); + } + return angles; +} + +VecList +rotor_body_offsets(const RealList& lengths, const RealList& angles_degrees) +{ + AMREX_ALWAYS_ASSERT(lengths.size() == angles_degrees.size()); + VecList offsets(lengths.size()); + for (int i = 0; i < static_cast(lengths.size()); ++i) { + const amrex::Real angle = + angles_degrees[i] * std::numbers::pi_v / 180.0_rt; + offsets[i] = { + lengths[i] * std::cos(angle), lengths[i] * std::sin(angle), 0.0_rt}; + } + return offsets; +} + +} // namespace drone + +template class ActModel; + +} // namespace kynema_sgf::actuator diff --git a/src/wind_energy/actuator/drone/drone_ops.H b/src/wind_energy/actuator/drone/drone_ops.H new file mode 100644 index 0000000000..2620515673 --- /dev/null +++ b/src/wind_energy/actuator/drone/drone_ops.H @@ -0,0 +1,71 @@ +#ifndef DRONE_OPS_H +#define DRONE_OPS_H + +#include "src/wind_energy/actuator/drone/Drone.H" +#include "src/wind_energy/actuator/actuator_ops.H" + +namespace kynema_sgf::actuator::ops { + +template <> +struct ReadInputsOp +{ + void operator()(Drone::DataType&, const utils::ActParser&); +}; + +template <> +struct InitDataOp +{ + void operator()(Drone::DataType&); +}; + +template <> +struct UpdatePosOp +{ + void operator()(Drone::DataType&); +}; + +template <> +struct UpdateVelOp +{ + void operator()(Drone::DataType&); +}; + +template <> +struct ComputeForceOp +{ + void operator()(Drone::DataType&); +}; + +template <> +struct ProcessOutputsOp +{ + explicit ProcessOutputsOp(Drone::DataType& data) : m_data(data) {} + void read_io_options(const utils::ActParser& pp) + { + pp.query("output_frequency", m_output_frequency); + } + void prepare_outputs(const std::string& out_dir); + void write_outputs(); + +private: + Drone::DataType& m_data; + std::string m_filename; + int m_output_frequency{1}; +}; + +template <> +class ActSrcOp +{ +public: + explicit ActSrcOp(Drone::DataType& data) : m_data(data) {} + void initialize(); + void setup_op(); + void operator()(int, const amrex::MFIter&, const amrex::Geometry&); + +private: + Drone::DataType& m_data; +}; + +} // namespace kynema_sgf::actuator::ops + +#endif diff --git a/src/wind_energy/actuator/drone/drone_ops.cpp b/src/wind_energy/actuator/drone/drone_ops.cpp new file mode 100644 index 0000000000..6ff6478923 --- /dev/null +++ b/src/wind_energy/actuator/drone/drone_ops.cpp @@ -0,0 +1,293 @@ +#include "src/wind_energy/actuator/drone/drone_ops.H" + +#include "src/CFDSim.H" +#include "src/utilities/constants.H" + +#include +#include +#include +#include + +#include "AMReX_Print.H" + +using namespace amrex::literals; + +namespace kynema_sgf::actuator::ops { +namespace { + +void input_error(const std::string& label, const std::string& message) +{ + amrex::Abort("Drone '" + label + "': " + message); +} + +void validate_distinct_angles(const std::string& label, const RealList& angles) +{ + for (int i = 0; i < static_cast(angles.size()); ++i) { + for (int j = i + 1; j < static_cast(angles.size()); ++j) { + amrex::Real delta = std::remainder(angles[i] - angles[j], 360.0_rt); + if (std::abs(delta) < constants::TIGHT_TOL) { + input_error(label, "arm angles must be distinct"); + } + } + } +} + +} // namespace + +void ReadInputsOp::operator()( + Drone::DataType& data, const utils::ActParser& pp) +{ + auto& meta = data.meta(); + const auto& label = data.info().label; + + pp.get("num_rotors", meta.num_rotors); + if (meta.num_rotors <= 0) { + input_error(label, "num_rotors must be positive"); + } + pp.get("center", meta.center); + pp.query("translation_velocity", meta.translation_velocity); + pp.query("body_orientation_degrees", meta.body_orientation_degrees); + meta.body_orientation = drone::body_rotation(meta.body_orientation_degrees); + + const auto& instance = pp.params(); + const auto& defaults = pp.default_params(); + const bool instance_length = instance.contains("arm_length"); + const bool instance_lengths = instance.contains("arm_lengths"); + if (instance_length && instance_lengths) { + input_error(label, "specify arm_length or arm_lengths, not both"); + } + if (instance_lengths) { + pp.getarr("arm_lengths", meta.arm_lengths); + } else if (instance_length) { + pp.get("arm_length", meta.arm_length); + meta.arm_lengths.assign(meta.num_rotors, meta.arm_length); + } else if ( + defaults.contains("arm_lengths") && defaults.contains("arm_length")) { + input_error(label, "specify arm_length or arm_lengths, not both"); + } else if (defaults.contains("arm_lengths")) { + pp.getarr("arm_lengths", meta.arm_lengths); + } else if (defaults.contains("arm_length")) { + pp.get("arm_length", meta.arm_length); + meta.arm_lengths.assign(meta.num_rotors, meta.arm_length); + } else { + input_error(label, "arm_length or arm_lengths is required"); + } + if (static_cast(meta.arm_lengths.size()) != meta.num_rotors) { + input_error(label, "arm_lengths must contain num_rotors values"); + } + if (std::ranges::any_of(meta.arm_lengths, [](const amrex::Real length) { + return length <= 0.0_rt; + })) { + input_error(label, "all arm lengths must be positive"); + } + + const bool instance_angles = instance.contains("arm_angles_degrees"); + const bool instance_phase = instance.contains("arm_phase_degrees"); + if (instance_angles && instance_phase) { + input_error( + label, + "arm_angles_degrees and arm_phase_degrees are mutually exclusive"); + } + if (instance_angles) { + pp.getarr("arm_angles_degrees", meta.arm_angles_degrees); + } else if (instance_phase) { + pp.get("arm_phase_degrees", meta.arm_phase_degrees); + meta.arm_angles_degrees = + drone::uniform_arm_angles(meta.num_rotors, meta.arm_phase_degrees); + } else if ( + defaults.contains("arm_angles_degrees") && + defaults.contains("arm_phase_degrees")) { + input_error( + label, + "arm_angles_degrees and arm_phase_degrees are mutually exclusive"); + } else if (defaults.contains("arm_angles_degrees")) { + pp.getarr("arm_angles_degrees", meta.arm_angles_degrees); + } else { + pp.query("arm_phase_degrees", meta.arm_phase_degrees); + meta.arm_angles_degrees = + drone::uniform_arm_angles(meta.num_rotors, meta.arm_phase_degrees); + } + if (static_cast(meta.arm_angles_degrees.size()) != meta.num_rotors) { + input_error(label, "arm_angles_degrees must contain num_rotors values"); + } + validate_distinct_angles(label, meta.arm_angles_degrees); + + pp.getarr("rotor_omegas", meta.rotor_omegas); + if (static_cast(meta.rotor_omegas.size()) != meta.num_rotors) { + input_error(label, "rotor_omegas must contain num_rotors values"); + } + + meta.mirror_blades.assign(meta.num_rotors, false); + if (pp.contains("mirror_blades")) { + amrex::Vector mirror_inputs; + pp.getarr("mirror_blades", mirror_inputs); + if (static_cast(mirror_inputs.size()) != meta.num_rotors) { + input_error(label, "mirror_blades must contain num_rotors values"); + } + for (int i = 0; i < meta.num_rotors; ++i) { + const auto value = amrex::toLower(mirror_inputs[i]); + if ((value != "true") && (value != "false")) { + input_error(label, "mirror_blades values must be true or false"); + } + meta.mirror_blades[i] = (value == "true"); + } + } + + const auto offsets = + drone::rotor_body_offsets(meta.arm_lengths, meta.arm_angles_degrees); + const auto rotor_normal = meta.body_orientation & vs::Vector::khat(); + meta.rotors.reserve(meta.num_rotors); + for (int i = 0; i < meta.num_rotors; ++i) { + const std::string rotor_label = label + ".R" + std::to_string(i + 1); + auto rotor = std::make_unique( + data.sim(), rotor_label, data.info().id * 100000 + i); + rotor->body_offset = offsets[i]; + + // A drone owns its rotor configuration. Reuse the same default and + // instance namespaces for both drone geometry and sector aerodynamics; + // each reader simply ignores inputs outside its responsibility. + utils::ActParser sector_pp("Actuator.Drone", "Actuator." + label); + rotor->data.meta().omega = meta.rotor_omegas[i]; + rotor->data.meta().user_omega = true; + ReadInputsOp()(rotor->data, sector_pp); + if (meta.mirror_blades[i]) { + for (auto& twist : rotor->data.meta().twist_inp) { + twist = -twist; + } + } + const auto rotor_center = + meta.center + (meta.body_orientation & rotor->body_offset); + sector::set_placement( + rotor->data, rotor_center, rotor_normal, meta.translation_velocity); + rotor->output.read_io_options(sector_pp); + meta.rotors.emplace_back(std::move(rotor)); + } + + amrex::GpuArray lo; + amrex::GpuArray hi; + for (int n = 0; n < AMREX_SPACEDIM; ++n) { + lo[n] = std::numeric_limits::max(); + hi[n] = std::numeric_limits::lowest(); + } + for (const auto& rotor : meta.rotors) { + const auto& box = rotor->data.info().bound_box; + for (int n = 0; n < AMREX_SPACEDIM; ++n) { + lo[n] = std::min(lo[n], box.lo(n)); + hi[n] = std::max(hi[n], box.hi(n)); + } + } + data.info().bound_box = amrex::RealBox(lo.data(), hi.data()); +} + +void InitDataOp::operator()(Drone::DataType& data) +{ + auto& meta = data.meta(); + int total_points = 0; + for (auto& rotor : meta.rotors) { + InitDataOp()(rotor->data); + total_points += static_cast(rotor->data.grid().vel_pos.size()); + } + data.grid().resize(0, total_points); + UpdatePosOp()(data); +} + +void UpdatePosOp::operator()(Drone::DataType& data) +{ + auto& grid = data.grid(); + int offset = 0; + for (auto& rotor : data.meta().rotors) { + UpdatePosOp()(rotor->data); + const auto& positions = rotor->data.grid().vel_pos; + std::copy( + positions.begin(), positions.end(), grid.vel_pos.begin() + offset); + offset += static_cast(positions.size()); + } +} + +void UpdateVelOp::operator()(Drone::DataType& data) +{ + auto& grid = data.grid(); + int offset = 0; + for (auto& rotor : data.meta().rotors) { + auto& rotor_grid = rotor->data.grid(); + const int npts = static_cast(rotor_grid.vel.size()); + std::copy_n(grid.vel.begin() + offset, npts, rotor_grid.vel.begin()); + std::copy_n( + grid.density.begin() + offset, npts, rotor_grid.density.begin()); + offset += npts; + } +} + +void ComputeForceOp::operator()(Drone::DataType& data) +{ + auto& meta = data.meta(); + const auto& time = data.sim().time(); + const amrex::Real midpoint_time = + time.current_time() + 0.5_rt * sector::timestep_width(data.sim()); + const auto drone_center = + meta.center + meta.translation_velocity * midpoint_time; + meta.total_force = vs::Vector::zero(); + meta.total_moment = vs::Vector::zero(); + for (auto& rotor : meta.rotors) { + ComputeForceOp()(rotor->data); + const auto& rotor_meta = rotor->data.meta(); + // Sector forces act on the fluid; report equal-and-opposite vehicle + // load. + meta.total_force = meta.total_force - rotor_meta.integrated_force; + meta.total_moment = + meta.total_moment - rotor_meta.integrated_moment - + ((rotor_meta.center - drone_center) ^ rotor_meta.integrated_force); + } +} + +void ProcessOutputsOp::prepare_outputs( + const std::string& out_dir) +{ + m_filename = out_dir + "/" + m_data.info().label + "_loads.csv"; + std::ofstream stream(m_filename); + stream << "time,force_x,force_y,force_z,moment_x,moment_y,moment_z\n"; + for (auto& rotor : m_data.meta().rotors) { + rotor->output.prepare_outputs(out_dir); + } +} + +void ProcessOutputsOp::write_outputs() +{ + const auto& time = m_data.sim().time(); + if ((m_output_frequency > 0) && + (time.time_index() % m_output_frequency == 0)) { + const auto& force = m_data.meta().total_force; + const auto& moment = m_data.meta().total_moment; + std::ofstream stream(m_filename, std::ios::app); + stream << time.new_time() << ',' << force.x() << ',' << force.y() << ',' + << force.z() << ',' << moment.x() << ',' << moment.y() << ',' + << moment.z() << '\n'; + } + for (auto& rotor : m_data.meta().rotors) { + rotor->output.write_outputs(); + } +} + +void ActSrcOp::initialize() +{ + for (auto& rotor : m_data.meta().rotors) { + rotor->source.initialize(); + } +} + +void ActSrcOp::setup_op() +{ + for (auto& rotor : m_data.meta().rotors) { + rotor->source.setup_op(); + } +} + +void ActSrcOp::operator()( + const int lev, const amrex::MFIter& mfi, const amrex::Geometry& geom) +{ + for (auto& rotor : m_data.meta().rotors) { + rotor->source(lev, mfi, geom); + } +} + +} // namespace kynema_sgf::actuator::ops diff --git a/src/wind_energy/actuator/sector/ActuatorSector.H b/src/wind_energy/actuator/sector/ActuatorSector.H index 0f3228c8bc..fa07028c7a 100644 --- a/src/wind_energy/actuator/sector/ActuatorSector.H +++ b/src/wind_energy/actuator/sector/ActuatorSector.H @@ -43,6 +43,7 @@ struct ActuatorSectorData //! Rotor spin angular speed [rad/s] amrex::Real omega{0.0_rt}; + bool user_omega{false}; //! Initial rotor hub center in the CFD frame [m] vs::Vector center0{0.0_rt, 0.0_rt, 0.0_rt}; @@ -178,6 +179,14 @@ struct ActuatorSectorData //! Integrated torque about the rotor normal [m^5/s^2] amrex::Real torque{0.0_rt}; + + //! Integrated density-normalized force applied to the fluid in the CFD + //! frame [m^4/s^2] + vs::Vector integrated_force{0.0_rt, 0.0_rt, 0.0_rt}; + + //! Integrated density-normalized moment applied to the fluid about the + //! rotor hub [m^5/s^2] + vs::Vector integrated_moment{0.0_rt, 0.0_rt, 0.0_rt}; }; struct ActuatorSector : public ActuatorType diff --git a/src/wind_energy/actuator/sector/actuator_sector_ops.H b/src/wind_energy/actuator/sector/actuator_sector_ops.H index 544c5515bd..8aa0c8fa00 100644 --- a/src/wind_energy/actuator/sector/actuator_sector_ops.H +++ b/src/wind_energy/actuator/sector/actuator_sector_ops.H @@ -43,6 +43,13 @@ void blade_basis( void update_midpoint_sample_points(ActuatorSector::DataType& data); +//! Set the initial placement and prescribed translation of a mounted rotor. +void set_placement( + ActuatorSector::DataType& data, + const vs::Vector& center, + const vs::Vector& rotor_normal, + const vs::Vector& translation_velocity); + void build_gaussian_table(ActuatorSectorData& meta); amrex::Real timestep_width(const CFDSim& sim); diff --git a/src/wind_energy/actuator/sector/actuator_sector_ops.cpp b/src/wind_energy/actuator/sector/actuator_sector_ops.cpp index 3eca97ffc9..1be1e8e537 100644 --- a/src/wind_energy/actuator/sector/actuator_sector_ops.cpp +++ b/src/wind_energy/actuator/sector/actuator_sector_ops.cpp @@ -255,6 +255,38 @@ void update_midpoint_sample_points(ActuatorSector::DataType& data) } } +void set_placement( + ActuatorSector::DataType& data, + const vs::Vector& center, + const vs::Vector& rotor_normal, + const vs::Vector& translation_velocity) +{ + auto& meta = data.meta(); + meta.center0 = center; + meta.center = center; + meta.rotor_normal = rotor_normal; + meta.translation_velocity = translation_velocity; + meta.rotor_angular_velocity = vs::Vector::zero(); + meta.user_rotor_angular_velocity = true; + + const amrex::Real max_eps = local_epsilon(meta, max_interp_chord(meta)); + const amrex::Real search_radius = + meta.rotor_radius + meta.support_radius_over_epsilon * max_eps; + if (vs::mag(translation_velocity) > + std::numeric_limits::epsilon()) { + const auto& geom = data.sim().mesh().Geom(0); + const auto plo = geom.ProbLoArray(); + const auto phi = geom.ProbHiArray(); + data.info().bound_box = + amrex::RealBox(plo[0], plo[1], plo[2], phi[0], phi[1], phi[2]); + } else { + data.info().bound_box = amrex::RealBox( + center.x() - search_radius, center.y() - search_radius, + center.z() - search_radius, center.x() + search_radius, + center.y() + search_radius, center.z() + search_radius); + } +} + void build_gaussian_table(ActuatorSectorData& meta) { if (meta.gaussian_table_error <= 0.0_rt) { @@ -425,7 +457,15 @@ void ReadInputsOp::operator()( { auto& meta = data.meta(); pp.get("rotor_diameter", meta.rotor_diameter); - pp.get("omega", meta.omega); + if (meta.user_omega) { + if (pp.contains("omega")) { + amrex::Abort( + "ActuatorSector omega was set by its owning model and must " + "not also be specified in the shared sector inputs"); + } + } else { + pp.get("omega", meta.omega); + } pp.get("airfoil_table", meta.airfoil_file); pp.query("airfoil_type", meta.airfoil_type); pp.query("num_blades", meta.num_blades); @@ -602,6 +642,8 @@ void ComputeForceOp::operator()( VecList force_theta_normal(nvel, vs::Vector::zero()); meta.thrust = 0.0_rt; meta.torque = 0.0_rt; + meta.integrated_force = vs::Vector::zero(); + meta.integrated_moment = vs::Vector::zero(); // First compute blade-section loads at the midpoint sample locations. The // force here is per unit span [m^3/s^2] in kinematic units and is converted @@ -720,6 +762,10 @@ void ComputeForceOp::operator()( meta.dr[ir] / static_cast(ntheta); grid.force[iq] = wt * ((e_theta * force_theta_normal[ip].x()) + (e_normal * force_theta_normal[ip].y())); + meta.integrated_force = meta.integrated_force + grid.force[iq]; + meta.integrated_moment = + meta.integrated_moment + + ((grid.pos[iq] - meta.center) ^ grid.force[iq]); grid.epsilon[iq] = vs::Vector::one() * meta.epsilon_profile[ir]; ++iq; } diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index c74bafe687..f5a2dfe774 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -247,6 +247,7 @@ add_test_r(act_moving_wing) add_test_r(act_pitching_wing_2D) add_test_r(act_sector) add_test_r(act_sector_quad) +add_test_r(act_drone_quad) add_test_r(box_refinement) add_test_r(burggraf_flow) add_test_r(channel_builder_multiphase_contraction) diff --git a/test/test_files/act_drone_quad/act_drone_quad.inp b/test/test_files/act_drone_quad/act_drone_quad.inp new file mode 100644 index 0000000000..3a1d735932 --- /dev/null +++ b/test/test_files/act_drone_quad/act_drone_quad.inp @@ -0,0 +1,83 @@ +# Time controls +time.max_step = 4 +time.fixed_dt = 1.7453292519943296e-4 # 25 deg rotor sweep per step [s] +time.use_force_cfl = false +time.plot_interval = -1 +time.checkpoint_interval = -1 + +# Domain extents are Cartesian coordinates [m]. +geometry.prob_lo = -0.1 -0.2 -0.2 +geometry.prob_hi = 0.2 0.2 0.2 +geometry.is_periodic = 0 0 0 +amr.n_cell = 48 64 64 +amr.max_level = 0 +amr.max_grid_size = 16 + +# Constant initial/background fields: density [kg/m^3], velocity [m/s]. +ConstValue.density.value = 1.0 +ConstValue.velocity.value = 1.0 0.0 0.0 + +incflo.use_godunov = 1 +incflo.diffusion_type = 2 +incflo.do_initial_proj = 1 +incflo.initial_iterations = 3 +transport.viscosity = 1.0e-5 +transport.laminar_prandtl = 0.7 +transport.turbulent_prandtl = 0.3333 +turbulence.model = Laminar + +incflo.physics = FreeStream Actuator +ICNS.source_terms = ActuatorForcing + +Actuator.labels = D1 +Actuator.D1.type = Drone + +# Shared drone layout and actuator-sector rotor definition. The arm length is +# the body-center to rotor-hub distance. A 45-degree phase creates an X layout. +Actuator.Drone.num_rotors = 4 +Actuator.Drone.arm_length = 0.10606601717798213 +Actuator.Drone.arm_phase_degrees = 45.0 +Actuator.Drone.rotor_omegas = 2500.0 -2500.0 2500.0 -2500.0 +Actuator.Drone.mirror_blades = false true false true + +Actuator.Drone.rotor_diameter = 0.10 +Actuator.Drone.root_radius_fraction = 0.18 +Actuator.Drone.num_blades = 2 + +# Rotate the default body xy rotor plane into the CFD yz plane. This maps the +# body +z rotor normal to CFD +x and reproduces act_sector_quad placement. +Actuator.D1.center = 0.0 0.0 0.0 +Actuator.D1.body_orientation_degrees = -90.0 0.0 -90.0 + +# Actuator-sector Gaussian width and generated sector grid +Actuator.Drone.epsilon_chord = 0.5 +Actuator.Drone.epsilon_dr = 1.5 +Actuator.Drone.epsilon_dl = 1.5 +Actuator.Drone.min_chord_dr = 2.0 + +# Actuator-sector blade section inputs +Actuator.Drone.span_locs = 0.0 1.0 +Actuator.Drone.chord = 0.01 0.006 +Actuator.Drone.twist = 12.0 4.0 +Actuator.Drone.airfoil_table = ../act_fixed_wing/DU21_A17.txt +Actuator.Drone.airfoil_type = openfast + +# Actuator-sector force projection and output +Actuator.Drone.gaussian_type = table +Actuator.Drone.gaussian_table_error = 1.0e-4 +Actuator.Drone.output_frequency = 1 + +# Boundary conditions: +x inflow, +x pressure outflow, transverse slip walls. +xlo.type = "mass_inflow" +xlo.density = 1.0 +xlo.velocity = 1.0 0.0 0.0 +xhi.type = "pressure_outflow" + +ylo.type = "slip_wall" +yhi.type = "slip_wall" + +zlo.type = "slip_wall" +zhi.type = "slip_wall" + +incflo.verbose = 0 +nodal_proj.verbose = 0 diff --git a/unit_tests/wind_energy/actuator/CMakeLists.txt b/unit_tests/wind_energy/actuator/CMakeLists.txt index baa158acc0..3c954ca640 100644 --- a/unit_tests/wind_energy/actuator/CMakeLists.txt +++ b/unit_tests/wind_energy/actuator/CMakeLists.txt @@ -5,6 +5,7 @@ target_sources(${kynema_sgf_unit_test_exe_name} PRIVATE test_airfoil.cpp test_actuator_free_functions.cpp test_actuator_sector.cpp + test_drone.cpp test_disk_uniform_ct.cpp test_FLLC.cpp test_actuator_joukowsky_disk.cpp diff --git a/unit_tests/wind_energy/actuator/test_drone.cpp b/unit_tests/wind_energy/actuator/test_drone.cpp new file mode 100644 index 0000000000..5b9b0d4d80 --- /dev/null +++ b/unit_tests/wind_energy/actuator/test_drone.cpp @@ -0,0 +1,327 @@ +#include "src/wind_energy/actuator/drone/Drone.H" +#include "src/wind_energy/actuator/drone/drone_ops.H" +#include "src/wind_energy/actuator/Actuator.H" +#include "src/wind_energy/actuator/ActuatorContainer.H" +#include "src/wind_energy/actuator/ActuatorModel.H" +#include "src/utilities/constants.H" +#include "ks_test_utils/MeshTest.H" + +#include "gtest/gtest.h" + +#include +#include + +#include "AMReX_REAL.H" + +using namespace amrex::literals; + +namespace kynema_sgf_tests { +namespace { + +using kynema_sgf::actuator::drone::rotor_body_offsets; +using kynema_sgf::actuator::drone::uniform_arm_angles; + +class DroneActuatorTest : public MeshTest +{ +protected: + void populate_parameters() override + { + MeshTest::populate_parameters(); + amrex::ParmParse pp_amr("amr"); + pp_amr.add("max_level", 0); + pp_amr.add("max_grid_size", 16); + pp_amr.addarr("n_cell", amrex::Vector{16, 16, 16}); + + amrex::ParmParse pp_time("time"); + pp_time.add("fixed_dt", 1.0e-4_rt); + + amrex::ParmParse pp_geom("geometry"); + pp_geom.addarr( + "prob_lo", amrex::Vector{-0.2_rt, -0.2_rt, -0.2_rt}); + pp_geom.addarr( + "prob_hi", amrex::Vector{0.2_rt, 0.2_rt, 0.2_rt}); + } + + void initialize_domain() + { + initialize_mesh(); + sim().repo().declare_field("actuator_src_term", 3, 0); + auto& vel = sim().repo().declare_field("velocity", 3, 3); + auto& density = sim().repo().declare_field("density", 1, 3); + vel.setVal(0.0_rt); + density.setVal(1.0_rt); + kynema_sgf::actuator::ActuatorContainer::ParticleType::NextID(1U); + } + + void populate_inputs() + { + amrex::ParmParse pp_a("Actuator"); + pp_a.add("labels", std::string("D1")); + + amrex::ParmParse pp_d("Actuator.Drone"); + pp_d.add("num_rotors", 4); + pp_d.add("arm_length", 0.075_rt); + pp_d.add("arm_phase_degrees", 45.0_rt); + pp_d.add("rotor_diameter", 0.05_rt); + pp_d.add("num_blades", 2); + pp_d.add("root_radius_fraction", 0.18_rt); + pp_d.add("epsilon_chord", 0.5_rt); + pp_d.add("airfoil_table", m_airfoil_file); + pp_d.add("airfoil_type", std::string("openfast")); + pp_d.addarr("span_locs", amrex::Vector{0.0_rt, 1.0_rt}); + pp_d.addarr("chord", amrex::Vector{0.01_rt, 0.006_rt}); + pp_d.addarr("twist", amrex::Vector{12.0_rt, 4.0_rt}); + + amrex::ParmParse pp_i("Actuator.D1"); + pp_i.add("type", std::string("Drone")); + pp_i.addarr( + "arm_lengths", + amrex::Vector{0.075_rt, 0.075_rt, 0.075_rt, 0.075_rt}); + pp_i.addarr( + "arm_angles_degrees", + amrex::Vector{45.0_rt, 135.0_rt, 225.0_rt, 315.0_rt}); + pp_i.addarr( + "center", amrex::Vector{0.0_rt, 0.0_rt, 0.0_rt}); + pp_i.addarr( + "translation_velocity", + amrex::Vector{0.1_rt, -0.2_rt, 0.3_rt}); + pp_i.addarr( + "rotor_omegas", + amrex::Vector{ + 2500.0_rt, -2500.0_rt, 2500.0_rt, -2500.0_rt}); + pp_i.addarr( + "mirror_blades", + amrex::Vector{"false", "true", "false", "true"}); + + amrex::ParmParse pp_s("Actuator.ActuatorSector"); + pp_s.add("rotor_diameter", 0.05_rt); + pp_s.add("num_blades", 2); + pp_s.add("root_radius_fraction", 0.18_rt); + pp_s.add("epsilon_chord", 0.5_rt); + pp_s.add("airfoil_table", m_airfoil_file); + pp_s.add("airfoil_type", std::string("openfast")); + pp_s.addarr("span_locs", amrex::Vector{0.0_rt, 1.0_rt}); + pp_s.addarr("chord", amrex::Vector{0.01_rt, 0.006_rt}); + pp_s.addarr("twist", amrex::Vector{12.0_rt, 4.0_rt}); + + const amrex::Real offset = 0.075_rt / std::sqrt(2.0_rt); + const amrex::Vector> centers{ + {offset, offset, 0.0_rt}, + {-offset, offset, 0.0_rt}, + {-offset, -offset, 0.0_rt}, + {offset, -offset, 0.0_rt}}; + for (int i = 0; i < 4; ++i) { + amrex::ParmParse pp_r("Actuator.R" + std::to_string(i + 1)); + pp_r.add("omega", (i % 2 == 0) ? 2500.0_rt : -2500.0_rt); + if (i % 2 != 0) { + pp_r.addarr( + "twist", + amrex::Vector{-12.0_rt, -4.0_rt}); + } + pp_r.addarr("center", centers[i]); + pp_r.addarr( + "translation_velocity", + amrex::Vector{0.1_rt, -0.2_rt, 0.3_rt}); + pp_r.addarr( + "rotor_normal", + amrex::Vector{0.0_rt, 0.0_rt, 1.0_rt}); + } + } + + void write_airfoil() const + { + std::ofstream os(m_airfoil_file); + os << "! test polar\n5 NumAlf\n! Alpha Cl Cd Cm\n! deg - - -\n" + "-180 0 0.04 0\n-10 -0.5 0.02 0\n0 0 0.01 0\n" + "10 0.5 0.02 0\n180 0 0.04 0\n"; + } + + const std::string m_airfoil_file{"drone_airfoil.txt"}; +}; + +class DronePhysicsTest : public kynema_sgf::actuator::Actuator +{ +public: + explicit DronePhysicsTest(kynema_sgf::CFDSim& sim) : Actuator(sim) {} + +protected: + void prepare_outputs() override {} +}; + +TEST(DroneGeometry, plus_layout) +{ + const auto angles = uniform_arm_angles(4, 0.0_rt); + const auto offsets = + rotor_body_offsets({2.0_rt, 2.0_rt, 2.0_rt, 2.0_rt}, angles); + + ASSERT_EQ(offsets.size(), 4U); + EXPECT_NEAR(offsets[0].x(), 2.0_rt, 1.0e-14_rt); + EXPECT_NEAR(offsets[0].y(), 0.0_rt, 1.0e-14_rt); + EXPECT_NEAR(offsets[1].x(), 0.0_rt, 1.0e-14_rt); + EXPECT_NEAR(offsets[1].y(), 2.0_rt, 1.0e-14_rt); + EXPECT_NEAR(offsets[2].x(), -2.0_rt, 1.0e-14_rt); + EXPECT_NEAR(offsets[3].y(), -2.0_rt, 1.0e-14_rt); +} + +TEST(DroneGeometry, unequal_irregular_arms) +{ + const kynema_sgf::actuator::RealList lengths{1.0_rt, 2.0_rt, 3.0_rt}; + const kynema_sgf::actuator::RealList angles{0.0_rt, 90.0_rt, 225.0_rt}; + const auto offsets = rotor_body_offsets(lengths, angles); + + EXPECT_NEAR(offsets[0].x(), 1.0_rt, 1.0e-14_rt); + EXPECT_NEAR(offsets[1].y(), 2.0_rt, 1.0e-14_rt); + EXPECT_NEAR(offsets[2].x(), -3.0_rt / std::sqrt(2.0_rt), 1.0e-14_rt); + EXPECT_NEAR(offsets[2].y(), -3.0_rt / std::sqrt(2.0_rt), 1.0e-14_rt); +} + +TEST(DroneGeometry, body_orientation_defaults_to_identity) +{ + const auto rotation = kynema_sgf::actuator::drone::body_rotation( + kynema_sgf::vs::Vector::zero()); + const auto value = + rotation & kynema_sgf::vs::Vector{1.0_rt, 2.0_rt, 3.0_rt}; + EXPECT_NEAR(value.x(), 1.0_rt, 1.0e-14_rt); + EXPECT_NEAR(value.y(), 2.0_rt, 1.0e-14_rt); + EXPECT_NEAR(value.z(), 3.0_rt, 1.0e-14_rt); +} + +TEST_F(DroneActuatorTest, composite_lifecycle) +{ + write_airfoil(); + initialize_domain(); + populate_inputs(); + + DronePhysicsTest actuator(sim()); + actuator.pre_init_actions(); + auto* drone = dynamic_cast*>( + &actuator.get_act(0)); + ASSERT_NE(drone, nullptr); + ASSERT_EQ(drone->meta().rotors.size(), 4U); + EXPECT_NEAR( + drone->meta().rotors[0]->data.meta().center.x(), + 0.075_rt / std::sqrt(2.0_rt), 1.0e-14_rt); + EXPECT_NEAR( + drone->meta().rotors[0]->data.meta().center.y(), + 0.075_rt / std::sqrt(2.0_rt), 1.0e-14_rt); + + actuator.post_init_actions(); + actuator.pre_advance_work(); + + remove(m_airfoil_file.c_str()); +} + +TEST_F(DroneActuatorTest, matches_equivalent_standalone_sectors) +{ + namespace act = kynema_sgf::actuator; + constexpr amrex::Real tol = kynema_sgf::constants::TIGHT_TOL; + + write_airfoil(); + initialize_domain(); + populate_inputs(); + + act::ActModel drone(sim(), "D1", 0); + drone.read_inputs(act::utils::ActParser("Actuator.Drone", "Actuator.D1")); + drone.init_actuator_source(); + ASSERT_EQ(drone.meta().rotors.size(), 4U); + + for (int i = 0; i < 4; ++i) { + const std::string label = "R" + std::to_string(i + 1); + act::ActModel sector( + sim(), label, i); + sector.read_inputs( + act::utils::ActParser( + "Actuator.ActuatorSector", "Actuator." + label)); + sector.init_actuator_source(); + + const auto& drone_data = drone.meta().rotors[i]->data; + const auto& drone_meta = drone_data.meta(); + const auto& sector_meta = sector.meta(); + const auto& drone_grid = drone_data.grid(); + const auto& sector_grid = sector.grid(); + + EXPECT_EQ(drone_meta.num_blades, sector_meta.num_blades); + EXPECT_NEAR(drone_meta.rotor_diameter, sector_meta.rotor_diameter, tol); + EXPECT_NEAR(drone_meta.rotor_radius, sector_meta.rotor_radius, tol); + EXPECT_NEAR(drone_meta.root_radius, sector_meta.root_radius, tol); + EXPECT_NEAR(drone_meta.omega, sector_meta.omega, tol); + + for (int n = 0; n < AMREX_SPACEDIM; ++n) { + EXPECT_NEAR(drone_meta.center[n], sector_meta.center[n], tol); + EXPECT_NEAR( + drone_meta.translation_velocity[n], + sector_meta.translation_velocity[n], tol); + EXPECT_NEAR( + drone_meta.rotor_normal[n], sector_meta.rotor_normal[n], tol); + } + + ASSERT_EQ(drone_meta.radius.size(), sector_meta.radius.size()); + EXPECT_EQ(drone_meta.dr.size(), sector_meta.dr.size()); + EXPECT_EQ(drone_meta.chord.size(), sector_meta.chord.size()); + EXPECT_EQ(drone_meta.twist.size(), sector_meta.twist.size()); + EXPECT_EQ( + drone_meta.epsilon_profile.size(), + sector_meta.epsilon_profile.size()); + EXPECT_EQ(drone_meta.theta_counts, sector_meta.theta_counts); + for (int j = 0; j < static_cast(drone_meta.radius.size()); ++j) { + EXPECT_NEAR(drone_meta.radius[j], sector_meta.radius[j], tol); + EXPECT_NEAR(drone_meta.dr[j], sector_meta.dr[j], tol); + EXPECT_NEAR(drone_meta.chord[j], sector_meta.chord[j], tol); + EXPECT_NEAR(drone_meta.twist[j], sector_meta.twist[j], tol); + EXPECT_NEAR( + drone_meta.epsilon_profile[j], sector_meta.epsilon_profile[j], + tol); + } + + EXPECT_EQ(drone_grid.pos.size(), sector_grid.pos.size()); + EXPECT_EQ(drone_grid.force.size(), sector_grid.force.size()); + EXPECT_EQ(drone_grid.epsilon.size(), sector_grid.epsilon.size()); + ASSERT_EQ(drone_grid.vel_pos.size(), sector_grid.vel_pos.size()); + EXPECT_EQ(drone_grid.vel.size(), sector_grid.vel.size()); + EXPECT_EQ(drone_grid.density.size(), sector_grid.density.size()); + for (int j = 0; j < static_cast(drone_grid.vel_pos.size()); ++j) { + for (int n = 0; n < AMREX_SPACEDIM; ++n) { + EXPECT_NEAR( + drone_grid.vel_pos[j][n], sector_grid.vel_pos[j][n], tol); + } + } + } + + remove(m_airfoil_file.c_str()); +} + +TEST_F(DroneActuatorTest, counter_rotating_rotors_have_same_axial_force_direction) +{ + namespace act = kynema_sgf::actuator; + + write_airfoil(); + initialize_domain(); + populate_inputs(); + + act::ActModel drone(sim(), "D1", 0); + drone.read_inputs(act::utils::ActParser("Actuator.Drone", "Actuator.D1")); + drone.init_actuator_source(); + ASSERT_EQ(drone.meta().rotors.size(), 4U); + + amrex::Real reference_force_normal = 0.0_rt; + for (int i = 0; i < 4; ++i) { + auto& rotor_data = drone.meta().rotors[i]->data; + act::ops::ComputeForceOp()( + rotor_data); + const auto& rotor_meta = rotor_data.meta(); + const amrex::Real force_normal = + rotor_meta.integrated_force & rotor_meta.rotor_normal; + ASSERT_NE(force_normal, 0.0_rt); + if (i == 0) { + reference_force_normal = force_normal; + } else { + EXPECT_GT(force_normal * reference_force_normal, 0.0_rt); + } + } + + remove(m_airfoil_file.c_str()); +} + +} // namespace +} // namespace kynema_sgf_tests From cd996b732447c55133db5badbaa51926be3299d9 Mon Sep 17 00:00:00 2001 From: Tony Martinez Date: Thu, 16 Jul 2026 08:26:37 -0600 Subject: [PATCH 02/16] Formatting --- src/wind_energy/actuator/drone/drone_ops.cpp | 3 ++- unit_tests/wind_energy/actuator/test_drone.cpp | 11 +++++------ 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/wind_energy/actuator/drone/drone_ops.cpp b/src/wind_energy/actuator/drone/drone_ops.cpp index 6ff6478923..2a61f691a5 100644 --- a/src/wind_energy/actuator/drone/drone_ops.cpp +++ b/src/wind_energy/actuator/drone/drone_ops.cpp @@ -127,7 +127,8 @@ void ReadInputsOp::operator()( for (int i = 0; i < meta.num_rotors; ++i) { const auto value = amrex::toLower(mirror_inputs[i]); if ((value != "true") && (value != "false")) { - input_error(label, "mirror_blades values must be true or false"); + input_error( + label, "mirror_blades values must be true or false"); } meta.mirror_blades[i] = (value == "true"); } diff --git a/unit_tests/wind_energy/actuator/test_drone.cpp b/unit_tests/wind_energy/actuator/test_drone.cpp index 5b9b0d4d80..70a34f145f 100644 --- a/unit_tests/wind_energy/actuator/test_drone.cpp +++ b/unit_tests/wind_energy/actuator/test_drone.cpp @@ -86,9 +86,8 @@ class DroneActuatorTest : public MeshTest "translation_velocity", amrex::Vector{0.1_rt, -0.2_rt, 0.3_rt}); pp_i.addarr( - "rotor_omegas", - amrex::Vector{ - 2500.0_rt, -2500.0_rt, 2500.0_rt, -2500.0_rt}); + "rotor_omegas", amrex::Vector{ + 2500.0_rt, -2500.0_rt, 2500.0_rt, -2500.0_rt}); pp_i.addarr( "mirror_blades", amrex::Vector{"false", "true", "false", "true"}); @@ -115,8 +114,7 @@ class DroneActuatorTest : public MeshTest pp_r.add("omega", (i % 2 == 0) ? 2500.0_rt : -2500.0_rt); if (i % 2 != 0) { pp_r.addarr( - "twist", - amrex::Vector{-12.0_rt, -4.0_rt}); + "twist", amrex::Vector{-12.0_rt, -4.0_rt}); } pp_r.addarr("center", centers[i]); pp_r.addarr( @@ -291,7 +289,8 @@ TEST_F(DroneActuatorTest, matches_equivalent_standalone_sectors) remove(m_airfoil_file.c_str()); } -TEST_F(DroneActuatorTest, counter_rotating_rotors_have_same_axial_force_direction) +TEST_F( + DroneActuatorTest, counter_rotating_rotors_have_same_axial_force_direction) { namespace act = kynema_sgf::actuator; From 9a4b9574fa520042df92da8668487c36ba3d3b94 Mon Sep 17 00:00:00 2001 From: Tony Martinez Date: Fri, 24 Jul 2026 10:17:12 -0600 Subject: [PATCH 03/16] Added drone capability and quaternion formulations --- src/core/vs/quaternion.H | 136 ++++++++++ src/core/vs/tensor.H | 1 + src/core/vs/tensorI.H | 28 -- src/wind_energy/actuator/CMakeLists.txt | 1 + src/wind_energy/actuator/drone/Drone.H | 7 +- src/wind_energy/actuator/drone/Drone.cpp | 7 +- src/wind_energy/actuator/drone/drone_ops.cpp | 51 +++- .../actuator/motion/CMakeLists.txt | 7 + .../actuator/motion/RigidBodyMotion.H | 50 ++++ .../actuator/motion/RigidBodyMotion.cpp | 250 ++++++++++++++++++ src/wind_energy/actuator/motion/RotorMotion.H | 37 +++ .../actuator/motion/RotorMotion.cpp | 107 ++++++++ src/wind_energy/actuator/motion/TimeTable.H | 50 ++++ src/wind_energy/actuator/motion/TimeTable.cpp | 178 +++++++++++++ .../actuator/sector/ActuatorSector.H | 14 + .../actuator/sector/actuator_sector_ops.cpp | 112 +++++--- .../wind_energy/actuator/CMakeLists.txt | 1 + .../actuator/test_actuator_motion.cpp | 203 ++++++++++++++ 18 files changed, 1163 insertions(+), 77 deletions(-) create mode 100644 src/core/vs/quaternion.H create mode 100644 src/wind_energy/actuator/motion/CMakeLists.txt create mode 100644 src/wind_energy/actuator/motion/RigidBodyMotion.H create mode 100644 src/wind_energy/actuator/motion/RigidBodyMotion.cpp create mode 100644 src/wind_energy/actuator/motion/RotorMotion.H create mode 100644 src/wind_energy/actuator/motion/RotorMotion.cpp create mode 100644 src/wind_energy/actuator/motion/TimeTable.H create mode 100644 src/wind_energy/actuator/motion/TimeTable.cpp create mode 100644 unit_tests/wind_energy/actuator/test_actuator_motion.cpp diff --git a/src/core/vs/quaternion.H b/src/core/vs/quaternion.H new file mode 100644 index 0000000000..223d9b7804 --- /dev/null +++ b/src/core/vs/quaternion.H @@ -0,0 +1,136 @@ +#ifndef VS_QUATERNION_H +#define VS_QUATERNION_H + +#include "src/core/vs/tensor.H" +#include "src/utilities/constants.H" +#include "src/utilities/trig_ops.H" + +#include +#include +#include + +#include "AMReX.H" + +namespace kynema_sgf::vs { + +/** Scalar-first unit quaternion used to interpolate and compose rotations. */ +struct Quaternion +{ + amrex::Real w{1.0}; + amrex::Real x{0.0}; + amrex::Real y{0.0}; + amrex::Real z{0.0}; +}; + +inline Quaternion normalized(Quaternion q) +{ + const amrex::Real norm = + std::sqrt(q.w * q.w + q.x * q.x + q.y * q.y + q.z * q.z); + if (norm <= std::numeric_limits::epsilon()) { + amrex::Abort("Cannot normalize a zero quaternion"); + } + q.w /= norm; + q.x /= norm; + q.y /= norm; + q.z /= norm; + return q; +} + +inline amrex::Real dot(const Quaternion& a, const Quaternion& b) +{ + return a.w * b.w + a.x * b.x + a.y * b.y + a.z * b.z; +} + +inline Quaternion conjugate(const Quaternion& q) +{ + return {q.w, -q.x, -q.y, -q.z}; +} + +//! Rotation composition: apply rhs first, followed by lhs +inline Quaternion operator*(const Quaternion& lhs, const Quaternion& rhs) +{ + return { + lhs.w * rhs.w - lhs.x * rhs.x - lhs.y * rhs.y - lhs.z * rhs.z, + lhs.w * rhs.x + lhs.x * rhs.w + lhs.y * rhs.z - lhs.z * rhs.y, + lhs.w * rhs.y - lhs.x * rhs.z + lhs.y * rhs.w + lhs.z * rhs.x, + lhs.w * rhs.z + lhs.x * rhs.y - lhs.y * rhs.x + lhs.z * rhs.w}; +} + +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Tensor tensor_unit(const Quaternion& q) +{ + return {1.0 - 2.0 * (q.y * q.y + q.z * q.z), 2.0 * (q.x * q.y - q.z * q.w), + 2.0 * (q.x * q.z + q.y * q.w), 2.0 * (q.x * q.y + q.z * q.w), + 1.0 - 2.0 * (q.x * q.x + q.z * q.z), 2.0 * (q.y * q.z - q.x * q.w), + 2.0 * (q.x * q.z - q.y * q.w), 2.0 * (q.y * q.z + q.x * q.w), + 1.0 - 2.0 * (q.x * q.x + q.y * q.y)}; +} + +inline Tensor tensor(const Quaternion& input) +{ + return tensor_unit(normalized(input)); +} + +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Quaternion +from_axis_angle(const Vector& axis, const amrex::Real angle_degrees) +{ + const amrex::Real half_angle = -0.5 * utils::radians(angle_degrees); + const amrex::Real axis_magnitude = mag(axis); + const amrex::Real scale = std::sin(half_angle) / axis_magnitude; + return { + std::cos(half_angle), scale * axis.x(), scale * axis.y(), + scale * axis.z()}; +} + +//! Compatibility interface for existing axis-angle tensor rotations +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Tensor +quaternion(const Vector& axis, const amrex::Real angle_degrees) +{ + return tensor_unit(from_axis_angle(axis, angle_degrees)); +} + +/** Construct Rz(yaw) Ry(pitch) Rx(roll), with angles in degrees. */ +inline Quaternion from_roll_pitch_yaw(const Vector& angles) +{ + // The existing vs rotation tensors use the opposite sign from the + // conventional active quaternion formula. + const amrex::Real roll = -0.5 * utils::radians(angles.x()); + const amrex::Real pitch = -0.5 * utils::radians(angles.y()); + const amrex::Real yaw = -0.5 * utils::radians(angles.z()); + const amrex::Real cr = std::cos(roll); + const amrex::Real sr = std::sin(roll); + const amrex::Real cp = std::cos(pitch); + const amrex::Real sp = std::sin(pitch); + const amrex::Real cy = std::cos(yaw); + const amrex::Real sy = std::sin(yaw); + return normalized( + {cr * cp * cy + sr * sp * sy, sr * cp * cy - cr * sp * sy, + cr * sp * cy + sr * cp * sy, cr * cp * sy - sr * sp * cy}); +} + +inline Quaternion slerp(Quaternion a, Quaternion b, const amrex::Real fraction) +{ + a = normalized(a); + b = normalized(b); + amrex::Real cosine = dot(a, b); + if (cosine < 0.0) { + b = {-b.w, -b.x, -b.y, -b.z}; + cosine = -cosine; + } + // Scale the near-parallel cutoff with the configured Real precision. + if (cosine > 1.0 - std::sqrt(constants::EPS)) { + return normalized( + {a.w + fraction * (b.w - a.w), a.x + fraction * (b.x - a.x), + a.y + fraction * (b.y - a.y), a.z + fraction * (b.z - a.z)}); + } + const amrex::Real angle = std::acos(std::clamp(cosine, -1.0, 1.0)); + const amrex::Real denom = std::sin(angle); + const amrex::Real wa = std::sin((1.0 - fraction) * angle) / denom; + const amrex::Real wb = std::sin(fraction * angle) / denom; + return { + wa * a.w + wb * b.w, wa * a.x + wb * b.x, wa * a.y + wb * b.y, + wa * a.z + wb * b.z}; +} + +} // namespace kynema_sgf::vs + +#endif diff --git a/src/core/vs/tensor.H b/src/core/vs/tensor.H index 540140da75..b4e1c1945d 100644 --- a/src/core/vs/tensor.H +++ b/src/core/vs/tensor.H @@ -160,5 +160,6 @@ using Tensor = TensorT; } // namespace kynema_sgf::vs #include "src/core/vs/tensorI.H" +#include "src/core/vs/quaternion.H" #endif /* VS_TENSOR_H */ diff --git a/src/core/vs/tensorI.H b/src/core/vs/tensorI.H index 7abc72ac9b..d9d71d5852 100644 --- a/src/core/vs/tensorI.H +++ b/src/core/vs/tensorI.H @@ -214,34 +214,6 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Tensor zrot(const amrex::Real angle) 0.0_rt, 0.0_rt, 0.0_rt, 1.0_rt}; } -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Tensor -quaternion(const Vector& axis, const amrex::Real angle) -{ - const amrex::Real ang = -1.0_rt * utils::radians(angle); - const amrex::Real cval = std::cos(0.5_rt * ang); - const amrex::Real sval = std::sin(0.5_rt * ang); - const auto vmag = mag(axis); - const amrex::Real q0 = cval; - const amrex::Real q1 = sval * axis.x() / vmag; - const amrex::Real q2 = sval * axis.y() / vmag; - const amrex::Real q3 = sval * axis.z() / vmag; - - Tensor t; - t.xx() = (q0 * q0) + (q1 * q1) - (q2 * q2) - (q3 * q3); - t.xy() = 2.0_rt * (q1 * q2 - q0 * q3); - t.xz() = 2.0_rt * (q0 * q2 + q1 * q3); - - t.yx() = 2.0_rt * (q1 * q2 + q0 * q3); - t.yy() = (q0 * q0) - (q1 * q1) + (q2 * q2) - (q3 * q3); - t.yz() = 2.0_rt * (q2 * q3 - q0 * q1); - - t.zx() = 2.0_rt * (q1 * q3 - q0 * q2); - t.zy() = 2.0_rt * (q0 * q1 + q2 * q3); - t.zz() = (q0 * q0) - (q1 * q1) - (q2 * q2) + (q3 * q3); - - return t; -} - } // namespace kynema_sgf::vs #endif /* VS_TENSORI_H */ diff --git a/src/wind_energy/actuator/CMakeLists.txt b/src/wind_energy/actuator/CMakeLists.txt index 2b115aa85a..9aadfa4b48 100644 --- a/src/wind_energy/actuator/CMakeLists.txt +++ b/src/wind_energy/actuator/CMakeLists.txt @@ -8,6 +8,7 @@ target_sources(${kynema_sgf_lib_name} ) add_subdirectory(aero) +add_subdirectory(motion) add_subdirectory(wing) add_subdirectory(sector) add_subdirectory(drone) diff --git a/src/wind_energy/actuator/drone/Drone.H b/src/wind_energy/actuator/drone/Drone.H index 329f240ebd..de2a6e1654 100644 --- a/src/wind_energy/actuator/drone/Drone.H +++ b/src/wind_energy/actuator/drone/Drone.H @@ -4,12 +4,15 @@ #include "src/wind_energy/actuator/actuator_types.H" #include "src/wind_energy/actuator/sector/ActuatorSector.H" #include "src/wind_energy/actuator/sector/actuator_sector_ops.H" +#include "src/wind_energy/actuator/motion/RigidBodyMotion.H" +#include "src/wind_energy/actuator/motion/RotorMotion.H" #include #include namespace kynema_sgf::actuator { +//! One child actuator sector and its fixed offset in the drone body frame struct DroneRotor { using SectorData = ActuatorSector::DataType; @@ -22,6 +25,7 @@ struct DroneRotor vs::Vector body_offset{vs::Vector::zero()}; }; +//! Composite actuator that gives all child rotors one rigid-body motion struct DroneData { int num_rotors{0}; @@ -33,7 +37,8 @@ struct DroneData vs::Vector translation_velocity{vs::Vector::zero()}; vs::Vector body_orientation_degrees{vs::Vector::zero()}; vs::Tensor body_orientation{vs::Tensor::identity()}; - RealList rotor_omegas; + std::shared_ptr body_motion; + std::shared_ptr rotor_motion; amrex::Vector mirror_blades; vs::Vector total_force{vs::Vector::zero()}; vs::Vector total_moment{vs::Vector::zero()}; diff --git a/src/wind_energy/actuator/drone/Drone.cpp b/src/wind_energy/actuator/drone/Drone.cpp index 09f8d935c1..07003b9f46 100644 --- a/src/wind_energy/actuator/drone/Drone.cpp +++ b/src/wind_energy/actuator/drone/Drone.cpp @@ -1,9 +1,9 @@ #include "src/wind_energy/actuator/drone/Drone.H" #include "src/wind_energy/actuator/drone/drone_ops.H" #include "src/wind_energy/actuator/ActuatorModel.H" +#include "src/utilities/trig_ops.H" #include -#include namespace kynema_sgf::actuator { @@ -15,6 +15,7 @@ namespace drone { vs::Tensor body_rotation(const vs::Vector& angles) { + // Apply intrinsic body x, then y, then z rotations to body-frame vectors. return vs::zrot(angles.z()) & vs::yrot(angles.y()) & vs::xrot(angles.x()); } @@ -32,10 +33,12 @@ VecList rotor_body_offsets(const RealList& lengths, const RealList& angles_degrees) { AMREX_ALWAYS_ASSERT(lengths.size() == angles_degrees.size()); + // Arms lie in the body x-y plane; rigid-body motion maps these fixed + // offsets into the CFD frame at runtime. VecList offsets(lengths.size()); for (int i = 0; i < static_cast(lengths.size()); ++i) { const amrex::Real angle = - angles_degrees[i] * std::numbers::pi_v / 180.0_rt; + ::kynema_sgf::utils::radians(angles_degrees[i]); offsets[i] = { lengths[i] * std::cos(angle), lengths[i] * std::sin(angle), 0.0_rt}; } diff --git a/src/wind_energy/actuator/drone/drone_ops.cpp b/src/wind_energy/actuator/drone/drone_ops.cpp index 2a61f691a5..a1a8fdcc5d 100644 --- a/src/wind_energy/actuator/drone/drone_ops.cpp +++ b/src/wind_energy/actuator/drone/drone_ops.cpp @@ -44,11 +44,29 @@ void ReadInputsOp::operator()( if (meta.num_rotors <= 0) { input_error(label, "num_rotors must be positive"); } - pp.get("center", meta.center); + if (!pp.contains("position_timetable")) { + pp.get("center", meta.center); + } pp.query("translation_velocity", meta.translation_velocity); pp.query("body_orientation_degrees", meta.body_orientation_degrees); + if (pp.contains("orientation_timetable") && + pp.contains("body_orientation_degrees")) { + input_error( + label, + "body_orientation_degrees cannot be combined with " + "orientation_timetable"); + } meta.body_orientation = drone::body_rotation(meta.body_orientation_degrees); + // All child sectors share these histories, keeping the rigid-body pose and + // rotor clocks synchronized across the composite actuator. + meta.body_motion = std::make_shared(); + meta.body_motion->read_inputs( + pp, meta.center, meta.body_orientation, "angular_velocity"); + meta.rotor_motion = std::make_shared(); + meta.rotor_motion->read_drone_inputs(pp, meta.num_rotors); + // Instance inputs take precedence over Drone defaults; a scalar arm length + // expands to every rotor while the list permits asymmetric layouts. const auto& instance = pp.params(); const auto& defaults = pp.default_params(); const bool instance_length = instance.contains("arm_length"); @@ -81,6 +99,8 @@ void ReadInputsOp::operator()( input_error(label, "all arm lengths must be positive"); } + // Explicit angles describe irregular layouts. Otherwise phase rotates a + // uniformly spaced layout about the body z axis. const bool instance_angles = instance.contains("arm_angles_degrees"); const bool instance_phase = instance.contains("arm_phase_degrees"); if (instance_angles && instance_phase) { @@ -112,11 +132,6 @@ void ReadInputsOp::operator()( } validate_distinct_angles(label, meta.arm_angles_degrees); - pp.getarr("rotor_omegas", meta.rotor_omegas); - if (static_cast(meta.rotor_omegas.size()) != meta.num_rotors) { - input_error(label, "rotor_omegas must contain num_rotors values"); - } - meta.mirror_blades.assign(meta.num_rotors, false); if (pp.contains("mirror_blades")) { amrex::Vector mirror_inputs; @@ -136,7 +151,6 @@ void ReadInputsOp::operator()( const auto offsets = drone::rotor_body_offsets(meta.arm_lengths, meta.arm_angles_degrees); - const auto rotor_normal = meta.body_orientation & vs::Vector::khat(); meta.rotors.reserve(meta.num_rotors); for (int i = 0; i < meta.num_rotors; ++i) { const std::string rotor_label = label + ".R" + std::to_string(i + 1); @@ -148,22 +162,34 @@ void ReadInputsOp::operator()( // instance namespaces for both drone geometry and sector aerodynamics; // each reader simply ignores inputs outside its responsibility. utils::ActParser sector_pp("Actuator.Drone", "Actuator." + label); - rotor->data.meta().omega = meta.rotor_omegas[i]; - rotor->data.meta().user_omega = true; + auto& rotor_meta = rotor->data.meta(); + rotor_meta.body_motion = meta.body_motion; + rotor_meta.rotor_motion = meta.rotor_motion; + rotor_meta.rotor_index = i; + rotor_meta.body_offset = offsets[i]; + rotor_meta.omega = meta.rotor_motion->omega(i, 0.0_rt); + rotor_meta.user_omega = true; ReadInputsOp()(rotor->data, sector_pp); if (meta.mirror_blades[i]) { + // Mirroring the blade geometry reverses twist without introducing + // a separate rotation-direction convention. for (auto& twist : rotor->data.meta().twist_inp) { twist = -twist; } } const auto rotor_center = - meta.center + (meta.body_orientation & rotor->body_offset); + meta.body_motion->position(0.0_rt) + + (meta.body_motion->orientation(0.0_rt) & rotor->body_offset); + const auto initial_normal = + meta.body_motion->orientation(0.0_rt) & vs::Vector::khat(); sector::set_placement( - rotor->data, rotor_center, rotor_normal, meta.translation_velocity); + rotor->data, rotor_center, initial_normal, + meta.body_motion->translation_velocity(0.0_rt)); rotor->output.read_io_options(sector_pp); meta.rotors.emplace_back(std::move(rotor)); } + // The composite search region is the union of its child-sector regions. amrex::GpuArray lo; amrex::GpuArray hi; for (int n = 0; n < AMREX_SPACEDIM; ++n) { @@ -225,8 +251,7 @@ void ComputeForceOp::operator()(Drone::DataType& data) const auto& time = data.sim().time(); const amrex::Real midpoint_time = time.current_time() + 0.5_rt * sector::timestep_width(data.sim()); - const auto drone_center = - meta.center + meta.translation_velocity * midpoint_time; + const auto drone_center = meta.body_motion->position(midpoint_time); meta.total_force = vs::Vector::zero(); meta.total_moment = vs::Vector::zero(); for (auto& rotor : meta.rotors) { diff --git a/src/wind_energy/actuator/motion/CMakeLists.txt b/src/wind_energy/actuator/motion/CMakeLists.txt new file mode 100644 index 0000000000..5bb1fb22d8 --- /dev/null +++ b/src/wind_energy/actuator/motion/CMakeLists.txt @@ -0,0 +1,7 @@ +target_sources(${kynema_sgf_lib_name} + PRIVATE + + TimeTable.cpp + RigidBodyMotion.cpp + RotorMotion.cpp + ) diff --git a/src/wind_energy/actuator/motion/RigidBodyMotion.H b/src/wind_energy/actuator/motion/RigidBodyMotion.H new file mode 100644 index 0000000000..80bdcec1e3 --- /dev/null +++ b/src/wind_energy/actuator/motion/RigidBodyMotion.H @@ -0,0 +1,50 @@ +#ifndef ACTUATOR_RIGID_BODY_MOTION_H +#define ACTUATOR_RIGID_BODY_MOTION_H + +#include "src/core/vs/quaternion.H" +#include "src/wind_energy/actuator/ActParser.H" +#include "src/wind_energy/actuator/motion/TimeTable.H" + +namespace kynema_sgf::actuator::motion { + +/** Prescribed translation and orientation of a rigid body. + * + * Positions and velocities are in the global CFD frame. The orientation is + * the active rotation that maps body-frame vectors into that frame. + */ +class RigidBodyMotion +{ +public: + void read_inputs( + const utils::ActParser& pp, + const vs::Vector& initial_position, + const vs::Tensor& initial_orientation, + const std::string& angular_velocity_key = "angular_velocity", + const vs::Vector& default_angular_velocity = vs::Vector::zero()); + + [[nodiscard]] vs::Vector position(amrex::Real time) const; + [[nodiscard]] vs::Vector translation_velocity(amrex::Real time) const; + [[nodiscard]] vs::Tensor orientation(amrex::Real time) const; + [[nodiscard]] vs::Vector angular_velocity(amrex::Real time) const; + + [[nodiscard]] bool moves() const; + +private: + enum class OrientationFormat { RollPitchYaw, Quaternion }; + + [[nodiscard]] vs::Quaternion orientation_quaternion(int row) const; + + vs::Vector m_initial_position{vs::Vector::zero()}; + vs::Vector m_constant_velocity{vs::Vector::zero()}; + vs::Tensor m_initial_orientation{vs::Tensor::identity()}; + vs::Vector m_constant_angular_velocity{vs::Vector::zero()}; + TimeTable m_position_table; + TimeTable m_velocity_table; + TimeTable m_orientation_table; + TimeTable m_angular_velocity_table; + OrientationFormat m_orientation_format{OrientationFormat::RollPitchYaw}; +}; + +} // namespace kynema_sgf::actuator::motion + +#endif diff --git a/src/wind_energy/actuator/motion/RigidBodyMotion.cpp b/src/wind_energy/actuator/motion/RigidBodyMotion.cpp new file mode 100644 index 0000000000..63a2f36bed --- /dev/null +++ b/src/wind_energy/actuator/motion/RigidBodyMotion.cpp @@ -0,0 +1,250 @@ +#include "src/wind_energy/actuator/motion/RigidBodyMotion.H" + +#include "src/core/vs/quaternion.H" +#include "src/wind_energy/actuator/sector/actuator_sector_ops.H" + +#include +#include +#include + +#include "AMReX.H" + +using namespace amrex::literals; + +namespace kynema_sgf::actuator::motion { +namespace { + +vs::Vector vector3(const RealList& values) +{ + return {values[0], values[1], values[2]}; +} + +vs::Tensor integrate_global_angular_velocity( + const TimeTable& table, + const vs::Tensor& initial_orientation, + const amrex::Real time) +{ + const auto& times = table.times(); + vs::Tensor result = initial_orientation; + amrex::Real start = times.front(); + const amrex::Real direction = (time >= start) ? 1.0_rt : -1.0_rt; + // Compose global-frame rotations one table interval at a time. Midpoint + // sampling gives a second-order update for a varying angular velocity. + while (direction * (time - start) > 0.0_rt) { + amrex::Real end = time; + if (direction > 0.0_rt) { + const auto next = + std::upper_bound(times.begin(), times.end(), start); + if (next != times.end()) { + end = std::min(time, *next); + } + } else { + end = time; + } + const amrex::Real midpoint = 0.5_rt * (start + end); + const auto omega = vector3(table.value(midpoint)); + result = + sector::rotation_matrix_from_vector(omega * (end - start)) & result; + start = end; + } + return result; +} + +} // namespace + +void RigidBodyMotion::read_inputs( + const utils::ActParser& pp, + const vs::Vector& initial_position, + const vs::Tensor& initial_orientation, + const std::string& angular_velocity_key, + const vs::Vector& default_angular_velocity) +{ + m_initial_position = initial_position; + m_initial_orientation = initial_orientation; + m_constant_angular_velocity = default_angular_velocity; + std::string extrapolation{"hold"}; + pp.query("timetable_extrapolation", extrapolation); + + // Translation may be prescribed by position or velocity, but never both. + const bool has_position = pp.contains("position_timetable"); + const bool has_velocity_table = pp.contains("velocity_timetable"); + const bool has_velocity = pp.contains("translation_velocity"); + if (static_cast(has_position) + static_cast(has_velocity_table) + + static_cast(has_velocity) > + 1) { + amrex::Abort( + "Specify only one of position_timetable, velocity_timetable, and " + "translation_velocity"); + } + if (has_position && pp.contains("center")) { + amrex::Abort("center cannot be combined with position_timetable"); + } + if (has_velocity_table && !pp.contains("center")) { + amrex::Abort("center is required with velocity_timetable"); + } + if (has_position) { + std::string filename; + pp.get("position_timetable", filename); + m_position_table.read(filename, 3, extrapolation); + } else if (has_velocity_table) { + std::string filename; + pp.get("velocity_timetable", filename); + m_velocity_table.read(filename, 3, extrapolation); + } else { + pp.query("translation_velocity", m_constant_velocity); + } + + // Rotation follows the same exclusive-source rule as translation. + const bool has_orientation = pp.contains("orientation_timetable"); + const bool has_angular_table = pp.contains("angular_velocity_timetable"); + const bool has_angular = pp.contains(angular_velocity_key); + if (static_cast(has_orientation) + + static_cast(has_angular_table) + + static_cast(has_angular) > + 1) { + amrex::Abort( + "Specify only one of orientation_timetable, " + "angular_velocity_timetable, and constant angular velocity"); + } + std::string angular_frame{"global"}; + pp.query("angular_velocity_frame", angular_frame); + if (amrex::toLower(angular_frame) != "global") { + amrex::Abort( + "The first motion implementation supports only global " + "angular_velocity_frame"); + } + if (has_orientation) { + std::string filename; + pp.get("orientation_timetable", filename); + std::string format{"roll_pitch_yaw"}; + pp.query("orientation_format", format); + format = amrex::toLower(format); + if (format == "roll_pitch_yaw") { + m_orientation_format = OrientationFormat::RollPitchYaw; + m_orientation_table.read(filename, 3, extrapolation); + } else if (format == "quaternion") { + m_orientation_format = OrientationFormat::Quaternion; + m_orientation_table.read(filename, 4, extrapolation); + } else { + amrex::Abort( + "orientation_format must be roll_pitch_yaw or quaternion"); + } + } else if (has_angular_table) { + std::string filename; + pp.get("angular_velocity_timetable", filename); + m_angular_velocity_table.read(filename, 3, extrapolation); + } else { + pp.query(angular_velocity_key, m_constant_angular_velocity); + } +} + +vs::Vector RigidBodyMotion::position(const amrex::Real time) const +{ + if (!m_position_table.empty()) { + return vector3(m_position_table.value(time)); + } + if (!m_velocity_table.empty()) { + return m_initial_position + vector3(m_velocity_table.integral(time)); + } + return m_initial_position + m_constant_velocity * time; +} + +vs::Vector RigidBodyMotion::translation_velocity(const amrex::Real time) const +{ + if (!m_position_table.empty()) { + return vector3(m_position_table.derivative(time)); + } + if (!m_velocity_table.empty()) { + return vector3(m_velocity_table.value(time)); + } + return m_constant_velocity; +} + +vs::Quaternion RigidBodyMotion::orientation_quaternion(const int row) const +{ + const auto values = m_orientation_table.row(row); + if (m_orientation_format == OrientationFormat::Quaternion) { + return vs::normalized({values[0], values[1], values[2], values[3]}); + } + return vs::from_roll_pitch_yaw({values[0], values[1], values[2]}); +} + +vs::Tensor RigidBodyMotion::orientation(const amrex::Real time) const +{ + if (!m_orientation_table.empty()) { + const auto& times = m_orientation_table.times(); + if (time <= times.front()) { + static_cast(m_orientation_table.value(time)); + return vs::tensor(orientation_quaternion(0)); + } + if (time >= times.back()) { + static_cast(m_orientation_table.value(time)); + return vs::tensor( + orientation_quaternion(static_cast(times.size()) - 1)); + } + const int upper = static_cast( + std::upper_bound(times.begin(), times.end(), time) - times.begin()); + const int lower = upper - 1; + const amrex::Real fraction = + (time - times[lower]) / (times[upper] - times[lower]); + return vs::tensor( + vs::slerp( + orientation_quaternion(lower), orientation_quaternion(upper), + fraction)); + } + if (!m_angular_velocity_table.empty()) { + return integrate_global_angular_velocity( + m_angular_velocity_table, m_initial_orientation, time); + } + return sector::rotation_matrix_from_vector( + m_constant_angular_velocity * time) & + m_initial_orientation; +} + +vs::Vector RigidBodyMotion::angular_velocity(const amrex::Real time) const +{ + if (!m_orientation_table.empty()) { + const auto& times = m_orientation_table.times(); + if (times.size() == 1 || time < times.front() || time > times.back()) { + static_cast(m_orientation_table.value(time)); + return vs::Vector::zero(); + } + int upper = static_cast( + std::upper_bound(times.begin(), times.end(), time) - times.begin()); + upper = std::min(upper, static_cast(times.size()) - 1); + const int lower = upper - 1; + auto a = orientation_quaternion(lower); + auto b = orientation_quaternion(upper); + if (vs::dot(a, b) < 0.0_rt) { + b = {-b.w, -b.x, -b.y, -b.z}; + } + const auto relative = vs::normalized(b * vs::conjugate(a)); + // Convert the relative quaternion over this interval into a constant + // body-frame angular velocity, then express it in the CFD frame. + const amrex::Real half_angle = + std::acos(std::clamp(relative.w, -1.0_rt, 1.0_rt)); + const amrex::Real sine = std::sin(half_angle); + if (std::abs(sine) <= std::numeric_limits::epsilon()) { + return vs::Vector::zero(); + } + const amrex::Real scale = + 2.0_rt * half_angle / (sine * (times[upper] - times[lower])); + const vs::Vector body_omega{ + relative.x * scale, relative.y * scale, relative.z * scale}; + return orientation(time) & body_omega; + } + if (!m_angular_velocity_table.empty()) { + return vector3(m_angular_velocity_table.value(time)); + } + return m_constant_angular_velocity; +} + +bool RigidBodyMotion::moves() const +{ + return !m_position_table.empty() || !m_velocity_table.empty() || + !m_orientation_table.empty() || !m_angular_velocity_table.empty() || + vs::mag(m_constant_velocity) > 0.0_rt || + vs::mag(m_constant_angular_velocity) > 0.0_rt; +} + +} // namespace kynema_sgf::actuator::motion diff --git a/src/wind_energy/actuator/motion/RotorMotion.H b/src/wind_energy/actuator/motion/RotorMotion.H new file mode 100644 index 0000000000..e73bdbdcd2 --- /dev/null +++ b/src/wind_energy/actuator/motion/RotorMotion.H @@ -0,0 +1,37 @@ +#ifndef ACTUATOR_ROTOR_MOTION_H +#define ACTUATOR_ROTOR_MOTION_H + +#include "src/wind_energy/actuator/ActParser.H" +#include "src/wind_energy/actuator/motion/TimeTable.H" + +namespace kynema_sgf::actuator::motion { + +/** Signed rotor-speed histories and their corresponding azimuths. */ +class RotorMotion +{ +public: + void read_drone_inputs(const utils::ActParser& pp, int num_rotors); + void read_sector_inputs( + const utils::ActParser& pp, + amrex::Real preset_omega, + bool omega_is_preset); + + [[nodiscard]] int num_rotors() const + { + return static_cast(m_constant_omegas.size()); + } + [[nodiscard]] amrex::Real omega(int rotor, amrex::Real time) const; + [[nodiscard]] amrex::Real azimuth(int rotor, amrex::Real time) const; + +private: + void read_initial_azimuths(const utils::ActParser& pp, int num_rotors); + + RealList m_constant_omegas; + //! Initial rotor angles stored internally in radians + RealList m_initial_azimuths; + TimeTable m_speed_table; +}; + +} // namespace kynema_sgf::actuator::motion + +#endif diff --git a/src/wind_energy/actuator/motion/RotorMotion.cpp b/src/wind_energy/actuator/motion/RotorMotion.cpp new file mode 100644 index 0000000000..5da7b1fa8d --- /dev/null +++ b/src/wind_energy/actuator/motion/RotorMotion.cpp @@ -0,0 +1,107 @@ +#include "src/wind_energy/actuator/motion/RotorMotion.H" + +#include "src/utilities/trig_ops.H" + +#include "AMReX.H" + +using namespace amrex::literals; + +namespace kynema_sgf::actuator::motion { + +void RotorMotion::read_initial_azimuths( + const utils::ActParser& pp, const int num_rotors) +{ + const bool shared = pp.contains("initial_azimuth_degrees"); + const bool individual = pp.contains("initial_azimuths_degrees"); + if (shared && individual) { + amrex::Abort( + "initial_azimuth_degrees and initial_azimuths_degrees are mutually " + "exclusive"); + } + m_initial_azimuths.assign(num_rotors, 0.0_rt); + if (shared) { + amrex::Real value = 0.0_rt; + pp.get("initial_azimuth_degrees", value); + m_initial_azimuths.assign( + num_rotors, ::kynema_sgf::utils::radians(value)); + } else if (individual) { + pp.getarr("initial_azimuths_degrees", m_initial_azimuths); + if (static_cast(m_initial_azimuths.size()) != num_rotors) { + amrex::Abort( + "initial_azimuths_degrees must contain num_rotors values"); + } + for (auto& value : m_initial_azimuths) { + value = ::kynema_sgf::utils::radians(value); + } + } +} + +void RotorMotion::read_drone_inputs( + const utils::ActParser& pp, const int num_rotors) +{ + const bool constant = pp.contains("rotor_omegas"); + const bool timetable = pp.contains("rotor_speed_timetable"); + if (constant == timetable) { + amrex::Abort( + "Specify exactly one of rotor_omegas and rotor_speed_timetable"); + } + m_constant_omegas.assign(num_rotors, 0.0_rt); + if (constant) { + pp.getarr("rotor_omegas", m_constant_omegas); + if (static_cast(m_constant_omegas.size()) != num_rotors) { + amrex::Abort("rotor_omegas must contain num_rotors values"); + } + } else { + std::string filename; + std::string extrapolation{"hold"}; + pp.get("rotor_speed_timetable", filename); + pp.query("timetable_extrapolation", extrapolation); + m_speed_table.read(filename, num_rotors, extrapolation); + } + read_initial_azimuths(pp, num_rotors); +} + +void RotorMotion::read_sector_inputs( + const utils::ActParser& pp, + const amrex::Real preset_omega, + const bool omega_is_preset) +{ + const bool constant = pp.contains("omega") || omega_is_preset; + const bool timetable = pp.contains("rotor_speed_timetable"); + if (constant == timetable) { + amrex::Abort( + "Specify exactly one of omega and rotor_speed_timetable for an " + "ActuatorSector"); + } + m_constant_omegas.assign(1, preset_omega); + if (!omega_is_preset && pp.contains("omega")) { + pp.get("omega", m_constant_omegas[0]); + } else if (timetable) { + std::string filename; + std::string extrapolation{"hold"}; + pp.get("rotor_speed_timetable", filename); + pp.query("timetable_extrapolation", extrapolation); + m_speed_table.read(filename, 1, extrapolation); + } + read_initial_azimuths(pp, 1); +} + +amrex::Real RotorMotion::omega(const int rotor, const amrex::Real time) const +{ + AMREX_ALWAYS_ASSERT(rotor >= 0 && rotor < num_rotors()); + return m_speed_table.empty() ? m_constant_omegas[rotor] + : m_speed_table.value(time)[rotor]; +} + +amrex::Real RotorMotion::azimuth(const int rotor, const amrex::Real time) const +{ + AMREX_ALWAYS_ASSERT(rotor >= 0 && rotor < num_rotors()); + // Integrating omega preserves phase when speed varies; omega(t) * time + // would only be correct for a constant rotor speed. + const amrex::Real angle = m_speed_table.empty() + ? m_constant_omegas[rotor] * time + : m_speed_table.integral(time)[rotor]; + return m_initial_azimuths[rotor] + angle; +} + +} // namespace kynema_sgf::actuator::motion diff --git a/src/wind_energy/actuator/motion/TimeTable.H b/src/wind_energy/actuator/motion/TimeTable.H new file mode 100644 index 0000000000..8da3c759c0 --- /dev/null +++ b/src/wind_energy/actuator/motion/TimeTable.H @@ -0,0 +1,50 @@ +#ifndef ACTUATOR_MOTION_TIME_TABLE_H +#define ACTUATOR_MOTION_TIME_TABLE_H + +#include "src/wind_energy/actuator/actuator_types.H" + +#include + +namespace kynema_sgf::actuator::motion { + +/** Piecewise-linear, multi-component time history. + * + * Values, derivatives, and integrals use the same interpolation intervals so + * that motion reconstructed from a velocity history remains consistent. + */ +class TimeTable +{ +public: + void read( + const std::string& filename, + int num_values, + const std::string& extrapolation = "hold"); + + [[nodiscard]] bool empty() const { return m_time.empty(); } + [[nodiscard]] int num_values() const { return m_num_values; } + [[nodiscard]] const RealList& times() const { return m_time; } + [[nodiscard]] RealList row(int index) const; + + [[nodiscard]] RealList value(amrex::Real time) const; + [[nodiscard]] RealList derivative(amrex::Real time) const; + [[nodiscard]] RealList integral(amrex::Real time) const; + +private: + [[nodiscard]] int interval(amrex::Real time) const; + void warn_hold(amrex::Real time) const; + + std::string m_filename; + int m_num_values{0}; + RealList m_time; + //! Row-major values used directly by interp::linear's ncomp interface + RealList m_values; + //! Exact trapezoidal integral of each completed interpolation interval + amrex::Vector m_prefix_integral; + mutable bool m_warned_below{false}; + mutable bool m_warned_above{false}; + bool m_error_on_extrapolation{false}; +}; + +} // namespace kynema_sgf::actuator::motion + +#endif diff --git a/src/wind_energy/actuator/motion/TimeTable.cpp b/src/wind_energy/actuator/motion/TimeTable.cpp new file mode 100644 index 0000000000..f25b1d8b07 --- /dev/null +++ b/src/wind_energy/actuator/motion/TimeTable.cpp @@ -0,0 +1,178 @@ +#include "src/wind_energy/actuator/motion/TimeTable.H" + +#include "src/utilities/linear_interpolation.H" + +#include +#include + +#include "AMReX.H" +#include "AMReX_Print.H" + +using namespace amrex::literals; + +namespace kynema_sgf::actuator::motion { + +void TimeTable::read( + const std::string& filename, + const int num_values, + const std::string& extrapolation) +{ + m_filename = filename; + m_num_values = num_values; + m_time.clear(); + m_values.clear(); + m_prefix_integral.clear(); + m_warned_below = false; + m_warned_above = false; + const auto extrapolation_lower = amrex::toLower(extrapolation); + if ((extrapolation_lower != "hold") && (extrapolation_lower != "error")) { + amrex::Abort("Actuator timetable extrapolation must be hold or error"); + } + m_error_on_extrapolation = (extrapolation_lower == "error"); + std::ifstream stream(filename); + if (!stream.good()) { + amrex::Abort("Cannot open actuator motion timetable: " + filename); + } + + stream.ignore(std::numeric_limits::max(), '\n'); + amrex::Real time; + while (stream >> time) { + RealList row(num_values); + for (auto& item : row) { + if (!(stream >> item)) { + amrex::Abort( + "Invalid actuator motion timetable row in: " + filename); + } + } + if (!m_time.empty() && time <= m_time.back()) { + amrex::Abort( + "Actuator motion timetable times must be strictly " + "increasing: " + + filename); + } + m_time.push_back(time); + m_values.insert(m_values.end(), row.begin(), row.end()); + } + if (m_time.empty()) { + amrex::Abort("Actuator motion timetable contains no data: " + filename); + } + + m_prefix_integral.assign(m_time.size(), RealList(num_values, 0.0_rt)); + for (int i = 1; i < static_cast(m_time.size()); ++i) { + const amrex::Real dt = m_time[i] - m_time[i - 1]; + for (int n = 0; n < num_values; ++n) { + m_prefix_integral[i][n] = + m_prefix_integral[i - 1][n] + + 0.5_rt * dt * + (m_values[(i - 1) * m_num_values + n] + + m_values[i * m_num_values + n]); + } + } +} + +RealList TimeTable::row(const int index) const +{ + AMREX_ALWAYS_ASSERT(index >= 0 && index < static_cast(m_time.size())); + const auto begin = m_values.begin() + index * m_num_values; + return RealList(begin, begin + m_num_values); +} + +void TimeTable::warn_hold(const amrex::Real time) const +{ + if (m_error_on_extrapolation && + ((time < m_time.front()) || (time > m_time.back()))) { + amrex::Abort( + "Actuator motion timetable evaluated outside its time range: " + + m_filename); + } + if (time < m_time.front() && !m_warned_below) { + amrex::Print() << "WARNING: Holding first value of actuator motion " + "timetable '" + << m_filename << "' before time " << m_time.front() + << ".\n"; + m_warned_below = true; + } + if (time > m_time.back() && !m_warned_above) { + amrex::Print() << "WARNING: Holding last value of actuator motion " + "timetable '" + << m_filename << "' after time " << m_time.back() + << ".\n"; + m_warned_above = true; + } +} + +int TimeTable::interval(const amrex::Real time) const +{ + if (m_time.size() == 1 || time <= m_time.front()) { + return 0; + } + if (time >= m_time.back()) { + return static_cast(m_time.size()) - 2; + } + return ::kynema_sgf::interp::bisection_search( + m_time.begin(), m_time.end(), time) + .idx; +} + +RealList TimeTable::value(const amrex::Real time) const +{ + warn_hold(time); + RealList result(m_num_values); + for (int n = 0; n < m_num_values; ++n) { + result[n] = ::kynema_sgf::interp::linear( + m_time, m_values, time, m_num_values, n); + } + return result; +} + +RealList TimeTable::derivative(const amrex::Real time) const +{ + warn_hold(time); + RealList result(m_num_values, 0.0_rt); + if (m_time.size() == 1 || time < m_time.front() || time > m_time.back()) { + return result; + } + const int i = interval(time); + const amrex::Real dt = m_time[i + 1] - m_time[i]; + for (int n = 0; n < m_num_values; ++n) { + result[n] = (m_values[(i + 1) * m_num_values + n] - + m_values[i * m_num_values + n]) / + dt; + } + return result; +} + +RealList TimeTable::integral(const amrex::Real time) const +{ + warn_hold(time); + RealList result(m_num_values, 0.0_rt); + if (time <= m_time.front()) { + const amrex::Real dt = time - m_time.front(); + for (int n = 0; n < m_num_values; ++n) { + result[n] = dt * m_values[n]; + } + return result; + } + if (time >= m_time.back()) { + result = m_prefix_integral.back(); + const amrex::Real dt = time - m_time.back(); + for (int n = 0; n < m_num_values; ++n) { + result[n] += dt * m_values[(m_time.size() - 1) * m_num_values + n]; + } + return result; + } + const int i = interval(time); + result = m_prefix_integral[i]; + const amrex::Real dt = time - m_time[i]; + const amrex::Real interval_dt = m_time[i + 1] - m_time[i]; + for (int n = 0; n < m_num_values; ++n) { + const amrex::Real slope = (m_values[(i + 1) * m_num_values + n] - + m_values[i * m_num_values + n]) / + interval_dt; + result[n] += + m_values[i * m_num_values + n] * dt + 0.5_rt * slope * dt * dt; + } + return result; +} + +} // namespace kynema_sgf::actuator::motion diff --git a/src/wind_energy/actuator/sector/ActuatorSector.H b/src/wind_energy/actuator/sector/ActuatorSector.H index fa07028c7a..893832fbd1 100644 --- a/src/wind_energy/actuator/sector/ActuatorSector.H +++ b/src/wind_energy/actuator/sector/ActuatorSector.H @@ -3,6 +3,8 @@ #include "src/wind_energy/actuator/actuator_types.H" #include "src/wind_energy/actuator/aero/AirfoilTable.H" +#include "src/wind_energy/actuator/motion/RigidBodyMotion.H" +#include "src/wind_energy/actuator/motion/RotorMotion.H" #include #include @@ -26,6 +28,18 @@ namespace kynema_sgf::actuator { */ struct ActuatorSectorData { + //! Prescribed rigid-body motion shared with an owning Drone when present + std::shared_ptr body_motion; + + //! Prescribed rotor-speed motion shared with an owning Drone when present + std::shared_ptr rotor_motion; + + //! Rotor index within rotor_motion + int rotor_index{0}; + + //! Hub offset from the rigid-body center in body coordinates [m] + vs::Vector body_offset{vs::Vector::zero()}; + //! Number of rotor blades [-] int num_blades{2}; diff --git a/src/wind_energy/actuator/sector/actuator_sector_ops.cpp b/src/wind_energy/actuator/sector/actuator_sector_ops.cpp index 1be1e8e537..35ed8d6ec8 100644 --- a/src/wind_energy/actuator/sector/actuator_sector_ops.cpp +++ b/src/wind_energy/actuator/sector/actuator_sector_ops.cpp @@ -222,14 +222,16 @@ void update_midpoint_sample_points(ActuatorSector::DataType& data) // convention used elsewhere in the code. The resulting sampled velocity is // then used by ComputeForceOp for this step. const amrex::Real tmid = 0.5_rt * (time.current_time() + time.new_time()); - const amrex::Real dt = timestep_width(data.sim()); - meta.center = meta.center0 + meta.translation_velocity * tmid; - meta.azimuth = meta.omega * tmid; - meta.delta_azimuth = meta.omega * dt; - - const auto orientation = orientation_at_time( - tmid, orientation_matrix_from_normal(meta.rotor_normal), - meta.rotor_angular_velocity); + const auto orientation = meta.body_motion->orientation(tmid); + meta.center = + meta.body_motion->position(tmid) + (orientation & meta.body_offset); + meta.translation_velocity = meta.body_motion->translation_velocity(tmid); + meta.rotor_angular_velocity = meta.body_motion->angular_velocity(tmid); + meta.omega = meta.rotor_motion->omega(meta.rotor_index, tmid); + meta.azimuth = meta.rotor_motion->azimuth(meta.rotor_index, tmid); + meta.delta_azimuth = + meta.rotor_motion->azimuth(meta.rotor_index, time.new_time()) - + meta.rotor_motion->azimuth(meta.rotor_index, time.current_time()); const int nr = static_cast(meta.radius.size()); const int nvel = meta.num_blades * nr; @@ -272,8 +274,9 @@ void set_placement( const amrex::Real max_eps = local_epsilon(meta, max_interp_chord(meta)); const amrex::Real search_radius = meta.rotor_radius + meta.support_radius_over_epsilon * max_eps; - if (vs::mag(translation_velocity) > - std::numeric_limits::epsilon()) { + if ((meta.body_motion && meta.body_motion->moves()) || + vs::mag(translation_velocity) > + std::numeric_limits::epsilon()) { const auto& geom = data.sim().mesh().Geom(0); const auto plo = geom.ProbLoArray(); const auto phi = geom.ProbHiArray(); @@ -463,7 +466,7 @@ void ReadInputsOp::operator()( "ActuatorSector omega was set by its owning model and must " "not also be specified in the shared sector inputs"); } - } else { + } else if (!pp.contains("rotor_speed_timetable")) { pp.get("omega", meta.omega); } pp.get("airfoil_table", meta.airfoil_file); @@ -473,6 +476,10 @@ void ReadInputsOp::operator()( pp.query("center", meta.center0); pp.query("translation_velocity", meta.translation_velocity); pp.query("rotor_normal", meta.rotor_normal); + if (pp.contains("orientation_timetable") && pp.contains("rotor_normal")) { + amrex::Abort( + "rotor_normal cannot be combined with orientation_timetable"); + } if (pp.contains("rotor_orientation") && !pp.contains("rotor_normal")) { amrex::Print() << "WARNING: ActuatorSector input 'rotor_orientation' is " @@ -498,6 +505,15 @@ void ReadInputsOp::operator()( pp.get("rotor_angular_velocity", meta.rotor_angular_velocity); meta.user_rotor_angular_velocity = true; } + if (pp.contains("angular_velocity")) { + if (pp.contains("rotor_angular_velocity")) { + amrex::Abort( + "angular_velocity and rotor_angular_velocity are mutually " + "exclusive"); + } + pp.get("angular_velocity", meta.rotor_angular_velocity); + meta.user_rotor_angular_velocity = true; + } pp.queryarr("span_locs", meta.span_locs); pp.queryarr("chord", meta.chord_inp); pp.queryarr("twist", meta.twist_inp); @@ -555,6 +571,23 @@ void ReadInputsOp::operator()( meta.rotor_rotation_degrees_per_revolution; } + // Standalone sectors own their histories. Drone sectors arrive with shared + // histories already injected so every rotor uses the same body motion. + if (!meta.body_motion) { + meta.body_motion = std::make_shared(); + const std::string angular_velocity_key = pp.contains("angular_velocity") + ? "angular_velocity" + : "rotor_angular_velocity"; + meta.body_motion->read_inputs( + pp, meta.center0, + sector::orientation_matrix_from_normal(meta.rotor_normal), + angular_velocity_key, meta.rotor_angular_velocity); + } + if (!meta.rotor_motion) { + meta.rotor_motion = std::make_shared(); + meta.rotor_motion->read_sector_inputs(pp, meta.omega, meta.user_omega); + } + const amrex::Real max_eps = sector::local_epsilon(meta, sector::max_interp_chord(meta)); const amrex::Real search_radius = @@ -564,11 +597,14 @@ void ReadInputsOp::operator()( const auto phi = geom.ProbHiArray(); const auto& c = meta.center0; - const bool moves = vs::mag(meta.translation_velocity) > - std::numeric_limits::epsilon(); + const bool moves = (meta.body_motion && meta.body_motion->moves()) || + vs::mag(meta.translation_velocity) > + std::numeric_limits::epsilon(); const bool rotates = vs::mag(meta.rotor_angular_velocity) > std::numeric_limits::epsilon(); if (moves || rotates) { + // A prescribed trajectory is not generally bounded by its initial + // placement, so retain the actuator throughout the computational box. data.info().bound_box = amrex::RealBox( moves ? plo[0] : c.x() - search_radius, moves ? plo[1] : c.y() - search_radius, @@ -630,12 +666,16 @@ void ComputeForceOp::operator()( const amrex::Real t0 = time.current_time(); const amrex::Real dt = sector::timestep_width(data.sim()); const amrex::Real tmid = t0 + 0.5_rt * dt; - const amrex::Real start_theta = meta.omega * t0; - const amrex::Real dtheta = meta.omega * dt; - const auto initial_orientation = - sector::orientation_matrix_from_normal(meta.rotor_normal); - const auto mid_orientation = sector::orientation_at_time( - tmid, initial_orientation, meta.rotor_angular_velocity); + // Aerodynamic loads use the same midpoint pose as the CFD velocity samples. + const amrex::Real mid_theta = + meta.rotor_motion->azimuth(meta.rotor_index, tmid); + const auto mid_orientation = meta.body_motion->orientation(tmid); + const auto body_position = meta.body_motion->position(tmid); + const auto body_velocity = meta.body_motion->translation_velocity(tmid); + const auto body_angular_velocity = meta.body_motion->angular_velocity(tmid); + const auto hub_offset = mid_orientation & meta.body_offset; + meta.center = body_position + hub_offset; + meta.omega = meta.rotor_motion->omega(meta.rotor_index, tmid); const int nr = static_cast(meta.radius.size()); const int nvel = meta.num_blades * nr; @@ -659,14 +699,15 @@ void ComputeForceOp::operator()( vs::Vector e_theta; vs::Vector e_normal; sector::blade_basis( - start_theta + 0.5_rt * dtheta + phase, mid_orientation, e_r, - e_theta, e_normal); + mid_theta + phase, mid_orientation, e_r, e_theta, e_normal); const auto rel_pos = e_r * meta.radius[ir]; - const auto frame_vel = meta.rotor_angular_velocity ^ rel_pos; + // Rigid-body rotation acts about the drone center, so the lever arm + // includes both the rotor hub offset and blade-section radius. + const auto frame_vel = + body_angular_velocity ^ (hub_offset + rel_pos); const auto spin_vel = e_theta * (meta.omega * meta.radius[ir]); - const auto blade_vel = - meta.translation_velocity + frame_vel + spin_vel; + const auto blade_vel = body_velocity + frame_vel + spin_vel; const auto rel_wind = grid.vel[ip] - blade_vel; const amrex::Real vtheta = rel_wind & e_theta; const amrex::Real vnormal = rel_wind & e_normal; @@ -718,11 +759,13 @@ void ComputeForceOp::operator()( // theta_counts is one and the sector naturally reduces to actuator-line // behavior. int nforce = 0; + const amrex::Real body_swept_speed = + vs::mag(body_velocity) + + vs::mag(body_angular_velocity) * + (vs::mag(meta.body_offset) + meta.rotor_radius); for (int ir = 0; ir < nr; ++ir) { const amrex::Real swept_speed = - vs::mag(meta.translation_velocity) + - std::abs(meta.omega) * meta.radius[ir] + - vs::mag(meta.rotor_angular_velocity) * meta.radius[ir]; + body_swept_speed + std::abs(meta.omega) * meta.radius[ir]; const int ntheta = amrex::max( 1, static_cast(std::ceil( swept_speed * dt * meta.epsilon_dl / @@ -747,15 +790,18 @@ void ComputeForceOp::operator()( const amrex::Real xi = (static_cast(it) + 0.5_rt) / static_cast(ntheta); const amrex::Real t = t0 + xi * dt; - const amrex::Real theta = start_theta + xi * dtheta + phase; - const auto orientation = sector::orientation_at_time( - t, initial_orientation, meta.rotor_angular_velocity); + // Re-evaluate the prescribed pose and azimuth at every swept + // quadrature time rather than approximating the trajectory. + const amrex::Real theta = + meta.rotor_motion->azimuth(meta.rotor_index, t) + phase; + const auto orientation = meta.body_motion->orientation(t); + const auto rotor_center = meta.body_motion->position(t) + + (orientation & meta.body_offset); vs::Vector e_r; vs::Vector e_theta; vs::Vector e_normal; sector::blade_basis(theta, orientation, e_r, e_theta, e_normal); - grid.pos[iq] = meta.center0 + meta.translation_velocity * t + - e_r * meta.radius[ir]; + grid.pos[iq] = rotor_center + e_r * meta.radius[ir]; // Split the radial section force uniformly across the swept // quadrature points so the integrated force remains unchanged. const amrex::Real wt = @@ -765,7 +811,7 @@ void ComputeForceOp::operator()( meta.integrated_force = meta.integrated_force + grid.force[iq]; meta.integrated_moment = meta.integrated_moment + - ((grid.pos[iq] - meta.center) ^ grid.force[iq]); + ((grid.pos[iq] - rotor_center) ^ grid.force[iq]); grid.epsilon[iq] = vs::Vector::one() * meta.epsilon_profile[ir]; ++iq; } diff --git a/unit_tests/wind_energy/actuator/CMakeLists.txt b/unit_tests/wind_energy/actuator/CMakeLists.txt index 3c954ca640..0d74dc2a73 100644 --- a/unit_tests/wind_energy/actuator/CMakeLists.txt +++ b/unit_tests/wind_energy/actuator/CMakeLists.txt @@ -5,6 +5,7 @@ target_sources(${kynema_sgf_unit_test_exe_name} PRIVATE test_airfoil.cpp test_actuator_free_functions.cpp test_actuator_sector.cpp + test_actuator_motion.cpp test_drone.cpp test_disk_uniform_ct.cpp test_FLLC.cpp diff --git a/unit_tests/wind_energy/actuator/test_actuator_motion.cpp b/unit_tests/wind_energy/actuator/test_actuator_motion.cpp new file mode 100644 index 0000000000..b62a76b6cc --- /dev/null +++ b/unit_tests/wind_energy/actuator/test_actuator_motion.cpp @@ -0,0 +1,203 @@ +#include "src/wind_energy/actuator/motion/RigidBodyMotion.H" +#include "src/wind_energy/actuator/motion/RotorMotion.H" +#include "src/wind_energy/actuator/motion/TimeTable.H" +#include "src/utilities/constants.H" + +#include "gtest/gtest.h" + +#include +#include +#include +#include + +#include "AMReX_ParmParse.H" + +using namespace amrex::literals; + +namespace kynema_sgf_tests { +namespace { + +using kynema_sgf::actuator::motion::RigidBodyMotion; +using kynema_sgf::actuator::motion::RotorMotion; +using kynema_sgf::actuator::motion::TimeTable; +using kynema_sgf::actuator::utils::ActParser; + +void write_table(const std::string& filename, const std::string& contents) +{ + std::ofstream stream(filename); + stream << contents; +} + +TEST(ActuatorMotion, timetable_interpolation_derivative_and_integral) +{ + const std::string filename{"motion_scalar_table.txt"}; + write_table(filename, "Time Value\n0 0\n2 2\n"); + + TimeTable table; + table.read(filename, 1); + EXPECT_NEAR(table.value(1.0_rt)[0], 1.0_rt, 1.0e-14_rt); + EXPECT_NEAR(table.derivative(1.0_rt)[0], 1.0_rt, 1.0e-14_rt); + EXPECT_NEAR(table.integral(1.0_rt)[0], 0.5_rt, 1.0e-14_rt); + EXPECT_NEAR(table.value(3.0_rt)[0], 2.0_rt, 1.0e-14_rt); + EXPECT_NEAR(table.integral(3.0_rt)[0], 4.0_rt, 1.0e-14_rt); + + std::remove(filename.c_str()); +} + +TEST(ActuatorMotion, timetable_error_extrapolation) +{ + const std::string filename{"motion_error_table.txt"}; + write_table(filename, "Time Value\n0 0\n1 1\n"); + TimeTable table; + table.read(filename, 1, "error"); + EXPECT_THROW(static_cast(table.value(2.0_rt)), amrex::RuntimeError); + std::remove(filename.c_str()); +} + +TEST(ActuatorMotion, position_and_quaternion_histories) +{ + constexpr amrex::Real tol = kynema_sgf::constants::TIGHT_TOL; + const std::string position_file{"motion_position_table.txt"}; + const std::string orientation_file{"motion_orientation_table.txt"}; + write_table(position_file, "Time X Y Z\n0 0 0 0\n1 2 0 0\n"); + write_table(orientation_file, "Time Qw Qx Qy Qz\n0 1 0 0 0\n1 0 0 0 1\n"); + + amrex::ParmParse pp("MotionHistory"); + pp.add("position_timetable", position_file); + pp.add("orientation_timetable", orientation_file); + pp.add("orientation_format", "quaternion"); + RigidBodyMotion motion; + motion.read_inputs( + ActParser("UnusedMotionDefaults", "MotionHistory"), + kynema_sgf::vs::Vector::zero(), kynema_sgf::vs::Tensor::identity()); + + const auto position = motion.position(0.5_rt); + const auto velocity = motion.translation_velocity(0.5_rt); + const auto rotated = + motion.orientation(0.5_rt) & kynema_sgf::vs::Vector::ihat(); + EXPECT_NEAR(position.x(), 1.0_rt, tol); + EXPECT_NEAR(velocity.x(), 2.0_rt, tol); + EXPECT_NEAR(rotated.x(), 0.0_rt, tol); + EXPECT_NEAR(rotated.y(), 1.0_rt, tol); + + std::remove(position_file.c_str()); + std::remove(orientation_file.c_str()); +} + +TEST(ActuatorMotion, roll_pitch_yaw_and_quaternion_histories_match) +{ + constexpr amrex::Real tol = kynema_sgf::constants::TIGHT_TOL; + const std::string rpy_file{"motion_rpy_table.txt"}; + const std::string quaternion_file{"motion_quaternion_table.txt"}; + write_table(rpy_file, "Time Roll Pitch Yaw\n0 0 0 0\n1 0 0 -180\n"); + write_table(quaternion_file, "Time Qw Qx Qy Qz\n0 1 0 0 0\n1 0 0 0 1\n"); + + amrex::ParmParse rpy_pp("MotionRpy"); + rpy_pp.add("orientation_timetable", rpy_file); + RigidBodyMotion rpy_motion; + rpy_motion.read_inputs( + ActParser("UnusedRpyDefaults", "MotionRpy"), + kynema_sgf::vs::Vector::zero(), kynema_sgf::vs::Tensor::identity()); + + amrex::ParmParse quaternion_pp("MotionQuaternion"); + quaternion_pp.add("orientation_timetable", quaternion_file); + quaternion_pp.add("orientation_format", "quaternion"); + RigidBodyMotion quaternion_motion; + quaternion_motion.read_inputs( + ActParser("UnusedQuaternionDefaults", "MotionQuaternion"), + kynema_sgf::vs::Vector::zero(), kynema_sgf::vs::Tensor::identity()); + + for (const amrex::Real time : {0.0_rt, 0.25_rt, 0.5_rt, 1.0_rt}) { + const auto rpy = rpy_motion.orientation(time); + const auto quaternion = quaternion_motion.orientation(time); + for (int n = 0; n < rpy.ncomp; ++n) { + EXPECT_NEAR(rpy[n], quaternion[n], tol); + } + } + + std::remove(rpy_file.c_str()); + std::remove(quaternion_file.c_str()); +} + +TEST(ActuatorMotion, velocity_history_integrates_position) +{ + const std::string filename{"motion_velocity_table.txt"}; + write_table(filename, "Time Ux Uy Uz\n0 0 0 0\n2 2 0 0\n"); + + amrex::ParmParse pp("MotionVelocity"); + pp.addarr("center", amrex::Vector{1.0_rt, 0.0_rt, 0.0_rt}); + pp.add("velocity_timetable", filename); + RigidBodyMotion motion; + motion.read_inputs( + ActParser("UnusedVelocityDefaults", "MotionVelocity"), + {1.0_rt, 0.0_rt, 0.0_rt}, kynema_sgf::vs::Tensor::identity()); + EXPECT_NEAR(motion.position(1.0_rt).x(), 1.5_rt, 1.0e-14_rt); + EXPECT_NEAR(motion.translation_velocity(1.0_rt).x(), 1.0_rt, 1.0e-14_rt); + + std::remove(filename.c_str()); +} + +TEST(ActuatorMotion, position_and_velocity_histories_are_exclusive) +{ + const std::string position_file{"motion_exclusive_position.txt"}; + const std::string velocity_file{"motion_exclusive_velocity.txt"}; + write_table(position_file, "Time X Y Z\n0 0 0 0\n1 1 0 0\n"); + write_table(velocity_file, "Time Ux Uy Uz\n0 0 0 0\n1 1 0 0\n"); + + amrex::ParmParse pp("MotionExclusive"); + pp.add("position_timetable", position_file); + pp.add("velocity_timetable", velocity_file); + RigidBodyMotion motion; + EXPECT_THROW( + motion.read_inputs( + ActParser("UnusedExclusiveDefaults", "MotionExclusive"), + kynema_sgf::vs::Vector::zero(), kynema_sgf::vs::Tensor::identity()), + amrex::RuntimeError); + + std::remove(position_file.c_str()); + std::remove(velocity_file.c_str()); +} + +TEST(ActuatorMotion, rotor_speed_integrates_azimuth_and_holds_rate) +{ + const std::string filename{"motion_rotor_table.txt"}; + write_table(filename, "Time R1 R2\n0 0 0\n2 2 -2\n"); + + amrex::ParmParse pp("RotorHistory"); + pp.add("rotor_speed_timetable", filename); + pp.add("initial_azimuth_degrees", 90.0_rt); + RotorMotion motion; + motion.read_drone_inputs( + ActParser("UnusedRotorDefaults", "RotorHistory"), 2); + + EXPECT_NEAR(motion.omega(0, 1.0_rt), 1.0_rt, 1.0e-14_rt); + EXPECT_NEAR( + motion.azimuth(0, 1.0_rt), + 0.5_rt + 0.5_rt * std::numbers::pi_v, 1.0e-14_rt); + EXPECT_NEAR( + motion.azimuth(0, 3.0_rt), + 4.0_rt + 0.5_rt * std::numbers::pi_v, 1.0e-14_rt); + EXPECT_NEAR( + motion.azimuth(1, 3.0_rt), + -4.0_rt + 0.5_rt * std::numbers::pi_v, 1.0e-14_rt); + + std::remove(filename.c_str()); +} + +TEST(ActuatorMotion, per_rotor_initial_azimuths) +{ + amrex::ParmParse pp("RotorInitialAzimuths"); + pp.addarr("rotor_omegas", amrex::Vector{1.0_rt, -1.0_rt}); + pp.addarr( + "initial_azimuths_degrees", + amrex::Vector{0.0_rt, 180.0_rt}); + RotorMotion motion; + motion.read_drone_inputs( + ActParser("UnusedAzimuthDefaults", "RotorInitialAzimuths"), 2); + EXPECT_NEAR(motion.azimuth(0, 0.0_rt), 0.0_rt, 1.0e-14_rt); + EXPECT_NEAR( + motion.azimuth(1, 0.0_rt), std::numbers::pi_v, 1.0e-14_rt); +} + +} // namespace +} // namespace kynema_sgf_tests From e3ea78aa48ed2c5cd8b922925d78fe04ca73a365 Mon Sep 17 00:00:00 2001 From: Tony Martinez Date: Fri, 24 Jul 2026 10:58:19 -0600 Subject: [PATCH 04/16] Fixed errors from github --- src/wind_energy/actuator/drone/Drone.H | 1 - src/wind_energy/actuator/drone/drone_ops.cpp | 50 +++++++------------ .../actuator/motion/RotorMotion.cpp | 29 +++++------ .../actuator/test_actuator_motion.cpp | 29 +++++------ .../wind_energy/actuator/test_drone.cpp | 33 ++++++------ 5 files changed, 63 insertions(+), 79 deletions(-) diff --git a/src/wind_energy/actuator/drone/Drone.H b/src/wind_energy/actuator/drone/Drone.H index de2a6e1654..602cef4074 100644 --- a/src/wind_energy/actuator/drone/Drone.H +++ b/src/wind_energy/actuator/drone/Drone.H @@ -29,7 +29,6 @@ struct DroneRotor struct DroneData { int num_rotors{0}; - amrex::Real arm_length{0.0_rt}; RealList arm_lengths; RealList arm_angles_degrees; amrex::Real arm_phase_degrees{0.0_rt}; diff --git a/src/wind_energy/actuator/drone/drone_ops.cpp b/src/wind_energy/actuator/drone/drone_ops.cpp index a1a8fdcc5d..16e21861d7 100644 --- a/src/wind_energy/actuator/drone/drone_ops.cpp +++ b/src/wind_energy/actuator/drone/drone_ops.cpp @@ -65,33 +65,20 @@ void ReadInputsOp::operator()( meta.rotor_motion = std::make_shared(); meta.rotor_motion->read_drone_inputs(pp, meta.num_rotors); - // Instance inputs take precedence over Drone defaults; a scalar arm length - // expands to every rotor while the list permits asymmetric layouts. - const auto& instance = pp.params(); - const auto& defaults = pp.default_params(); - const bool instance_length = instance.contains("arm_length"); - const bool instance_lengths = instance.contains("arm_lengths"); - if (instance_length && instance_lengths) { - input_error(label, "specify arm_length or arm_lengths, not both"); + // One arm length applies to every rotor; num_rotors values permit an + // asymmetric layout. MultiParser preserves instance-over-default priority. + if (!pp.contains("arm_length")) { + input_error(label, "arm_length is required"); } - if (instance_lengths) { - pp.getarr("arm_lengths", meta.arm_lengths); - } else if (instance_length) { - pp.get("arm_length", meta.arm_length); - meta.arm_lengths.assign(meta.num_rotors, meta.arm_length); - } else if ( - defaults.contains("arm_lengths") && defaults.contains("arm_length")) { - input_error(label, "specify arm_length or arm_lengths, not both"); - } else if (defaults.contains("arm_lengths")) { - pp.getarr("arm_lengths", meta.arm_lengths); - } else if (defaults.contains("arm_length")) { - pp.get("arm_length", meta.arm_length); - meta.arm_lengths.assign(meta.num_rotors, meta.arm_length); - } else { - input_error(label, "arm_length or arm_lengths is required"); + pp.getarr("arm_length", meta.arm_lengths); + if (meta.arm_lengths.size() == 1) { + const auto arm_length = meta.arm_lengths.front(); + meta.arm_lengths.assign(meta.num_rotors, arm_length); } if (static_cast(meta.arm_lengths.size()) != meta.num_rotors) { - input_error(label, "arm_lengths must contain num_rotors values"); + input_error( + label, + "arm_length must contain either one value or num_rotors values"); } if (std::ranges::any_of(meta.arm_lengths, [](const amrex::Real length) { return length <= 0.0_rt; @@ -101,9 +88,16 @@ void ReadInputsOp::operator()( // Explicit angles describe irregular layouts. Otherwise phase rotates a // uniformly spaced layout about the body z axis. + const auto& instance = pp.params(); + const auto& defaults = pp.default_params(); const bool instance_angles = instance.contains("arm_angles_degrees"); const bool instance_phase = instance.contains("arm_phase_degrees"); - if (instance_angles && instance_phase) { + const bool instance_layout = instance_angles || instance_phase; + const bool conflicting_layout = + (instance_angles && instance_phase) || + (!instance_layout && defaults.contains("arm_angles_degrees") && + defaults.contains("arm_phase_degrees")); + if (conflicting_layout) { input_error( label, "arm_angles_degrees and arm_phase_degrees are mutually exclusive"); @@ -114,12 +108,6 @@ void ReadInputsOp::operator()( pp.get("arm_phase_degrees", meta.arm_phase_degrees); meta.arm_angles_degrees = drone::uniform_arm_angles(meta.num_rotors, meta.arm_phase_degrees); - } else if ( - defaults.contains("arm_angles_degrees") && - defaults.contains("arm_phase_degrees")) { - input_error( - label, - "arm_angles_degrees and arm_phase_degrees are mutually exclusive"); } else if (defaults.contains("arm_angles_degrees")) { pp.getarr("arm_angles_degrees", meta.arm_angles_degrees); } else { diff --git a/src/wind_energy/actuator/motion/RotorMotion.cpp b/src/wind_energy/actuator/motion/RotorMotion.cpp index 5da7b1fa8d..f94d7dcb22 100644 --- a/src/wind_energy/actuator/motion/RotorMotion.cpp +++ b/src/wind_energy/actuator/motion/RotorMotion.cpp @@ -11,28 +11,23 @@ namespace kynema_sgf::actuator::motion { void RotorMotion::read_initial_azimuths( const utils::ActParser& pp, const int num_rotors) { - const bool shared = pp.contains("initial_azimuth_degrees"); - const bool individual = pp.contains("initial_azimuths_degrees"); - if (shared && individual) { - amrex::Abort( - "initial_azimuth_degrees and initial_azimuths_degrees are mutually " - "exclusive"); - } m_initial_azimuths.assign(num_rotors, 0.0_rt); - if (shared) { - amrex::Real value = 0.0_rt; - pp.get("initial_azimuth_degrees", value); - m_initial_azimuths.assign( - num_rotors, ::kynema_sgf::utils::radians(value)); - } else if (individual) { - pp.getarr("initial_azimuths_degrees", m_initial_azimuths); - if (static_cast(m_initial_azimuths.size()) != num_rotors) { + if (pp.contains("initial_azimuth_degrees")) { + RealList input_azimuths; + pp.getarr("initial_azimuth_degrees", input_azimuths); + if (input_azimuths.size() == 1) { + const auto initial_azimuth = input_azimuths.front(); + input_azimuths.assign(num_rotors, initial_azimuth); + } + if (static_cast(input_azimuths.size()) != num_rotors) { amrex::Abort( - "initial_azimuths_degrees must contain num_rotors values"); + "initial_azimuth_degrees must contain either one value or " + "num_rotors values"); } - for (auto& value : m_initial_azimuths) { + for (auto& value : input_azimuths) { value = ::kynema_sgf::utils::radians(value); } + m_initial_azimuths = input_azimuths; } } diff --git a/unit_tests/wind_energy/actuator/test_actuator_motion.cpp b/unit_tests/wind_energy/actuator/test_actuator_motion.cpp index b62a76b6cc..5504929263 100644 --- a/unit_tests/wind_energy/actuator/test_actuator_motion.cpp +++ b/unit_tests/wind_energy/actuator/test_actuator_motion.cpp @@ -21,6 +21,7 @@ using kynema_sgf::actuator::motion::RigidBodyMotion; using kynema_sgf::actuator::motion::RotorMotion; using kynema_sgf::actuator::motion::TimeTable; using kynema_sgf::actuator::utils::ActParser; +constexpr amrex::Real test_tol = kynema_sgf::constants::TIGHT_TOL; void write_table(const std::string& filename, const std::string& contents) { @@ -35,11 +36,11 @@ TEST(ActuatorMotion, timetable_interpolation_derivative_and_integral) TimeTable table; table.read(filename, 1); - EXPECT_NEAR(table.value(1.0_rt)[0], 1.0_rt, 1.0e-14_rt); - EXPECT_NEAR(table.derivative(1.0_rt)[0], 1.0_rt, 1.0e-14_rt); - EXPECT_NEAR(table.integral(1.0_rt)[0], 0.5_rt, 1.0e-14_rt); - EXPECT_NEAR(table.value(3.0_rt)[0], 2.0_rt, 1.0e-14_rt); - EXPECT_NEAR(table.integral(3.0_rt)[0], 4.0_rt, 1.0e-14_rt); + EXPECT_NEAR(table.value(1.0_rt)[0], 1.0_rt, test_tol); + EXPECT_NEAR(table.derivative(1.0_rt)[0], 1.0_rt, test_tol); + EXPECT_NEAR(table.integral(1.0_rt)[0], 0.5_rt, test_tol); + EXPECT_NEAR(table.value(3.0_rt)[0], 2.0_rt, test_tol); + EXPECT_NEAR(table.integral(3.0_rt)[0], 4.0_rt, test_tol); std::remove(filename.c_str()); } @@ -131,8 +132,8 @@ TEST(ActuatorMotion, velocity_history_integrates_position) motion.read_inputs( ActParser("UnusedVelocityDefaults", "MotionVelocity"), {1.0_rt, 0.0_rt, 0.0_rt}, kynema_sgf::vs::Tensor::identity()); - EXPECT_NEAR(motion.position(1.0_rt).x(), 1.5_rt, 1.0e-14_rt); - EXPECT_NEAR(motion.translation_velocity(1.0_rt).x(), 1.0_rt, 1.0e-14_rt); + EXPECT_NEAR(motion.position(1.0_rt).x(), 1.5_rt, test_tol); + EXPECT_NEAR(motion.translation_velocity(1.0_rt).x(), 1.0_rt, test_tol); std::remove(filename.c_str()); } @@ -170,16 +171,16 @@ TEST(ActuatorMotion, rotor_speed_integrates_azimuth_and_holds_rate) motion.read_drone_inputs( ActParser("UnusedRotorDefaults", "RotorHistory"), 2); - EXPECT_NEAR(motion.omega(0, 1.0_rt), 1.0_rt, 1.0e-14_rt); + EXPECT_NEAR(motion.omega(0, 1.0_rt), 1.0_rt, test_tol); EXPECT_NEAR( motion.azimuth(0, 1.0_rt), - 0.5_rt + 0.5_rt * std::numbers::pi_v, 1.0e-14_rt); + 0.5_rt + 0.5_rt * std::numbers::pi_v, test_tol); EXPECT_NEAR( motion.azimuth(0, 3.0_rt), - 4.0_rt + 0.5_rt * std::numbers::pi_v, 1.0e-14_rt); + 4.0_rt + 0.5_rt * std::numbers::pi_v, test_tol); EXPECT_NEAR( motion.azimuth(1, 3.0_rt), - -4.0_rt + 0.5_rt * std::numbers::pi_v, 1.0e-14_rt); + -4.0_rt + 0.5_rt * std::numbers::pi_v, test_tol); std::remove(filename.c_str()); } @@ -189,14 +190,14 @@ TEST(ActuatorMotion, per_rotor_initial_azimuths) amrex::ParmParse pp("RotorInitialAzimuths"); pp.addarr("rotor_omegas", amrex::Vector{1.0_rt, -1.0_rt}); pp.addarr( - "initial_azimuths_degrees", + "initial_azimuth_degrees", amrex::Vector{0.0_rt, 180.0_rt}); RotorMotion motion; motion.read_drone_inputs( ActParser("UnusedAzimuthDefaults", "RotorInitialAzimuths"), 2); - EXPECT_NEAR(motion.azimuth(0, 0.0_rt), 0.0_rt, 1.0e-14_rt); + EXPECT_NEAR(motion.azimuth(0, 0.0_rt), 0.0_rt, test_tol); EXPECT_NEAR( - motion.azimuth(1, 0.0_rt), std::numbers::pi_v, 1.0e-14_rt); + motion.azimuth(1, 0.0_rt), std::numbers::pi_v, test_tol); } } // namespace diff --git a/unit_tests/wind_energy/actuator/test_drone.cpp b/unit_tests/wind_energy/actuator/test_drone.cpp index 70a34f145f..502f83d3a7 100644 --- a/unit_tests/wind_energy/actuator/test_drone.cpp +++ b/unit_tests/wind_energy/actuator/test_drone.cpp @@ -20,6 +20,7 @@ namespace { using kynema_sgf::actuator::drone::rotor_body_offsets; using kynema_sgf::actuator::drone::uniform_arm_angles; +constexpr amrex::Real test_tol = kynema_sgf::constants::TIGHT_TOL; class DroneActuatorTest : public MeshTest { @@ -75,7 +76,7 @@ class DroneActuatorTest : public MeshTest amrex::ParmParse pp_i("Actuator.D1"); pp_i.add("type", std::string("Drone")); pp_i.addarr( - "arm_lengths", + "arm_length", amrex::Vector{0.075_rt, 0.075_rt, 0.075_rt, 0.075_rt}); pp_i.addarr( "arm_angles_degrees", @@ -153,12 +154,12 @@ TEST(DroneGeometry, plus_layout) rotor_body_offsets({2.0_rt, 2.0_rt, 2.0_rt, 2.0_rt}, angles); ASSERT_EQ(offsets.size(), 4U); - EXPECT_NEAR(offsets[0].x(), 2.0_rt, 1.0e-14_rt); - EXPECT_NEAR(offsets[0].y(), 0.0_rt, 1.0e-14_rt); - EXPECT_NEAR(offsets[1].x(), 0.0_rt, 1.0e-14_rt); - EXPECT_NEAR(offsets[1].y(), 2.0_rt, 1.0e-14_rt); - EXPECT_NEAR(offsets[2].x(), -2.0_rt, 1.0e-14_rt); - EXPECT_NEAR(offsets[3].y(), -2.0_rt, 1.0e-14_rt); + EXPECT_NEAR(offsets[0].x(), 2.0_rt, test_tol); + EXPECT_NEAR(offsets[0].y(), 0.0_rt, test_tol); + EXPECT_NEAR(offsets[1].x(), 0.0_rt, test_tol); + EXPECT_NEAR(offsets[1].y(), 2.0_rt, test_tol); + EXPECT_NEAR(offsets[2].x(), -2.0_rt, test_tol); + EXPECT_NEAR(offsets[3].y(), -2.0_rt, test_tol); } TEST(DroneGeometry, unequal_irregular_arms) @@ -167,10 +168,10 @@ TEST(DroneGeometry, unequal_irregular_arms) const kynema_sgf::actuator::RealList angles{0.0_rt, 90.0_rt, 225.0_rt}; const auto offsets = rotor_body_offsets(lengths, angles); - EXPECT_NEAR(offsets[0].x(), 1.0_rt, 1.0e-14_rt); - EXPECT_NEAR(offsets[1].y(), 2.0_rt, 1.0e-14_rt); - EXPECT_NEAR(offsets[2].x(), -3.0_rt / std::sqrt(2.0_rt), 1.0e-14_rt); - EXPECT_NEAR(offsets[2].y(), -3.0_rt / std::sqrt(2.0_rt), 1.0e-14_rt); + EXPECT_NEAR(offsets[0].x(), 1.0_rt, test_tol); + EXPECT_NEAR(offsets[1].y(), 2.0_rt, test_tol); + EXPECT_NEAR(offsets[2].x(), -3.0_rt / std::sqrt(2.0_rt), test_tol); + EXPECT_NEAR(offsets[2].y(), -3.0_rt / std::sqrt(2.0_rt), test_tol); } TEST(DroneGeometry, body_orientation_defaults_to_identity) @@ -179,9 +180,9 @@ TEST(DroneGeometry, body_orientation_defaults_to_identity) kynema_sgf::vs::Vector::zero()); const auto value = rotation & kynema_sgf::vs::Vector{1.0_rt, 2.0_rt, 3.0_rt}; - EXPECT_NEAR(value.x(), 1.0_rt, 1.0e-14_rt); - EXPECT_NEAR(value.y(), 2.0_rt, 1.0e-14_rt); - EXPECT_NEAR(value.z(), 3.0_rt, 1.0e-14_rt); + EXPECT_NEAR(value.x(), 1.0_rt, test_tol); + EXPECT_NEAR(value.y(), 2.0_rt, test_tol); + EXPECT_NEAR(value.z(), 3.0_rt, test_tol); } TEST_F(DroneActuatorTest, composite_lifecycle) @@ -199,10 +200,10 @@ TEST_F(DroneActuatorTest, composite_lifecycle) ASSERT_EQ(drone->meta().rotors.size(), 4U); EXPECT_NEAR( drone->meta().rotors[0]->data.meta().center.x(), - 0.075_rt / std::sqrt(2.0_rt), 1.0e-14_rt); + 0.075_rt / std::sqrt(2.0_rt), test_tol); EXPECT_NEAR( drone->meta().rotors[0]->data.meta().center.y(), - 0.075_rt / std::sqrt(2.0_rt), 1.0e-14_rt); + 0.075_rt / std::sqrt(2.0_rt), test_tol); actuator.post_init_actions(); actuator.pre_advance_work(); From a2659252721b2c27936ab65ab29a8768725bc924 Mon Sep 17 00:00:00 2001 From: Tony Martinez Date: Fri, 24 Jul 2026 12:46:56 -0600 Subject: [PATCH 05/16] updated documentation and outputs --- docs/sphinx/spelling-wordlist.txt | 7 +- docs/sphinx/user/features.rst | 3 +- docs/sphinx/user/inputs_Actuator.rst | 286 +++++++++++++++++- src/wind_energy/actuator/drone/drone_ops.H | 2 +- src/wind_energy/actuator/drone/drone_ops.cpp | 128 +++++--- .../actuator/sector/actuator_sector_ops.H | 24 ++ .../actuator/sector/actuator_sector_ops.cpp | 97 ++++-- .../wind_energy/actuator/test_drone.cpp | 43 ++- 8 files changed, 516 insertions(+), 74 deletions(-) diff --git a/docs/sphinx/spelling-wordlist.txt b/docs/sphinx/spelling-wordlist.txt index 1c9a93264e..2824062a2d 100644 --- a/docs/sphinx/spelling-wordlist.txt +++ b/docs/sphinx/spelling-wordlist.txt @@ -116,6 +116,9 @@ Prandtl pre priori proj +quadcopter +quaternion +quaternions reStructuredText regrid regridding @@ -139,6 +142,8 @@ teardown timestep timestepping timesteps +timetable +timetables tke vel vof @@ -153,4 +158,4 @@ ylo yhi yt zlo -zhi \ No newline at end of file +zhi diff --git a/docs/sphinx/user/features.rst b/docs/sphinx/user/features.rst index 8ae7c69fc1..241d6c9c77 100644 --- a/docs/sphinx/user/features.rst +++ b/docs/sphinx/user/features.rst @@ -76,6 +76,8 @@ Flow physics * Actuator turbine representations: Joukowsky disks, uniform disks, actuator line [:ref:`doc `, :ref:`inp `] + * Prescribed-motion rotor and multi-rotor Drone actuator representations [:ref:`inp `] + * Coupling with OpenFAST * Coupling with Nalu-Wind for blade resolved simulations @@ -164,4 +166,3 @@ High performance computing * native AMReX solvers such as MLMG [:ref:`inp `] * hypre - diff --git a/docs/sphinx/user/inputs_Actuator.rst b/docs/sphinx/user/inputs_Actuator.rst index 9d00698fa5..63d292d11d 100644 --- a/docs/sphinx/user/inputs_Actuator.rst +++ b/docs/sphinx/user/inputs_Actuator.rst @@ -27,7 +27,7 @@ turbines as actuator disks and actuator line models. This string identifies the type of actuator to use. The ones currently supported are: ``UniformCtDisk``, ``JoukowskyDisk``, ``TurbineFastLine``, ``TurbineFastDisk``, ``TurbineKynemaLine``, ``FixedWingLine``, and - ``ActuatorSector``. + ``ActuatorSector``, and ``Drone``. It is recommended to group common parameters across actuators using the ``Actuator.[type].[param]``. For example:: @@ -272,6 +272,111 @@ Example for ``FixedWingLine``:: This input allows actuator force coordinate directions to be deactivated by specifying a 0.0 in for the x, y, or z component of this vector. +Prescribed actuator motion +"""""""""""""""""""""""""" + +``ActuatorSector`` and ``Drone`` support constant or time-dependent rigid-body +motion. Timetable files are whitespace-delimited text files. The first line is +a header and each subsequent line begins with time in seconds. Times must be +strictly increasing. Values are linearly interpolated; velocity histories are +integrated with the corresponding piecewise-linear, trapezoidal rule. + +Translation must be specified using no more than one of +``translation_velocity``, ``velocity_timetable``, and ``position_timetable``. +When a velocity source is used, ``center`` is the position at time zero. +``center`` must not be combined with ``position_timetable`` because the +position history supplies the complete position. + +Similarly, orientation must be specified using no more than one of a constant +angular velocity, ``angular_velocity_timetable``, and +``orientation_timetable``. Angular velocities are expressed in the global CFD +frame; this is currently the only supported angular-velocity frame. + +For example, a position history has three coordinates in meters:: + + Time X Y Z + 0.0 0.0 0.0 1.0 + 1.0 0.0 0.0 2.0 + +The recommended orientation format is roll, pitch, and yaw in degrees:: + + Time Roll Pitch Yaw + 0.0 0.0 0.0 0.0 + 1.0 0.0 5.0 10.0 + +These angles define the body-to-CFD rotation +``Rz(yaw) Ry(pitch) Rx(roll)``. A quaternion history is also supported using +scalar-first ``W X Y Z`` unit quaternions. Input quaternions are normalized +before they are interpolated. Both orientation formats are interpolated with +quaternion spherical linear interpolation. + +.. input_param:: Actuator.ActuatorSector.position_timetable + + **type:** String, optional + + File containing ``Time X Y Z``, where position is in meters in the CFD + frame. + +.. input_param:: Actuator.ActuatorSector.velocity_timetable + + **type:** String, optional + + File containing ``Time U V W``, where translational velocity is in m/s in + the CFD frame. ``center`` is required and supplies the initial position. + +.. input_param:: Actuator.ActuatorSector.orientation_timetable + + **type:** String, optional + + File containing either ``Time Roll Pitch Yaw`` in degrees or + ``Time W X Y Z``, as selected by ``orientation_format``. + +.. input_param:: Actuator.ActuatorSector.orientation_format + + **type:** String, optional, default = ``roll_pitch_yaw`` + + Format of ``orientation_timetable``. Valid values are + ``roll_pitch_yaw`` and ``quaternion``. + +.. input_param:: Actuator.ActuatorSector.angular_velocity_timetable + + **type:** String, optional + + File containing ``Time OmegaX OmegaY OmegaZ``, where angular velocity is in + rad/s in the global CFD frame. + +.. input_param:: Actuator.ActuatorSector.angular_velocity_frame + + **type:** String, optional, default = ``global`` + + Coordinate frame for constant and tabulated body angular velocity. Only + ``global`` is currently supported. + +.. input_param:: Actuator.ActuatorSector.rotor_speed_timetable + + **type:** String, optional + + Rotor-speed history in rad/s. For ``ActuatorSector`` the file contains + ``Time Omega`` and replaces ``omega``. For ``Drone`` it contains one + angular-speed column per rotor and replaces ``rotor_omegas``. Rotor azimuth + is obtained by integrating this history. + +.. input_param:: Actuator.ActuatorSector.initial_azimuth_degrees + + **type:** Real number, optional, default = 0.0 + + Initial blade azimuth in degrees. A Drone accepts either one value shared by + every rotor or ``num_rotors`` values. + +.. input_param:: Actuator.ActuatorSector.timetable_extrapolation + + **type:** String, optional, default = ``hold`` + + Behavior outside a timetable's time range. ``hold`` holds the nearest + endpoint value and prints a warning the first time each bound is crossed. + ``error`` terminates the simulation instead. + + ActuatorSector """""""""""""" @@ -516,6 +621,185 @@ Example for ``ActuatorSector``:: airfoil coefficients. +Drone +""""" + +The ``Drone`` actuator is a rigid-body composite of ``ActuatorSector`` rotors. +All rotors share the drone position and orientation, while each rotor retains +its own signed speed, azimuth, hub location, and blade-mirroring choice. Rotor +aerodynamic and force-projection inputs use the ``Actuator.Drone`` namespace +and are passed directly to every child sector. + +The drone center is the origin of its body frame. Arms lie in the body x-y +plane, and each rotor axis initially points along body +z. Viewed from body +z, +positive arm angles proceed from body +x toward body +y:: + + body +y + ^ + | + R2 | R1 + \ | / + body -x <-----+-----> body +x + / | \ + R3 | R4 + | + + body +z: out of page + +The sketch shows a four-rotor X layout with ``arm_phase_degrees = 45``. +``R1`` is placed at the phase angle and subsequent rotors are placed at +increasing, uniformly spaced angles. A positive arm angle therefore denotes a +counterclockwise geometric placement in this view; it does not specify rotor +spin direction. + +Example for a stationary X-layout quadcopter:: + + incflo.physics = FreeStream Actuator + ICNS.source_terms = ActuatorForcing + Actuator.labels = Q1 + Actuator.Q1.type = Drone + + Actuator.Drone.num_rotors = 4 + Actuator.Drone.arm_length = 0.075 + Actuator.Drone.arm_phase_degrees = 45.0 + Actuator.Drone.rotor_omegas = 2500.0 -2500.0 2500.0 -2500.0 + Actuator.Drone.initial_azimuth_degrees = 0.0 90.0 180.0 270.0 + Actuator.Drone.mirror_blades = false true false true + + Actuator.Q1.center = 0.0 0.0 1.0 + Actuator.Q1.body_orientation_degrees = 0.0 0.0 0.0 + + Actuator.Drone.rotor_diameter = 0.10 + Actuator.Drone.root_radius_fraction = 0.18 + Actuator.Drone.num_blades = 2 + Actuator.Drone.epsilon_chord = 0.5 + Actuator.Drone.epsilon_dr = 1.5 + Actuator.Drone.epsilon_dl = 1.5 + Actuator.Drone.min_chord_dr = 2.0 + Actuator.Drone.span_locs = 0.0 1.0 + Actuator.Drone.chord = 0.01 0.006 + Actuator.Drone.twist = 12.0 4.0 + Actuator.Drone.airfoil_table = airfoil.txt + Actuator.Drone.airfoil_type = openfast + +For unequal arm lengths, replace the scalar with one value per rotor:: + + Actuator.Drone.arm_length = 0.075 0.080 0.075 0.080 + +For an irregular layout, provide one base angle per rotor. The phase can then +rotate the complete pattern without changing its relative geometry:: + + Actuator.Drone.arm_angles_degrees = 0.0 90.0 190.0 295.0 + Actuator.Drone.arm_phase_degrees = 20.0 + +.. input_param:: Actuator.Drone.num_rotors + + **type:** int, mandatory + + Number of rotors. The value must be positive. + +.. input_param:: Actuator.Drone.arm_length + + **type:** Real number or list of real numbers, mandatory + + Body-center to rotor-hub distance in meters. One value is applied to every + rotor. Alternatively, provide exactly ``num_rotors`` positive values in + rotor order. + +.. input_param:: Actuator.Drone.arm_phase_degrees + + **type:** Real number, optional, default = 0.0 + + Body-frame angle of the first rotor, measured in degrees from body +x toward + body +y when the default uniform pattern is used. More generally, this + value is added to every base angle in ``arm_angles_degrees``, rotating the + complete arm pattern without changing the relative angles. + +.. input_param:: Actuator.Drone.arm_angles_degrees + + **type:** List of real numbers, optional + + Base body-frame arm angles in degrees. The list must contain ``num_rotors`` + distinct angles. If omitted, the base pattern is uniformly spaced at + ``0, 360 / num_rotors, ...``. ``arm_phase_degrees`` is added to every base + angle. + +.. input_param:: Actuator.Drone.center + + **type:** List of 3 real numbers, conditionally mandatory + + Initial body-center position in meters in the CFD frame. It is required + unless ``position_timetable`` supplies the complete position history, and + cannot be combined with that history. + +.. input_param:: Actuator.Drone.body_orientation_degrees + + **type:** List of 3 real numbers, optional, default = 0.0 0.0 0.0 + + Initial body-to-CFD roll, pitch, and yaw angles in degrees. The rotations + are composed as ``Rz(yaw) Ry(pitch) Rx(roll)``. This input cannot be combined + with ``orientation_timetable``. + +.. input_param:: Actuator.Drone.translation_velocity + + **type:** List of 3 real numbers, optional, default = 0.0 0.0 0.0 + + Constant body translational velocity in m/s in the CFD frame. This input is + mutually exclusive with ``position_timetable`` and ``velocity_timetable``. + +.. input_param:: Actuator.Drone.angular_velocity + + **type:** List of 3 real numbers, optional, default = 0.0 0.0 0.0 + + Constant body angular velocity in rad/s in the global CFD frame. This input + is mutually exclusive with ``orientation_timetable`` and + ``angular_velocity_timetable``. + +.. input_param:: Actuator.Drone.rotor_omegas + + **type:** List of real numbers, conditionally mandatory + + Signed rotor angular speeds in rad/s, with exactly one value per rotor. + Specify exactly one of ``rotor_omegas`` and ``rotor_speed_timetable``. + Positive values follow the clockwise-positive blade convention used by + ``ActuatorSector``. + +.. input_param:: Actuator.Drone.initial_azimuth_degrees + + **type:** Real number or list of real numbers, optional, default = 0.0 + + Initial blade azimuth in degrees. One value is applied to every rotor, or + exactly ``num_rotors`` values may be supplied. + +.. input_param:: Actuator.Drone.mirror_blades + + **type:** List of booleans, optional, default = ``false`` for every rotor + + One value per rotor. A mirrored rotor reverses the sign of its blade twist. + This permits geometrically paired counter-rotating propellers without a + separate rotation-direction input. + +The Drone accepts the prescribed-motion parameters described above using +either the ``Actuator.Drone`` defaults or the individual +``Actuator.