Line data Source code
1 : /*
2 : This program is free software: you can redistribute it and/or modify
3 : it under the terms of the GNU General Public License as published by
4 : the Free Software Foundation, either version 3 of the License, or
5 : (at your option) any later version.
6 :
7 : This program is distributed in the hope that it will be useful,
8 : but WITHOUT ANY WARRANTY; without even the implied warranty of
9 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
10 : GNU General Public License for more details.
11 :
12 : You should have received a copy of the GNU General Public License
13 : along with this program. If not, see <http://www.gnu.org/licenses/>.
14 : */
15 : /*
16 : Submarine simulator class
17 : */
18 :
19 : #include "SIM_Submarine.h"
20 : #include <AP_Motors/AP_Motors.h>
21 :
22 : #include <stdio.h>
23 :
24 : using namespace SITL;
25 :
26 : static Thruster vectored_thrusters[] =
27 : { // Motor # Roll Factor Pitch Factor Yaw Factor Throttle Factor Forward Factor Lateral Factor
28 : Thruster(0, 0, 0, 1.0f, 0, -1.0f, 1.0f),
29 : Thruster(1, 0, 0, -1.0f, 0, -1.0f, -1.0f),
30 : Thruster(2, 0, 0, -1.0f, 0, 1.0f, 1.0f),
31 : Thruster(3, 0, 0, 1.0f, 0, 1.0f, -1.0f),
32 : Thruster(4, 1.0f, 0, 0, -1.0f, 0, 0),
33 : Thruster(5, -1.0f, 0, 0, -1.0f, 0, 0)
34 : };
35 :
36 :
37 : static Thruster vectored_6dof_thrusters[] =
38 : {
39 : // Motor # Roll Factor Pitch Factor Yaw Factor Throttle Factor Forward Factor Lateral Factor
40 : Thruster(0, 0, 0, 1.0f, 0, -1.0f, 1.0f),
41 : Thruster(1, 0, 0, -1.0f, 0, -1.0f, -1.0f),
42 : Thruster(2, 0, 0, -1.0f, 0, 1.0f, 1.0f),
43 : Thruster(3, 0, 0, 1.0f, 0, 1.0f, -1.0f),
44 : Thruster(4, 1.0f, -1.0f, 0, -1.0f, 0, 0),
45 : Thruster(5, -1.0f, -1.0f, 0, -1.0f, 0, 0),
46 : Thruster(6, 1.0f, 1.0f, 0, -1.0f, 0, 0),
47 : Thruster(7, -1.0f, 1.0f, 0, -1.0f, 0, 0)
48 : };
49 :
50 83 : Submarine::Submarine(const char *frame_str) :
51 : Aircraft(frame_str),
52 83 : frame(NULL)
53 : {
54 83 : frame_height = 0.0;
55 83 : ground_behavior = GROUND_BEHAVIOR_NONE;
56 :
57 : // default to vectored frame
58 83 : thrusters = vectored_thrusters;
59 83 : n_thrusters = 6;
60 :
61 83 : if (strstr(frame_str, "vectored_6dof")) {
62 1 : thrusters = vectored_6dof_thrusters;
63 1 : n_thrusters = 8;
64 : }
65 83 : lock_step_scheduled = true;
66 :
67 83 : constexpr float default_battery_resistance_ohm = 0.033;
68 166 : battery.setup(sitl->batt_capacity_ah,
69 : default_battery_resistance_ohm,
70 83 : sitl->batt_voltage,
71 : ambient_outside_temperature_degC());
72 83 : }
73 :
74 65467018 : float Submarine::perpendicular_distance_to_rangefinder_surface() const
75 : {
76 65467018 : const float floor_depth = calculate_sea_floor_depth(position);
77 65467018 : return floor_depth - position.z;
78 : }
79 :
80 : // calculate rotational and linear accelerations
81 6753197 : void Submarine::calculate_forces(const struct sitl_input &input, Vector3f &rot_accel, Vector3f &body_accel)
82 : {
83 6753197 : rot_accel = Vector3f(0,0,0);
84 :
85 : // slight positive buoyancy
86 6753197 : body_accel = dcm.transposed() * Vector3f(0, 0, -calculate_buoyancy_acceleration());
87 :
88 47450111 : for (int i = 0; i < n_thrusters; i++) {
89 40696914 : Thruster t = thrusters[i];
90 40696914 : int16_t pwm = input.servos[t.servo];
91 40696914 : float output = 0;
92 : // if valid pwm and not in the esc deadzone
93 : // TODO: extract deadzone from parameters/vehicle code
94 40696914 : if (!battery_is_empty() && pwm < 2000 && pwm > 1000 && (pwm < 1475 || pwm > 1525)) {
95 16683433 : output = (pwm - 1500) / 400.0; // range -1~1
96 : }
97 :
98 40696914 : float thrust = output * fabs(output) * frame_property.thrust; // approximate pwm to thrust function using a quadratic curve
99 40696914 : body_accel += t.linear * thrust / frame_property.weight;
100 40696914 : rot_accel += t.rotational * thrust * frame_property.thruster_mount_radius / frame_property.moment_of_inertia;
101 : }
102 :
103 : // Limit movement at the sea floor
104 6753197 : const float floor_depth = calculate_sea_floor_depth(position);
105 6753197 : if (position.z > floor_depth && body_accel.z > -GRAVITY_MSS) {
106 0 : body_accel.z = -GRAVITY_MSS;
107 : }
108 :
109 : // Calculate linear drag forces
110 6753197 : Vector3f linear_drag_forces;
111 6753197 : calculate_drag_force(velocity_air_bf, frame_property.linear_drag_coefficient, linear_drag_forces);
112 : // Add forces in body frame accel
113 6753197 : body_accel -= linear_drag_forces / frame_property.weight;
114 :
115 : // Calculate angular drag forces
116 : // TODO: This results in the wrong units. Fix the math.
117 6753197 : Vector3f angular_drag_torque;
118 6753197 : calculate_angular_drag_torque(gyro, frame_property.angular_drag_coefficient, angular_drag_torque);
119 :
120 : // Calculate torque induced by buoyancy foams on the frame
121 6753197 : Vector3f buoyancy_torque;
122 6753197 : calculate_buoyancy_torque(buoyancy_torque);
123 : // Add forces in body frame accel
124 6753197 : rot_accel -= angular_drag_torque / frame_property.moment_of_inertia;
125 6753197 : rot_accel += buoyancy_torque / frame_property.moment_of_inertia;
126 6753197 : add_shove_forces(rot_accel, body_accel);
127 6753197 : }
128 :
129 6753197 : void Submarine::update_battery(const struct sitl_input &input)
130 : {
131 6753197 : battery.maybe_reset(sitl->batt_voltage, sitl->batt_capacity_ah);
132 6753197 : battery_current = 0.0f;
133 6753197 : constexpr float current_draw_scaler_amps = 15.0f;
134 47272379 : for (uint8_t i=0; i<6; i++) {
135 40519182 : const float pwm = input.servos[i];
136 40519182 : const float fraction = fabsf(pwm - 1500) / 500.0f;
137 40519182 : battery_current += fraction * current_draw_scaler_amps;
138 : }
139 6753197 : battery.consume_energy(battery_current, AP_HAL::micros64());
140 6753197 : battery_voltage = battery.get_voltage();
141 6753197 : battery_temperature_degC = battery.get_temperature_degC();
142 6753197 : }
143 :
144 : /**
145 : * @brief Calculate the torque induced by buoyancy foam
146 : *
147 : * @param torque Output torques
148 : */
149 6753197 : void Submarine::calculate_buoyancy_torque(Vector3f &torque)
150 : {
151 : // Let's assume 2 Liters water displacement at the top, and ~ 2kg of weight at the bottom.
152 6753197 : const Vector3f force_up(0,0,-40); // 40 N upwards
153 6753197 : const Vector3f force_position = dcm.transposed() * Vector3f(0, 0, 0.15); // offset in meters
154 6753197 : torque = force_position % force_up;
155 6753197 : }
156 :
157 :
158 : /**
159 : * @brief Calculate sea floor depth from submarine position
160 : * This creates a non planar floor for rangefinder sensor test
161 : * TODO: Create a better sea floor with procedural generatation
162 : *
163 : * @param position
164 : * @return float
165 : */
166 72220215 : float Submarine::calculate_sea_floor_depth(const Vector3d &/*position*/) const
167 : {
168 72220215 : return 50;
169 : }
170 :
171 : /**
172 : * @brief Calculate drag force against body
173 : *
174 : * @param velocity Body frame velocity of fluid
175 : * @param drag_coefficient Drag coefficient of body
176 : * @param force Output forces
177 : * $ F_D = rho * v^2 * A * C_D / 2 $
178 : * rho = water density (kg/m^3), V = velocity (m/s), A = area (m^2), C_D = drag_coefficient
179 : */
180 6753197 : void Submarine::calculate_drag_force(const Vector3f &velocity, const Vector3f &drag_coefficient, Vector3f &force) const
181 : {
182 : /**
183 : * @brief It's necessary to keep the velocity orientation from the body frame.
184 : * To do so, a mathematical artifice is used to do velocity square but without loosing the direction.
185 : * $(|V|/V)*V^2$ = $|V|*V$
186 : */
187 : const Vector3f velocity_2(
188 6753197 : fabsf(velocity.x) * velocity.x,
189 6753197 : fabsf(velocity.y) * velocity.y,
190 6753197 : fabsf(velocity.z) * velocity.z
191 6753197 : );
192 :
193 6753197 : force = (velocity_2 * water_density) * frame_property.equivalent_sphere_area / 2.0f;
194 6753197 : force *= drag_coefficient;
195 6753197 : }
196 :
197 : /**
198 : * @brief Calculate angular drag torque using the equivalente sphere area and assuming a laminar external flow.
199 : *
200 : * $F_D = C_D*A*\rho*V^2/2$
201 : * where:
202 : * $F_D$ is the drag force
203 : * $C_D$ is the drag coefficient
204 : * $A$ is the surface area in contact with the fluid
205 : * $/rho$ is the fluid density (1000kg/m³ for water)
206 : * $V$ is the fluid velocity velocity relative to the surface
207 : *
208 : * @param angular_velocity Body frame velocity of fluid
209 : * @param drag_coefficient Rotational drag coefficient of body
210 : */
211 6753197 : void Submarine::calculate_angular_drag_torque(const Vector3f &angular_velocity, const Vector3f &drag_coefficient, Vector3f &torque) const
212 : {
213 : /**
214 : * @brief It's necessary to keep the velocity orientation from the body frame.
215 : * To do so, a mathematical artifice is used to do velocity square but without loosing the direction.
216 : * $(|V|/V)*V^2$ = $|V|*V$
217 : */
218 : Vector3f v_2(
219 6753197 : fabsf(angular_velocity.x) * angular_velocity.x,
220 6753197 : fabsf(angular_velocity.y) * angular_velocity.y,
221 6753197 : fabsf(angular_velocity.z) * angular_velocity.z
222 6753197 : );
223 6753197 : Vector3f f_d = v_2 *= drag_coefficient * frame_property.equivalent_sphere_area * 1000 / 2;
224 6753197 : torque = f_d * frame_property.equivalent_sphere_radius;
225 6753197 : }
226 :
227 :
228 : /**
229 : * @brief Calculate buoyancy force of the frame
230 : *
231 : * @return float
232 : */
233 6753197 : float Submarine::calculate_buoyancy_acceleration()
234 : {
235 : // position.z is the depth of the centre of the frame, so the frame starts to enter the
236 : // water at -height/2 and is completely below the water level at +height/2
237 6753197 : const float submerged_proportion = constrain_float(position.z + frame_property.height/2, 0, frame_property.height) / frame_property.height;
238 :
239 : // Gravity is applied separately in Aircraft::update_dynamics. This return is buoyant
240 : // acceleration only (upward in NED), which scales with displaced volume. Fully
241 : // submerged with SIM_BUOYANCY=0 yields GRAVITY_MSS, cancelling gravity.
242 6753197 : return submerged_proportion * (GRAVITY_MSS + sitl->buoyancy / frame_property.mass);
243 : };
244 :
245 : /*
246 : update the Submarine simulation by one time step
247 : */
248 6753197 : void Submarine::update(const struct sitl_input &input)
249 : {
250 : // get wind vector setup
251 6753197 : update_wind(input);
252 :
253 6753197 : Vector3f rot_accel;
254 :
255 6753197 : calculate_forces(input, rot_accel, accel_body);
256 6753197 : update_battery(input);
257 :
258 6753197 : update_dynamics(rot_accel);
259 6753197 : update_external_payload(input);
260 :
261 : // update lat/lon/altitude
262 6753197 : update_position();
263 6753197 : time_advance();
264 :
265 : // update magnetic field
266 6753197 : update_mag_field_bf();
267 6753197 : }
268 :
269 : /*
270 : return true if we are on the ground
271 : */
272 20259591 : bool Submarine::on_ground() const
273 : {
274 20259591 : return false;
275 : }
|
