From e31f98157bc7ec3fb358970c386e7fdf634d7087 Mon Sep 17 00:00:00 2001 From: IamPete1 <33176108+IamPete1@users.noreply.github.com> Date: Fri, 5 Oct 2018 23:23:08 +0100 Subject: [PATCH] SITL: sailboat add matlab VPP tool --- libraries/SITL/tools/Sailboat_VPP.m | 159 ++++++++++++++++++++++++++++ 1 file changed, 159 insertions(+) create mode 100644 libraries/SITL/tools/Sailboat_VPP.m diff --git a/libraries/SITL/tools/Sailboat_VPP.m b/libraries/SITL/tools/Sailboat_VPP.m new file mode 100644 index 0000000000..3229f49ccd --- /dev/null +++ b/libraries/SITL/tools/Sailboat_VPP.m @@ -0,0 +1,159 @@ +clc +clear all +close all + +%{ +Very basic Velocity Prediction Program, this duplicates the physics used in the sailboat SITL model. It can be used to generate polar plots to check +STIL physics is sensible. + +In the future it is hoped this program can be used as a comparison for a built in polar generator that is yet to be added to the ArduRover code. +Alternatively a drag polar generated here can be used directly in the code to select optimal sailing angles. + +Sail lift, drag and Hull drag not considered properly, boat heel not considered. The polar generated is however realistic for a 'generic' sailboat. +%} + +sail_aoa_in = [0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150, 160, 170, 179.9]; % deg +lift_curve_in = [0.00, 0.50, 1.00, 1.10, 0.95, 0.75, 0.60, 0.40, 0.20, 0.00, -0.20, -0.40, -0.60, -0.75, -0.95, -1.10, -1.00, -0.50, 0]; +drag_curve_in = [0.10, 0.10, 0.20, 0.40, 0.80, 1.20, 1.50, 1.70, 1.90, 1.95, 1.90, 1.70, 1.50, 1.20, 0.80, 0.40, 0.20, 0.10, 0.10]; + + +min_sail_angle = 0; % deg +max_sail_angle = 90; % deg + +% intpolate to more detailed curves +sail_aoa = sail_aoa_in(1):0.1:sail_aoa_in(end)-0.1; +lift_curve = interp1(sail_aoa_in, lift_curve_in, sail_aoa); +drag_curve = interp1(sail_aoa_in, drag_curve_in, sail_aoa); + +% Plot lift curve +figure +subplot(2,1,1) +hold all +title('Lift and Drag Curves') +scatter(sail_aoa_in,lift_curve_in,'*','r') +plot(sail_aoa,lift_curve,'r') +ylabel('Cl') + +subplot(2,1,2) +hold all +scatter(sail_aoa_in,drag_curve_in,'*','b') +plot(sail_aoa,drag_curve,'b') +xlabel('AoA (deg)') +ylabel('Cd') + +%% VPP +% Do for range of sailing angles and wind speeds, increse speed until foward force equals zero +wind_speed = 3:2:15; % m/s +heading = 0:1:179; % deg +speed_step = 0.01; % m/s + +for n = 1:length(wind_speed) + for m = 1:length(heading) + + thrust = 1; + boat_speed{n,m}(1) = 0; + j = 0; + while thrust > 0 && boat_speed{n,m}(end) < 10 + j = j + 1; + if j == 1 + boat_speed{n,m}(j) = speed_step; + else + boat_speed{n,m}(j) = boat_speed{n,m}(j-1) + speed_step; + end + + % calculate apparent wind in earth-frame (this is the direction the wind is coming from) + wind_apparent_speed{n,m}(j) = sqrt(wind_speed(n)^2 + boat_speed{n,m}(j)^2 + 2 * wind_speed(n) * boat_speed{n,m}(j) * cosd(heading(m))); + + % Aparent wind angle in body frame + wind_apparent_dir_bf{n,m}(j) = acosd((wind_speed(n) * cosd(heading(m)) + boat_speed{n,m}(j))/wind_apparent_speed{n,m}(j)); + + % calculate lift at all angles-of-attack from wind to mainsail + % calculate Lift force (perpendicular to wind direction) and Drag force (parallel to wind direction) + lift = lift_curve * wind_apparent_speed{n,m}(j); + drag = drag_curve * wind_apparent_speed{n,m}(j); + + % rotate lift and drag from wind frame into body frame + sin_rot = sind(wind_apparent_dir_bf{n,m}(j)); + cos_rot = cosd(wind_apparent_dir_bf{n,m}(j)); + force_fwd = (lift * sin_rot) - (drag * cos_rot); + + % accel in body frame due acceleration from sail and deceleration from hull friction + hull_drag = 0.5 * boat_speed{n,m}(j)^2; + thrust_range = force_fwd - hull_drag; + + % angle of sail to boat + sail_boat_angle = wind_apparent_dir_bf{n,m}(j) - sail_aoa; + % discard imposible sail angles + index = find(sail_boat_angle > max_sail_angle | sail_boat_angle < min_sail_angle); + thrust_range(index) = NaN; + + % select sail angle that gives best thrust + index = find(thrust_range == max(thrust_range),1); + + if ~isempty(index) + thrust = thrust_range(index); + sail_force{n,m}(j) = force_fwd(index); + sail_angle_trim{n,m}(j) = sail_boat_angle(index); + hull_drag_trim{n,m}(j) = hull_drag; + else + thrust = 0; + end + end + + max_boat_speed(n,m) = boat_speed{n,m}(j); + % Set minumum values to NaN to make plot neater + if max_boat_speed(n,m) == speed_step + max_boat_speed(n,m) = NaN; + end + end +end + + +%% Polar plot +figure +for n = 1:length(wind_speed) + polarplot(deg2rad(heading),max_boat_speed(n,:)) + hold on + legend_val{n} = sprintf('%g m/s',wind_speed(n)); +end +ax = gca; +d = ax.ThetaDir; +ax.ThetaDir = 'clockwise'; +ax.ThetaZeroLocation = 'top'; +legend(legend_val) +title('Polar Plot, boat speed (m/s) vs heading (deg)') + + +%% Plot sail thrust and hull drag for a given heading +plot_heading = 50; +plot_wind_speed = 5; + +% convet from values to indexs +plot_heading = find(heading == plot_heading); +plot_wind_speed = find(wind_speed == plot_wind_speed); + +if isempty(plot_heading) || isempty(plot_wind_speed) + fprintf('Plot headings and or windspeed not within the range caulated') + return +end + +figure +subplot(2,1,1) +hold all +title(sprintf('Thrust vs drag for heading of %g deg and wind speed of %g m/s',heading(plot_heading),wind_speed(plot_wind_speed))) +yyaxis left +plot(boat_speed{plot_wind_speed,plot_heading},wind_apparent_speed{plot_wind_speed,plot_heading}) +ylabel('Apparent wind speed (m/s)') +yyaxis right +plot(boat_speed{plot_wind_speed,plot_heading},wind_apparent_dir_bf{plot_wind_speed,plot_heading}) +ylabel('Apparent wind direction (deg)') + + +subplot(2,1,2) +hold all +plot(boat_speed{plot_wind_speed,plot_heading},sail_force{plot_wind_speed,plot_heading}) +plot(boat_speed{plot_wind_speed,plot_heading},hull_drag_trim{plot_wind_speed,plot_heading}) + +legend('Sail Forwards force','Hull drag','Location','northwest') +ylabel('Force (N)') +xlabel('Boat speed (m/s)')