2016-04-21 06:56:44 -03:00
|
|
|
/*
|
|
|
|
This program is free software: you can redistribute it and/or modify
|
|
|
|
it under the terms of the GNU General Public License as published by
|
|
|
|
the Free Software Foundation, either version 3 of the License, or
|
|
|
|
(at your option) any later version.
|
|
|
|
|
|
|
|
This program is distributed in the hope that it will be useful,
|
|
|
|
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
|
|
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
|
|
GNU General Public License for more details.
|
|
|
|
|
|
|
|
You should have received a copy of the GNU General Public License
|
|
|
|
along with this program. If not, see <http://www.gnu.org/licenses/>.
|
|
|
|
*/
|
|
|
|
/*
|
|
|
|
simple electric motor simulator class
|
|
|
|
*/
|
|
|
|
|
|
|
|
#include "SIM_Motor.h"
|
|
|
|
#include <AP_Motors/AP_Motors.h>
|
|
|
|
|
|
|
|
using namespace SITL;
|
|
|
|
|
|
|
|
// calculate rotational accel and thrust for a motor
|
2018-07-31 09:33:23 -03:00
|
|
|
void Motor::calculate_forces(const struct sitl_input &input,
|
2016-04-21 06:56:44 -03:00
|
|
|
uint8_t motor_offset,
|
2022-04-13 11:03:53 -03:00
|
|
|
Vector3f &torque,
|
2020-10-20 04:19:54 -03:00
|
|
|
Vector3f &thrust,
|
|
|
|
const Vector3f &velocity_air_bf,
|
2022-04-20 14:53:11 -03:00
|
|
|
const Vector3f &gyro,
|
2020-10-20 04:19:54 -03:00
|
|
|
float air_density,
|
2022-04-22 20:30:47 -03:00
|
|
|
float voltage,
|
|
|
|
bool use_drag)
|
2016-04-21 06:56:44 -03:00
|
|
|
{
|
2016-05-01 05:06:49 -03:00
|
|
|
// fudge factors
|
2021-06-07 23:47:30 -03:00
|
|
|
const float yaw_scale = radians(40);
|
2016-05-01 05:06:49 -03:00
|
|
|
|
2020-10-20 04:19:54 -03:00
|
|
|
const float pwm = input.servos[motor_offset+servo];
|
|
|
|
float command = pwm_to_command(pwm);
|
|
|
|
float voltage_scale = voltage / voltage_max;
|
|
|
|
|
2020-10-24 17:38:02 -03:00
|
|
|
if (voltage_scale < 0.1) {
|
|
|
|
// battery is dead
|
2022-04-13 11:03:53 -03:00
|
|
|
torque.zero();
|
2020-10-24 17:38:02 -03:00
|
|
|
thrust.zero();
|
|
|
|
current = 0;
|
|
|
|
return;
|
|
|
|
}
|
|
|
|
|
2020-10-20 04:19:54 -03:00
|
|
|
// apply slew limiter to command
|
|
|
|
uint64_t now_us = AP_HAL::micros64();
|
|
|
|
if (last_calc_us != 0 && slew_max > 0) {
|
|
|
|
float dt = (now_us - last_calc_us)*1.0e-6;
|
|
|
|
float slew_max_change = slew_max * dt;
|
|
|
|
command = constrain_float(command, last_command-slew_max_change, last_command+slew_max_change);
|
|
|
|
}
|
|
|
|
last_calc_us = now_us;
|
|
|
|
last_command = command;
|
2016-05-01 05:06:49 -03:00
|
|
|
|
|
|
|
// the yaw torque of the motor
|
2022-04-13 16:58:16 -03:00
|
|
|
Vector3f rotor_torque = thrust_vector * yaw_factor * command * yaw_scale * voltage_scale * -1.0;
|
2020-10-20 04:19:54 -03:00
|
|
|
|
2022-04-20 14:53:11 -03:00
|
|
|
// velocity of motor through air
|
|
|
|
Vector3f motor_vel = velocity_air_bf;
|
|
|
|
|
|
|
|
// add velocity of motor about center due to vehicle rotation
|
|
|
|
motor_vel += -(position % gyro);
|
|
|
|
|
|
|
|
// calculate velocity into prop, clipping at zero
|
|
|
|
float velocity_in = MAX(0, -motor_vel.projected(thrust_vector).z);
|
2016-05-01 05:06:49 -03:00
|
|
|
|
|
|
|
// get thrust for untilted motor
|
2022-04-13 09:50:36 -03:00
|
|
|
float motor_thrust = calc_thrust(command, air_density, velocity_in, voltage_scale);
|
2020-10-20 04:19:54 -03:00
|
|
|
|
2022-04-20 14:53:11 -03:00
|
|
|
// thrust in bodyframe NED
|
2022-04-13 16:58:16 -03:00
|
|
|
thrust = thrust_vector * motor_thrust;
|
2016-05-01 05:06:49 -03:00
|
|
|
|
|
|
|
// work out roll and pitch of motor relative to it pointing straight up
|
|
|
|
float roll = 0, pitch = 0;
|
2016-07-03 09:02:09 -03:00
|
|
|
|
|
|
|
uint64_t now = AP_HAL::micros64();
|
2016-05-01 05:06:49 -03:00
|
|
|
|
|
|
|
// possibly roll and/or pitch the motor
|
2016-04-21 07:29:49 -03:00
|
|
|
if (roll_servo >= 0) {
|
2016-07-03 09:02:09 -03:00
|
|
|
uint16_t servoval = update_servo(input.servos[roll_servo+motor_offset], now, last_roll_value);
|
2016-04-21 20:55:22 -03:00
|
|
|
if (roll_min < roll_max) {
|
2016-04-21 21:26:37 -03:00
|
|
|
roll = constrain_float(roll_min + (servoval-1000)*0.001*(roll_max-roll_min), roll_min, roll_max);
|
2016-04-21 20:55:22 -03:00
|
|
|
} else {
|
2016-04-21 21:26:37 -03:00
|
|
|
roll = constrain_float(roll_max + (2000-servoval)*0.001*(roll_min-roll_max), roll_max, roll_min);
|
2016-04-21 20:55:22 -03:00
|
|
|
}
|
2016-04-21 07:29:49 -03:00
|
|
|
}
|
|
|
|
if (pitch_servo >= 0) {
|
2016-07-03 09:02:09 -03:00
|
|
|
uint16_t servoval = update_servo(input.servos[pitch_servo+motor_offset], now, last_pitch_value);
|
2016-04-21 20:55:22 -03:00
|
|
|
if (pitch_min < pitch_max) {
|
2016-04-21 21:26:37 -03:00
|
|
|
pitch = constrain_float(pitch_min + (servoval-1000)*0.001*(pitch_max-pitch_min), pitch_min, pitch_max);
|
2016-04-21 20:55:22 -03:00
|
|
|
} else {
|
2016-04-21 21:26:37 -03:00
|
|
|
pitch = constrain_float(pitch_max + (2000-servoval)*0.001*(pitch_min-pitch_max), pitch_max, pitch_min);
|
2016-04-21 20:55:22 -03:00
|
|
|
}
|
2016-05-01 05:06:49 -03:00
|
|
|
}
|
2016-07-03 09:02:09 -03:00
|
|
|
last_change_usec = now;
|
2016-05-01 05:06:49 -03:00
|
|
|
|
2021-06-07 23:47:30 -03:00
|
|
|
// calculate torque in newton-meters
|
2022-04-13 16:58:16 -03:00
|
|
|
torque = (position % thrust) + rotor_torque;
|
2021-06-07 23:47:30 -03:00
|
|
|
|
2016-05-01 05:06:49 -03:00
|
|
|
// possibly rotate the thrust vector and the rotor torque
|
|
|
|
if (!is_zero(roll) || !is_zero(pitch)) {
|
2016-04-21 07:29:49 -03:00
|
|
|
Matrix3f rotation;
|
2016-05-01 05:06:49 -03:00
|
|
|
rotation.from_euler(radians(roll), radians(pitch), 0);
|
2016-04-21 07:29:49 -03:00
|
|
|
thrust = rotation * thrust;
|
2021-06-07 23:47:30 -03:00
|
|
|
torque = rotation * torque;
|
2016-04-21 07:29:49 -03:00
|
|
|
}
|
2016-05-01 05:06:49 -03:00
|
|
|
|
2022-04-22 20:30:47 -03:00
|
|
|
if (use_drag) {
|
|
|
|
// calculate momentum drag per motor
|
|
|
|
const float momentum_drag_factor = momentum_drag_coefficient * sqrtf(air_density * true_prop_area);
|
|
|
|
Vector3f momentum_drag;
|
|
|
|
momentum_drag.x = momentum_drag_factor * motor_vel.x * (sqrtf(fabsf(thrust.y)) + sqrtf(fabsf(thrust.z)));
|
|
|
|
momentum_drag.y = momentum_drag_factor * motor_vel.y * (sqrtf(fabsf(thrust.x)) + sqrtf(fabsf(thrust.z)));
|
|
|
|
// The application of momentum drag to the Z axis is a 'hack' to compensate for incorrect modelling
|
|
|
|
// of the variation of thust with inflow velocity. If not applied, the vehicle will
|
|
|
|
// climb at an unrealistic rate during operation in STABILIZE. TODO replace prop and motor model in
|
|
|
|
// with one based on DC motor, momentum disc and blade element theory.
|
|
|
|
momentum_drag.z = momentum_drag_factor * motor_vel.z * (sqrtf(fabsf(thrust.x)) + sqrtf(fabsf(thrust.y)) + sqrtf(fabsf(thrust.z)));
|
|
|
|
|
|
|
|
thrust -= momentum_drag;
|
|
|
|
}
|
|
|
|
|
2020-10-20 04:19:54 -03:00
|
|
|
// calculate current
|
|
|
|
float power = power_factor * fabsf(motor_thrust);
|
|
|
|
current = power / MAX(voltage, 0.1);
|
2016-04-21 06:56:44 -03:00
|
|
|
}
|
|
|
|
|
2016-07-03 09:02:09 -03:00
|
|
|
/*
|
|
|
|
update and return current value of a servo. Calculated as 1000..2000
|
|
|
|
*/
|
2021-02-01 12:37:57 -04:00
|
|
|
uint16_t Motor::update_servo(uint16_t demand, uint64_t time_usec, float &last_value) const
|
2016-07-03 09:02:09 -03:00
|
|
|
{
|
|
|
|
if (servo_rate <= 0) {
|
|
|
|
return demand;
|
|
|
|
}
|
|
|
|
if (servo_type == SERVO_RETRACT) {
|
|
|
|
// handle retract servos
|
|
|
|
if (demand > 1700) {
|
|
|
|
demand = 2000;
|
|
|
|
} else if (demand < 1300) {
|
|
|
|
demand = 1000;
|
|
|
|
} else {
|
|
|
|
demand = last_value;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
demand = constrain_int16(demand, 1000, 2000);
|
|
|
|
float dt = (time_usec - last_change_usec) * 1.0e-6f;
|
|
|
|
// assume servo moves through 90 degrees over 1000 to 2000
|
|
|
|
float max_change = 1000 * (dt / servo_rate) * 60.0f / 90.0f;
|
|
|
|
last_value = constrain_float(demand, last_value-max_change, last_value+max_change);
|
|
|
|
return uint16_t(last_value+0.5);
|
|
|
|
}
|
2019-01-12 23:07:36 -04:00
|
|
|
|
|
|
|
|
|
|
|
// calculate current and voltage
|
2020-10-20 04:19:54 -03:00
|
|
|
float Motor::get_current(void) const
|
|
|
|
{
|
|
|
|
return current;
|
|
|
|
}
|
|
|
|
|
|
|
|
// setup PWM ranges for this motor
|
|
|
|
void Motor::setup_params(uint16_t _pwm_min, uint16_t _pwm_max, float _spin_min, float _spin_max, float _expo, float _slew_max,
|
2022-04-13 11:03:53 -03:00
|
|
|
float _diagonal_size, float _power_factor, float _voltage_max, float _effective_prop_area,
|
2022-04-22 20:30:47 -03:00
|
|
|
float _velocity_max, Vector3f _position, Vector3f _thrust_vector, float _yaw_factor,
|
|
|
|
float _true_prop_area, float _momentum_drag_coefficient)
|
2019-01-12 23:07:36 -04:00
|
|
|
{
|
2020-10-20 04:19:54 -03:00
|
|
|
mot_pwm_min = _pwm_min;
|
|
|
|
mot_pwm_max = _pwm_max;
|
|
|
|
mot_spin_min = _spin_min;
|
|
|
|
mot_spin_max = _spin_max;
|
|
|
|
mot_expo = _expo;
|
|
|
|
slew_max = _slew_max;
|
|
|
|
power_factor = _power_factor;
|
|
|
|
voltage_max = _voltage_max;
|
2022-04-13 09:50:36 -03:00
|
|
|
effective_prop_area = _effective_prop_area;
|
|
|
|
max_outflow_velocity = _velocity_max;
|
2022-04-22 20:30:47 -03:00
|
|
|
true_prop_area = _true_prop_area;
|
|
|
|
momentum_drag_coefficient = _momentum_drag_coefficient;
|
2022-04-13 16:58:16 -03:00
|
|
|
|
2022-04-13 18:15:57 -03:00
|
|
|
if (!_position.is_zero()) {
|
|
|
|
position = _position;
|
|
|
|
} else {
|
|
|
|
position.x = cosf(radians(angle)) * _diagonal_size;
|
|
|
|
position.y = sinf(radians(angle)) * _diagonal_size;
|
|
|
|
position.z = 0;
|
|
|
|
}
|
|
|
|
|
|
|
|
if (!_thrust_vector.is_zero()) {
|
|
|
|
thrust_vector = _thrust_vector;
|
|
|
|
}
|
2022-04-20 14:52:31 -03:00
|
|
|
if (!is_zero(_yaw_factor)) {
|
|
|
|
yaw_factor = _yaw_factor;
|
|
|
|
}
|
2020-10-20 04:19:54 -03:00
|
|
|
}
|
2019-01-12 23:07:36 -04:00
|
|
|
|
2020-10-20 04:19:54 -03:00
|
|
|
/*
|
|
|
|
convert a PWM value to a command value from 0 to 1
|
|
|
|
*/
|
|
|
|
float Motor::pwm_to_command(float pwm) const
|
|
|
|
{
|
|
|
|
const float pwm_thrust_max = mot_pwm_min + mot_spin_max * (mot_pwm_max - mot_pwm_min);
|
|
|
|
const float pwm_thrust_min = mot_pwm_min + mot_spin_min * (mot_pwm_max - mot_pwm_min);
|
|
|
|
const float pwm_thrust_range = pwm_thrust_max - pwm_thrust_min;
|
|
|
|
return constrain_float((pwm-pwm_thrust_min)/pwm_thrust_range, 0, 1);
|
|
|
|
}
|
2019-01-12 23:07:36 -04:00
|
|
|
|
2020-10-20 04:19:54 -03:00
|
|
|
/*
|
|
|
|
calculate thrust given a command value
|
|
|
|
*/
|
2022-04-13 09:50:36 -03:00
|
|
|
float Motor::calc_thrust(float command, float air_density, float velocity_in, float voltage_scale) const
|
2020-10-20 04:19:54 -03:00
|
|
|
{
|
2022-04-13 09:50:36 -03:00
|
|
|
float velocity_out = voltage_scale * max_outflow_velocity * sqrtf((1-mot_expo)*command + mot_expo*sq(command));
|
2020-10-20 04:19:54 -03:00
|
|
|
float ret = 0.5 * air_density * effective_prop_area * (sq(velocity_out) - sq(velocity_in));
|
|
|
|
#if 0
|
|
|
|
if (command > 0) {
|
|
|
|
::printf("air_density=%f effective_prop_area=%f velocity_in=%f velocity_max=%f\n",
|
2022-04-13 09:50:36 -03:00
|
|
|
air_density, effective_prop_area, velocity_in, voltage_scale * max_outflow_velocity);
|
2020-10-20 04:19:54 -03:00
|
|
|
::printf("calc_thrust %.3f %.3f\n", command, ret);
|
2019-03-01 07:06:47 -04:00
|
|
|
}
|
2020-10-20 04:19:54 -03:00
|
|
|
#endif
|
|
|
|
return ret;
|
2019-01-12 23:07:36 -04:00
|
|
|
}
|