From a8273753087315a4766fd6ad078396dcf71ad69a Mon Sep 17 00:00:00 2001 From: Mario Bras Date: Mon, 6 Jul 2020 12:32:46 -0700 Subject: [PATCH] Rework MTOW calculations and automate design point calculation process - Fix time calculation bugs in MTOW - Rework the way battery and fuel weight is calculated in MTOW function - Add transition mission segment type - Fix MATLAB's JSON output with `prettyjson` function - Update mission plotting - Update drag and lift slope calculation - Automate design point calculation - Update example project file --- adt.m | 31 ++-- build_mission.m | 49 ++++++ design_point.m | 350 +++++++++++++++++++++++++++++------------- drag_buildup.m | 45 ++++-- estimate_mf_prop.m | 19 --- estimate_mf_subs.m | 9 -- examples/project.json | 249 ++++++++++++++++++------------ friction_coeff.m | 2 +- fuselage.m | 18 --- lift_slope.m | 29 ++++ lifting_surface.m | 24 --- load_project.m | 12 +- mtow.m | 271 +++++++++++++++++--------------- plot_mission.m | 49 ++++-- prettyjson.m | 71 +++++++++ save_project.m | 5 +- 16 files changed, 781 insertions(+), 452 deletions(-) create mode 100644 build_mission.m delete mode 100644 estimate_mf_prop.m delete mode 100644 estimate_mf_subs.m delete mode 100644 fuselage.m create mode 100644 lift_slope.m delete mode 100644 lifting_surface.m create mode 100644 prettyjson.m diff --git a/adt.m b/adt.m index b4bcb01..533f6bc 100644 --- a/adt.m +++ b/adt.m @@ -4,11 +4,9 @@ % % This file is subject to the license terms in the LICENSE file included in this distribution -function data = adt(open_filename, varargin) -close all; - +function data = adt(filename, varargin) %% Load project file -data = load_project(open_filename); +data = load_project(filename); %% Concept data.concept = ahp(data.concept); @@ -18,30 +16,27 @@ constants.g = 9.81; % m/s^2 constants.e = 0.85; % Oswald efficiency number. -data.aircraft.propulsion.total_efficiency = data.aircraft.propulsion.prop_efficiency * data.aircraft.propulsion.gear_efficiency * data.aircraft.propulsion.em_efficiency * data.aircraft.propulsion.esc_efficiency * data.aircraft.propulsion.dist_efficiency; -data.mission.mf_prop = estimate_mf_prop(1.2, 0.04, 1.5, 0.3, 0.1, 0.5, data.aircraft.propulsion.config); -data.mission.mf_subs = estimate_mf_subs(data.mission.mf_prop, data.mission.mf_struct); - -% Take-off mass estimation -[data.mission, data.aircraft] = mtow(data.mission, data.aircraft, constants); +% Build mission profile +data.mission = build_mission(data.mission); %% Plot mission profile plot_mission(data.mission); +% Take-off mass estimation +[data.mission, data.aircraft] = mtow(data.mission, data.aircraft, constants); + %% Design point calculation data.aircraft = design_point(data.mission, data.aircraft, constants); -%% Engine selection -data.mission.mf_prop = data.aircraft.propulsion.mass / data.aircraft.mass_to; - -%% Lifting surface design -data.aircraft = lifting_surface(data.aircraft, data.mission.segments{3}); +% Recalculate take-off mass based on updated values from design point +% [data.mission, data.aircraft] = mtow(data.mission, data.aircraft, constants); +% data.aircraft = design_point(data.mission, data.aircraft, constants); -%% Fuselage design -data.aircraft = fuselage(data.aircraft, data.mission.segments{3}); +%% Lift curve slope +data.aircraft = lift_slope(data.aircraft, data.mission); %% Drag buildup -data.aircraft = drag_buildup(data.aircraft); +data.aircraft = drag_buildup(data.aircraft, data.mission); %% Save new project file if ~isempty(varargin) diff --git a/build_mission.m b/build_mission.m new file mode 100644 index 0000000..09c442a --- /dev/null +++ b/build_mission.m @@ -0,0 +1,49 @@ +% Aircraft design tool +% +% Mario Bras (mbras@uvic.ca) and Ricardo Marques (ricardoemarques@uvic.ca) 2019 +% +% This file is subject to the license terms in the LICENSE file included in this distribution + +function mission = build_mission(mission) + +% Iterate over mission segments +for i = 1 : length(mission.segments) + [mission.segments{i}.temperature, mission.segments{i}.speed_sound, mission.segments{i}.pressure, mission.segments{i}.density] = atmosisa(mission.segments{i}.altitude); + + if ~isfield(mission.segments{i}, 'velocity') + mission.segments{i}.velocity = 0; + end + + if strcmp(mission.segments{i}.type, 'taxi') % Taxi segment + mission.segments{i}.range = 0; + elseif strcmp(mission.segments{i}.type, 'hover') % Hover segment + mission.segments{i}.range = 0; + elseif strcmp(mission.segments{i}.type, 'transition') % Transition segment + mission.segments{i}.range = mean(mission.segments{i}.velocity) * mission.segments{i}.time; + elseif strcmp(mission.segments{i}.type, 'climb') % Climb segment + mission.segments{i}.time = (mission.segments{i}.altitude(2) - mission.segments{i}.altitude(1)) / mission.segments{i}.velocity / sind(mission.segments{i}.angle); + mission.segments{i}.range = mission.segments{i}.velocity * mission.segments{i}.time * cosd(mission.segments{i}.angle); + elseif strcmp(mission.segments{i}.type, 'vertical_climb') % Vertical climb segment + mission.segments{i}.time = (mission.segments{i}.altitude(2) - mission.segments{i}.altitude(1)) / mission.segments{i}.velocity; + mission.segments{i}.range = 0; + elseif strcmp(mission.segments{i}.type, 'acceleration') % Acceleration segment + mission.segments{i}.range = mission.segments{i}.velocity * mission.segments{i}.time; + elseif strcmp(mission.segments{i}.type, 'cruise') % Cruise segment + mission.segments{i}.time = mission.segments{i}.range / mission.segments{i}.velocity; + elseif strcmp(mission.segments{i}.type, 'hold') % Hold segment + mission.segments{i}.range = mission.segments{i}.velocity * mission.segments{i}.time; + % elseif strcmp(mission.segments{i}.type, 'combat') % Combat segment + elseif strcmp(mission.segments{i}.type, 'descent') % Descent segment + mission.segments{i}.time = abs(mission.segments{i}.altitude(2) - mission.segments{i}.altitude(1)) / abs(mission.segments{i}.velocity) / abs(sind(mission.segments{i}.angle)); + mission.segments{i}.range = abs(mission.segments{i}.velocity) * mission.segments{i}.time * abs(cosd(mission.segments{i}.angle)); + elseif strcmp(mission.segments{i}.type, 'vertical_descent') % Vertical descent segment + mission.segments{i}.time = abs(mission.segments{i}.altitude(2) - mission.segments{i}.altitude(1)) / abs(mission.segments{i}.velocity); + mission.segments{i}.range = 0; + elseif strcmp(mission.segments{i}.type, 'landing') % Landing segment + mission.segments{i}.range = 0; + elseif strcmp(mission.segments{i}.type, 'drop') % Drop segment + mission.segments{i}.range = 0; + elseif strcmp(mission.segments{i}.type, 'load') % Load segment + mission.segments{i}.range = 0; + end +end diff --git a/design_point.m b/design_point.m index 5e3e2ca..6056436 100644 --- a/design_point.m +++ b/design_point.m @@ -6,72 +6,139 @@ function aircraft = design_point(mission, aircraft, constants) -wl = 0:10:10000; -ptl = 0:0.001:0.1; -dl = 0:10:10000; - -f(1) = figure; -legend; -hold on; -a(1) = gca; - -f(2) = figure; +wl = 0:5:1000; +dl = 0:5:10000; +pl = 0:0.0005:0.3; +[plf_grid, wl_grid] = meshgrid(pl, wl); +[plv_grid, dl_grid] = meshgrid(pl, dl); +cf = ones(length(wl), length(pl)); +cv = ones(length(dl), length(pl)); + +% Configure plot +colors = {'#0072BD','#D95319','#EDB120','#7E2F8E','#77AC30','#4DBEEE','#A2142F'}; +fig = figure('DeleteFcn', 'doc datacursormode'); +yyaxis right; legend; hold on; -a(2) = gca; - -a(1).XLim = [0 wl(end)]; -a(1).YLim = [0 ptl(end)]; -a(1).XLabel.String = 'W/S'; -a(1).YLabel.String = 'W/P'; -a(1).Title.String = 'Forward Flight Design Point'; - -a(2).XLim = [0 dl(end)]; -a(2).YLim = [0 ptl(end)]; -a(2).XLabel.String = 'W/A'; -a(2).YLabel.String = 'W/P'; -a(2).Title.String = 'Vertical Flight Design Point'; +a = gca; +a.Title.String = 'Design Point'; +a.XLim = [0 pl(end)]; +a.XLabel.String = 'W/P'; +a.YLim = [0 dl(end)]; +a.YLabel.String = 'W/A'; +a.LineStyleOrder = '-'; +colororder(colors) +yyaxis left; +a.YLim = [0 wl(end)]; +a.YLabel.String = 'W/S'; +a.LineStyleOrder = '-'; +colororder(colors) k = 1 / pi / aircraft.performance.aspect_ratio / constants.e; -%% Forward flight constraints -% EDIT NUMBER/TYPE OF CONSTRAINTS TO PLOT -r = range(mission.segments{3}.density, mission.segments{3}.velocity, aircraft.performance.c_d0, k, mission.segments{3}.propulsion); -e = endurance(mission.segments{6}.density, mission.segments{6}.velocity, aircraft.performance.c_d0, k, mission.segments{6}.propulsion); -cr = cruise_speed(wl, mission.segments{3}.density, mission.segments{3}.velocity, aircraft.performance.c_d0, k, aircraft.propulsion.prop_efficiency, mission.segments{3}.propulsion); -cl1 = climb_angle(wl, mission.segments{8}.density(1), aircraft.performance.c_d0, k, mission.segments{8}.angle, aircraft.propulsion.prop_efficiency, mission.segments{8}.propulsion); -cl2 = climb_angle(wl, mission.segments{8}.density(2), aircraft.performance.c_d0, k, mission.segments{8}.angle, aircraft.propulsion.prop_efficiency, mission.segments{8}.propulsion); - -figure(f(1)); -plot(a(1), [r r], [ptl(1) ptl(end)], 'DisplayName', 'Range'); -plot(a(1), [e e], [ptl(1) ptl(end)], 'DisplayName', 'Endurance'); -plot(a(1), wl, cr, 'DisplayName', 'Cruise Speed'); -plot(a(1), wl, cl1, 'DisplayName', strcat('Climb Angle @ \rho = ', num2str(mission.segments{8}.density(1)))); -plot(a(1), wl, cl2, 'DisplayName', strcat('Climb Angle @ \rho = ', num2str(mission.segments{8}.density(2)))); +% Iterate over horizontal flight mission segments +for i = 1 : length(mission.segments) +% if strcmp(mission.segments{i}.type, 'hover') % Hover segment +% pti = find_propulsion_type_index(aircraft, mission.segments{i}.propulsion_type); +% cv = cv .* hover_region(plv_grid, dl_grid, mission.segments{i}.density, aircraft.propulsion_types{pti}.fm); +% h = hover(dl, mission.segments{i}.density, aircraft.propulsion_types{pti}.fm); +% cv = cv .* transition_region(plv_grid, wl_grid, dl_grid, mission.segments{i}.density, k, aircraft.performance.cd_0, aircraft.propulsion_types{pti}.tip_velocity, aircraft.propulsion_types{pti}.ss, aircraft.propulsion_types{pti}.cd, aircraft.propulsion_types{pti}.k_i, mission.segments{i+1}.velocity, mission.segments{i}.transition_angle); +% tr = transition(wl, dl, mission.segments{i}.density, k, aircraft.performance.cd_0, aircraft.propulsion_types{pti}.tip_velocity, aircraft.propulsion_types{pti}.ss, aircraft.propulsion_types{pti}.cd, aircraft.propulsion_types{pti}.k_i, mission.segments{i+1}.velocity, mission.segments{i}.transition_angle); +% yyaxis right; +% plot(h, dl, 'DisplayName', 'Hover', strcat("Hover constraint, segment ", int2str(i), ' (hover)')); +% plot(tr, dl, 'DisplayName', strcat("Transition constraint, \theta = ", num2str(mission.segments{i}.transition_angle), ", segment ", int2str(i), ' (hover)')); + if strcmp(mission.segments{i}.type, 'climb') % Climb segment + pti = find_propulsion_type_index(aircraft, mission.segments{i}.propulsion_type); + cf = cf .* climb_region(plf_grid, wl_grid, mission.segments{i}.density(1), mission.segments{i}.velocity, aircraft.performance.cd_0, k, mission.segments{i}.angle, aircraft.propulsion_types{pti}); + cl = climb(wl, mission.segments{i}.density(1), mission.segments{i}.velocity, aircraft.performance.cd_0, k, mission.segments{i}.angle, aircraft.propulsion_types{pti}); + yyaxis left; + plot(cl, wl, 'DisplayName', strcat("Climb constraint, segment ", int2str(i), ' (climb)')); +% elseif strcmp(mission.segments{i}.type, 'vertical_climb') % Vertical climb segment +% pti = find_propulsion_type_index(aircraft, mission.segments{i}.propulsion_type); +% cv = cv .* vertical_climb_region(plv_grid, dl_grid, mission.segments{i}.density(1), aircraft.propulsion_types{pti}.tip_velocity, aircraft.propulsion_types{pti}.ss, aircraft.propulsion_types{pti}.cd, aircraft.propulsion_types{pti}.k_i, mission.segments{i}.velocity); +% vcl = vertical_climb(dl, mission.segments{i}.density(1), aircraft.propulsion_types{pti}.tip_velocity, aircraft.propulsion_types{pti}.ss, aircraft.propulsion_types{pti}.cd, aircraft.propulsion_types{pti}.k_i, mission.segments{i}.velocity); +% yyaxis right; +% plot(vcl, dl, 'DisplayName', 'Vertical Climb', 'DisplayName', strcat("Vertical climb constraint, segment ", int2str(i), ' (vertical climb)')); + elseif strcmp(mission.segments{i}.type, 'cruise') % Cruise segment + pti = find_propulsion_type_index(aircraft, mission.segments{i}.propulsion_type); + cf = cf .* range_region(wl_grid, mission.segments{i}.density, mission.segments{i}.velocity, aircraft.performance.cd_0, k, aircraft.propulsion_types{pti}); + r = range(mission.segments{i}.density, mission.segments{i}.velocity, aircraft.performance.cd_0, k, aircraft.propulsion_types{pti}); + cf = cf .* cruise_speed_region(plf_grid, wl_grid, mission.segments{i}.density, mission.segments{i}.velocity, aircraft.performance.cd_0, k, aircraft.propulsion_types{pti}); + cr = cruise_speed(wl, mission.segments{i}.density, mission.segments{i}.velocity, aircraft.performance.cd_0, k, aircraft.propulsion_types{pti}); + yyaxis left; + plot([pl(1) pl(end)], [r r], 'DisplayName', strcat("Range constraint, segment ", int2str(i), ' (cruise)')); + plot(cr, wl, 'DisplayName', strcat("Cruise speed constraint, segment ", int2str(i), ' (cruise)')); + elseif strcmp(mission.segments{i}.type, 'hold') % Hold segment + pti = find_propulsion_type_index(aircraft, mission.segments{i}.propulsion_type); + cf = cf .* endurance_region(wl_grid, mission.segments{i}.density, mission.segments{i}.velocity, aircraft.performance.cd_0, k, aircraft.propulsion_types{pti}); + e = endurance(mission.segments{i}.density, mission.segments{i}.velocity, aircraft.performance.cd_0, k, aircraft.propulsion_types{pti}); + yyaxis left; + plot([pl(1) pl(end)], [e e], 'DisplayName', strcat("Endurance constraint, segment ", int2str(i), ' (hold)')); + elseif strcmp(mission.segments{i}.type, 'descent') % Descent segment +% elseif strcmp(mission.segments{i}.type, 'vertical_descent') % Vertical descent segment +% cv = cv .* vertical_descent_region(plv_grid, dl_grid, mission.segments{i}.density(1), aircraft.propulsion_types{pti}.tip_velocity, aircraft.propulsion_types{pti}.ss, aircraft.propulsion_types{pti}.cd, aircraft.propulsion_types{pti}.k_i, mission.segments{i}.velocity); +% vd = vertical_descent(dl, mission.segments{i}.density(1), aircraft.propulsion_types{pti}.tip_velocity, aircraft.propulsion_types{pti}.ss, aircraft.propulsion_types{pti}.cd, aircraft.propulsion_types{pti}.k_i, mission.segments{i}.velocity); +% yyaxis right; +% plot(vd, dl, ''DisplayName', strcat("Vertical descent constraint, segment ", int2str(i), ' (vertical descent)')); + end +end + +cf(~cf) = NaN; +yyaxis left; +sf = surf(pl, wl, cf, 'FaceAlpha', 0.2, 'FaceColor', '#0072BD', 'EdgeColor', 'none', 'DisplayName', 'Forward Flight Design Space'); % Pick forward flight wing and power loading -[aircraft.performance.wl, aircraft.performance.pl_fwd] = ginput(1); -aircraft.performance.area_ref = aircraft.weight_to / aircraft.performance.wl; -aircraft.performance.b_ref = sqrt(aircraft.performance.aspect_ratio * aircraft.performance.area_ref); -aircraft.performance.c_ref = aircraft.performance.area_ref / aircraft.performance.b_ref; -aircraft.propulsion.fwd_power = aircraft.weight_to / aircraft.performance.pl_fwd; - -%% Vertical flight constraints -% EDIT NUMBER/TYPE OF CONSTRAINTS TO PLOT -h = hover(dl, mission.segments{3}.density, aircraft.propulsion.fm); -vcl = vertical_climb(dl, mission.segments{2}.density(1), aircraft.propulsion.tip_velocity, aircraft.propulsion.ss, aircraft.propulsion.c_d, aircraft.propulsion.k_i, mission.segments{2}.velocity); -tr = transition(aircraft.performance.wl, dl, mission.segments{3}.density, k, aircraft.performance.c_d0, aircraft.propulsion.tip_velocity, aircraft.propulsion.ss, aircraft.propulsion.c_d, aircraft.propulsion.k_i, mission.segments{3}.velocity, aircraft.propulsion.tt_tilt); - -figure(f(2)); -plot(a(2), [dl(1) dl(end)], [aircraft.performance.pl_fwd aircraft.performance.pl_fwd], 'k--', 'DisplayName', 'Forward Flight W/P'); -plot(a(2), dl, h, 'DisplayName', 'Hover'); -plot(a(2), dl, vcl, 'DisplayName', 'Vertical Climb'); -plot(a(2), dl, tr, 'DisplayName', strcat('Transition @ \theta = ', num2str(aircraft.propulsion.tt_tilt))); +yyaxis left; +[x, y] = ginput(1); +scatter(x, y, 'filled', 'MarkerEdgeColor', '#0072BD', 'MarkerFaceColor', '#0072BD', 'DisplayName', 'Forward Flight Design Point'); + +aircraft.performance.power_loading_fwd = x; +aircraft.performance.wing_loading = y; +aircraft.performance.wing_area_ref = aircraft.mass * constants.g / aircraft.performance.wing_loading; +aircraft.performance.wing_span_ref = sqrt(aircraft.performance.aspect_ratio * aircraft.performance.wing_area_ref); +aircraft.performance.wing_chord_ref = aircraft.performance.wing_area_ref / aircraft.performance.wing_span_ref; +aircraft.performance.power_fwd = aircraft.mass * constants.g / aircraft.performance.power_loading_fwd; + +% Iterate over vertical flight mission segments +for i = 1 : length(mission.segments) + if strcmp(mission.segments{i}.type, 'hover') % Hover segment + pti = find_propulsion_type_index(aircraft, mission.segments{i}.propulsion_type); + cv = cv .* hover_region(plv_grid, dl_grid, mission.segments{i}.density, aircraft.propulsion_types{pti}.fm); + h = hover(dl, mission.segments{i}.density, aircraft.propulsion_types{pti}.fm); + yyaxis right; + plot(h, dl, 'DisplayName', strcat("Hover constraint, segment ", int2str(i), ' (hover)')); + elseif strcmp(mission.segments{i}.type, 'transition') % Transition segment + cv = cv .* transition_region(plv_grid, aircraft.performance.wing_loading, dl_grid, mission.segments{i}.density, k, aircraft.performance.cd_0, aircraft.propulsion_types{pti}.tip_velocity, aircraft.propulsion_types{pti}.ss, aircraft.propulsion_types{pti}.cd, aircraft.propulsion_types{pti}.k_i, mission.segments{i+1}.velocity, mission.segments{i}.transition_angle); + tr = transition(aircraft.performance.wing_loading, dl, mission.segments{i}.density, k, aircraft.performance.cd_0, aircraft.propulsion_types{pti}.tip_velocity, aircraft.propulsion_types{pti}.ss, aircraft.propulsion_types{pti}.cd, aircraft.propulsion_types{pti}.k_i, mission.segments{i+1}.velocity, mission.segments{i}.transition_angle); + yyaxis right; + plot(tr, dl, 'DisplayName', strcat("Transition constraint, \theta = ", num2str(mission.segments{i}.transition_angle), ", segment ", int2str(i), ' (hover)')); + elseif strcmp(mission.segments{i}.type, 'vertical_climb') % Vertical climb segment + pti = find_propulsion_type_index(aircraft, mission.segments{i}.propulsion_type); + cv = cv .* vertical_climb_region(plv_grid, dl_grid, mission.segments{i}.density(1), aircraft.propulsion_types{pti}.tip_velocity, aircraft.propulsion_types{pti}.ss, aircraft.propulsion_types{pti}.cd, aircraft.propulsion_types{pti}.k_i, mission.segments{i}.velocity); + vcl = vertical_climb(dl, mission.segments{i}.density(1), aircraft.propulsion_types{pti}.tip_velocity, aircraft.propulsion_types{pti}.ss, aircraft.propulsion_types{pti}.cd, aircraft.propulsion_types{pti}.k_i, mission.segments{i}.velocity); + yyaxis right; + plot(vcl, dl, 'DisplayName', strcat("Vertical climb constraint, segment ", int2str(i), ' (vertical climb)')); + elseif strcmp(mission.segments{i}.type, 'vertical_descent') % Vertical descent segment +% cv = cv .* vertical_descent_region(plv_grid, dl_grid, mission.segments{i}.density(1), aircraft.propulsion_types{pti}.tip_velocity, aircraft.propulsion_types{pti}.ss, aircraft.propulsion_types{pti}.cd, aircraft.propulsion_types{pti}.k_i, mission.segments{i}.velocity); +% vd = vertical_descent(dl, mission.segments{i}.density(1), aircraft.propulsion_types{pti}.tip_velocity, aircraft.propulsion_types{pti}.ss, aircraft.propulsion_types{pti}.cd, aircraft.propulsion_types{pti}.k_i, mission.segments{i}.velocity); +% yyaxis right; +% plot(vd, dl, 'DisplayName', strcat("Vertical descent constraint, segment ", int2str(i), ' (vertical descent)')); + end +end + +cv(~cv) = NaN; +yyaxis right; +sv = surf(pl, dl, cv, 'FaceAlpha', 0.2, 'FaceColor', '#D95319', 'EdgeColor', 'none', 'DisplayName', 'Vertical Flight Design Space'); % Pick vertical flight disk and power loading -[aircraft.performance.dl, aircraft.performance.pl_vert] = ginput(1); -aircraft.propulsion.a = aircraft.weight_to / aircraft.performance.dl; -aircraft.propulsion.vert_power = aircraft.weight_to / aircraft.performance.pl_vert; +yyaxis right; +[x, y] = ginput(1); +scatter(x, y, 'filled', 'MarkerEdgeColor', '#D95319', 'MarkerFaceColor', '#D95319', 'DisplayName', 'Vertical Flight Design Point'); + +aircraft.performance.power_loading_vert = x; +aircraft.performance.disc_loading = y; +aircraft.performance.disc_area_ref = aircraft.mass * constants.g / aircraft.performance.disc_loading; +aircraft.performance.power_vert = aircraft.mass * constants.g / aircraft.performance.power_loading_vert; %% Performance functions function test = is_jet(propulsion) @@ -80,101 +147,174 @@ function test = is_prop(propulsion) test = strcmp(propulsion.type, 'prop'); -function v = v_min_thrust(wl, rho, k, c_d0) -v = sqrt(2 * wl / rho * sqrt(k / c_d0)); +function v = v_min_thrust(wl, rho, k, cd_0) +v = sqrt(2 * wl / rho * sqrt(k / cd_0)); -function v = v_min_power(wl, rho, k, c_d0) -v = sqrt(2 * wl / rho * sqrt(k / 3 / c_d0)); +function v = v_min_power(wl, rho, k, cd_0) +v = sqrt(2 * wl / rho * sqrt(k / 3 / cd_0)); -function v = v_best_climb_rate_jet(tl, wl, rho, k, c_d0) -v = sqrt(wl / 3 / rho / c_d0 * (1 / tl + sqrt(1 / tl^2 + 12 * c_d0 * k))); +function v = v_best_climb_rate_jet(tl, wl, rho, k, cd_0) +v = sqrt(wl / 3 / rho / cd_0 * (1 / tl + sqrt(1 / tl^2 + 12 * cd_0 * k))); -function v = v_best_climb_rate_prop(wl, rho, k, c_d0) -v = v_min_power(wl, rho, k, c_d0); +function v = v_best_climb_rate_prop(wl, rho, k, cd_0) +v = v_min_power(wl, rho, k, cd_0); -function v = v_best_climb_angle_jet(wl, rho, k, c_d0) -v = v_min_thrust(wl, rho, k, c_d0); +function v = v_best_climb_angle_jet(wl, rho, k, cd_0) +v = v_min_thrust(wl, rho, k, cd_0); -function v = v_best_climb_angle_prop(wl, rho, k, c_d0) -v = 0.875 * v_best_climb_rate_prop(wl, rho, k, c_d0); % Raymer pp. 466 +function v = v_best_climb_angle_prop(wl, rho, k, cd_0) +v = 0.875 * v_best_climb_rate_prop(wl, rho, k, cd_0); % Raymer pp. 466 -function c_l = cl_min_thrust(k, c_d0) -c_l = sqrt(c_d0 / k); +function c_l = cl_min_thrust(k, cd_0) +c_l = sqrt(cd_0 / k); -function c_l = cl_min_power(k, c_d0) -c_l = sqrt(3 * c_d0 / k); +function c_l = cl_min_power(k, cd_0) +c_l = sqrt(3 * cd_0 / k); -function c_d = cd_min_thrust(c_d0) -c_d = 2 * c_d0; +function cd = cd_min_thrust(cd_0) +cd = 2 * cd_0; -function c_d = cd_min_power(c_d0) -c_d = 4 * c_d0; +function cd = cd_min_power(cd_0) +cd = 4 * cd_0; %% Vertical flight constraint functions function pl = hover(dl, rho, fm) pl = fm .* sqrt(2 .* rho ./ dl); -function pl = vertical_climb(dl, rho, v_tip, ss, c_d, k_i, v_y) -pl = 1 ./ (v_y - k_i .* v_y ./ 2 + k_i .* sqrt(v_y.^2 + 2 .* dl ./ rho) ./ 2 + rho .* v_tip.^3 .* ss .* c_d ./ dl ./ 8); +function pl = vertical_climb(dl, rho, v_tip, ss, cd, k_i, v_y) +pl = 1 ./ (v_y - k_i .* v_y ./ 2 + k_i .* sqrt(v_y.^2 + 2 .* dl ./ rho) ./ 2 + rho .* v_tip.^3 .* ss .* cd ./ dl ./ 8); -function pl = vertical_descent(dl, rho, v_tip, ss, c_d, k_i, v_y, v_i) +function pl = vertical_descent(dl, rho, v_tip, ss, cd, k_i, v_y, v_i) if v_y / v_i <= -2 % If this condition is met, the vertical climb equation is used for descent, else, an empirical equation is employed - pl = 1 ./ (v_y - k_i ./ 2 * (v_y + sqrt(v_y.^2 - 2 .* dl ./ rho)) + rho .* v_tip.^3 .* ss .* c_d ./ dl ./ 8); + pl = 1 ./ (v_y - k_i ./ 2 * (v_y + sqrt(v_y.^2 - 2 .* dl ./ rho)) + rho .* v_tip.^3 .* ss .* cd ./ dl ./ 8); else v_d = v_i * (k_i - 1.125 * v_y / v_i - 1.372 * (v_y / v_i)^2 - 1.718 * (v_y / v_i)^3 - 0.655 * (v_y / v_i)^4); % Induced velocity in descent according to an empirical relation (see lecture slides) - pl = 1 ./ (v_y + k_i .* v_d + rho .* v_tip.^3 ./ dl .* ss .* c_d ./ 8); + pl = 1 ./ (v_y + k_i .* v_d + rho .* v_tip.^3 ./ dl .* ss .* cd ./ 8); end -function pl = transition(wl, dl, rho, k, c_d0, v_tip, ss, c_d, k_i, v, tt_tilt) +function pl = transition(wl, dl, rho, k, cd_0, v_tip, ss, cd, k_i, v, tt_tilt) aa = 0; % Assuming zero angle of attack of the blades mm = v * cosd(aa) / v_tip; -pl = 1 ./ (k_i ./ sind(tt_tilt) .* sqrt(-v.^2 ./ 2 + sqrt((v.^2 ./ 2).^2 + (dl ./ 2 ./ rho ./ sind(tt_tilt)).^2)) + rho .* v_tip.^3 ./ dl .* (ss .* c_d ./ 8 .* (1 + 4.6 .* mm.^2)) + 0.5 .* rho .* v^3 .* c_d0 ./ wl + 2 .* wl .* k ./ rho ./ v); +pl = 1 ./ (k_i ./ sind(tt_tilt) .* sqrt(-v.^2 ./ 2 + sqrt((v.^2 ./ 2).^2 + (dl ./ 2 ./ rho ./ sind(tt_tilt)).^2)) + rho .* v_tip.^3 ./ dl .* (ss .* cd ./ 8 .* (1 + 4.6 .* mm.^2)) + 0.5 .* rho .* v^3 .* cd_0 ./ wl + 2 .* wl .* k ./ rho ./ v); + +%% Vertical flight constraint regions +function c = hover_region(pl, dl, rho, fm) +c = pl < fm .* sqrt(2 .* rho ./ dl); -%% Forward flight constraint functions -function wl = range(rho, v, c_d0, k, propulsion) +function c = vertical_climb_region(pl, dl, rho, v_tip, ss, cd, k_i, v_y) +c = pl < 1 ./ (v_y - k_i .* v_y ./ 2 + k_i .* sqrt(v_y.^2 + 2 .* dl ./ rho) ./ 2 + rho .* v_tip.^3 .* ss .* cd ./ dl ./ 8); + +function c = vertical_descent_region(pl, dl, rho, v_tip, ss, cd, k_i, v_y, v_i) +if v_y / v_i <= -2 % If this condition is met, the vertical climb equation is used for descent, else, an empirical equation is employed + c = pl < 1 ./ (v_y - k_i ./ 2 * (v_y + sqrt(v_y.^2 - 2 .* dl ./ rho)) + rho .* v_tip.^3 .* ss .* cd ./ dl ./ 8); +else + v_d = v_i * (k_i - 1.125 * v_y / v_i - 1.372 * (v_y / v_i)^2 - 1.718 * (v_y / v_i)^3 - 0.655 * (v_y / v_i)^4); % Induced velocity in descent according to an empirical relation (see lecture slides) + c = pl < 1 ./ (v_y + k_i .* v_d + rho .* v_tip.^3 ./ dl .* ss .* cd ./ 8); +end + +function c = transition_region(pl, wl, dl, rho, k, cd_0, v_tip, ss, cd, k_i, v, tt_tilt) +aa = 0; % Assuming zero angle of attack of the blades +mm = v * cosd(aa) / v_tip; +c = pl < 1 ./ (k_i ./ sind(tt_tilt) .* sqrt(-v.^2 ./ 2 + sqrt((v.^2 ./ 2).^2 + (dl ./ 2 ./ rho ./ sind(tt_tilt)).^2)) + rho .* v_tip.^3 ./ dl .* (ss .* cd ./ 8 .* (1 + 4.6 .* mm.^2)) + 0.5 .* rho .* v^3 .* cd_0 ./ wl + 2 .* wl .* k ./ rho ./ v); + +% Forward flight constraint functions +function wl = range(rho, v, cd_0, k, propulsion) if is_jet(propulsion) - wl = 0.5 * rho * v^2 * sqrt(c_d0 / 3 / k); + wl = 0.5 * rho * v^2 * sqrt(cd_0 / 3 / k); elseif is_prop(propulsion) - wl = 0.5 * rho * v^2 * sqrt(c_d0 / k); + wl = 0.5 * rho * v^2 * sqrt(cd_0 / k); end -function wl = endurance(rho, v, c_d0, k, propulsion) +function wl = endurance(rho, v, cd_0, k, propulsion) if is_jet(propulsion) - wl = 0.5 * rho * v^2 * sqrt(c_d0 / k); + wl = 0.5 * rho * v^2 * sqrt(cd_0 / k); elseif is_prop(propulsion) - wl = 0.5 * rho * v^2 * sqrt(3 * c_d0 / k); + wl = 0.5 * rho * v^2 * sqrt(3 * cd_0 / k); end function wl = stall_speed(rho, v_s, c_lmax) wl = 0.5 * rho * v_s^2 * c_lmax; -function ptl = cruise_speed(wl, rho, v, c_d0, k, ee_prop, propulsion) +function ptl = cruise_speed(wl, rho, v, cd_0, k, propulsion) if is_jet(propulsion) - ptl = 1 ./ (rho .* v.^2 .* c_d0 ./ 2 ./ wl + 2 .* k .* wl ./ rho ./ v.^2); + ptl = 1 ./ (rho .* v.^2 .* cd_0 ./ 2 ./ wl + 2 .* k .* wl ./ rho ./ v.^2); elseif is_prop(propulsion) - ptl = ee_prop ./ (rho .* v.^3 .* c_d0 ./ 2 ./ wl + 2 .* k .* wl ./ rho ./ v); + ptl = propulsion.efficiency ./ (rho .* v.^3 .* cd_0 ./ 2 ./ wl + 2 .* k .* wl ./ rho ./ v); end -function ptl = climb(wl, rho, v, c_d0, k, gg, ee_prop, propulsion) +function pl = climb(wl, rho, v, cd_0, k, gg, propulsion) if is_jet(propulsion) - ptl = 1 ./ (sind(gg) + rho .* v.^2 .* c_d0 ./ 2 ./ wl + 2 .* k .* wl ./ rho ./ v.^2); + pl = 1 ./ (sind(gg) + rho .* v.^2 .* cd_0 ./ 2 ./ wl + 2 .* k .* wl ./ rho ./ v.^2); elseif is_prop(propulsion) - ptl = ee_prop ./ (v .* sind(gg) + rho .* v.^3 .* c_d0 ./ 2 ./ wl + 2 .* k .* wl ./ rho ./ v); + pl = propulsion.efficiency ./ (v .* sind(gg) + rho .* v.^3 .* cd_0 ./ 2 ./ wl + 2 .* k .* wl ./ rho ./ v); end -function ptl = climb_angle(wl, rho, c_d0, k, gg, ee_prop, propulsion) +function pl = climb_angle(wl, rho, cd_0, k, gg, propulsion) if is_jet(propulsion) - ptl = climb(wl, rho, v_best_climb_angle_jet(wl, rho, k, c_d0), c_d0, k, gg, ee_prop, propulsion); + pl = climb(wl, rho, v_best_climb_angle_jet(wl, rho, k, cd_0), cd_0, k, gg, propulsion); % TODO: Replace with segment speed elseif is_prop(propulsion) - ptl = climb(wl, rho, v_best_climb_angle_prop(wl, rho, k, c_d0), c_d0, k, gg, ee_prop, propulsion); + pl = climb(wl, rho, v_best_climb_angle_prop(wl, rho, k, cd_0), cd_0, k, gg, propulsion); % TODO: Replace with segment speed end -% function ptl = climb_rate(wl, rho, c_d0, k, gg, ee_prop, propulsion) +% function pl = climb_rate(wl, rho, cd_0, k, gg, propulsion) % if is_jet(propulsion) -% % ptl = fsolve(@(x)climb_rate_jet_error(x, wl, rho, c_d0, k, gg, ee_prop, propulsion), 0.01, optimoptions('fsolve', 'Display','iter')); +% % pl = fsolve(@(x)climb_rate_jet_error(x, wl, rho, cd_0, k, gg, propulsion), 0.01, optimoptions('fsolve', 'Display','iter')); % elseif is_prop(propulsion) -% ptl = climb(wl, rho, v_best_climb_rate_prop(wl, rho, k, c_d0), c_d0, k, gg, ee_prop, propulsion); +% pl = climb(wl, rho, v_best_climb_rate_prop(wl, rho, k, cd_0), cd_0, k, gg, propulsion); % end -% function err = climb_rate_jet_error(tl, wl, rho, c_d0, k, gg, ee_prop, propulsion) -% err = climb(wl, rho, v_best_climb_rate_jet(tl, wl, rho, k, c_d0), c_d0, k, gg, ee_prop, propulsion) - tl; +% function err = climb_rate_jet_error(tl, wl, rho, cd_0, k, gg, propulsion) +% err = climb(wl, rho, v_best_climb_rate_jet(tl, wl, rho, k, cd_0), cd_0, k, gg, propulsion) - tl; + +% Forward flight constraint regions +function c = range_region(wl, rho, v, cd_0, k, propulsion) +if is_jet(propulsion) + c = wl < 0.5 * rho * v^2 * sqrt(cd_0 / 3 / k); +elseif is_prop(propulsion) + c = wl < 0.5 * rho * v^2 * sqrt(cd_0 / k); +end + +function c = endurance_region(wl, rho, v, cd_0, k, propulsion) +if is_jet(propulsion) + c = wl < 0.5 * rho * v^2 * sqrt(cd_0 / k); +elseif is_prop(propulsion) + c = wl < 0.5 * rho * v^2 * sqrt(3 * cd_0 / k); +end + +function c = stall_speed_region(rho, v_s, c_lmax) +c = wl < 0.5 * rho * v_s^2 * c_lmax; + +function c = cruise_speed_region(pl, wl, rho, v, cd_0, k, propulsion) +if is_jet(propulsion) + c = pl < 1 ./ (rho .* v.^2 .* cd_0 ./ 2 ./ wl + 2 .* k .* wl ./ rho ./ v.^2); +elseif is_prop(propulsion) + c = pl < propulsion.efficiency ./ (rho .* v.^3 .* cd_0 ./ 2 ./ wl + 2 .* k .* wl ./ rho ./ v); +end + +function c = climb_region(pl, wl, rho, v, cd_0, k, gg, propulsion) +if is_jet(propulsion) + c = pl < 1 ./ (sind(gg) + rho .* v.^2 .* cd_0 ./ 2 ./ wl + 2 .* k .* wl ./ rho ./ v.^2); +elseif is_prop(propulsion) + c = pl < propulsion.efficiency ./ (v .* sind(gg) + rho .* v.^3 .* cd_0 ./ 2 ./ wl + 2 .* k .* wl ./ rho ./ v); +end + +function c = climb_angle_region(pl, wl, rho, cd_0, k, gg, propulsion) +if is_jet(propulsion) + c = climb_region(pl, wl, rho, v_best_climb_angle_jet(wl, rho, k, cd_0), cd_0, k, gg, propulsion); +elseif is_prop(propulsion) + c = climb_region(pl, wl, rho, v_best_climb_angle_prop(wl, rho, k, cd_0), cd_0, k, gg, propulsion); +end + +function i = find_energy_source_index(aircraft, energy_source) +for i = 1 : length(aircraft.energy_sources) + if strcmp(aircraft.energy_sources{i}.name, energy_source) + break; + end +end + +function i = find_propulsion_type_index(aircraft, propulsion_type) +for i = 1 : length(aircraft.propulsion_types) + if strcmp(aircraft.propulsion_types{i}.name, propulsion_type) + break; + end +end diff --git a/drag_buildup.m b/drag_buildup.m index 0e8a244..e21f17f 100644 --- a/drag_buildup.m +++ b/drag_buildup.m @@ -4,18 +4,43 @@ % % This file is subject to the license terms in the LICENSE file included in this distribution -function aircraft = drag_buildup(aircraft) +function aircraft = drag_buildup(aircraft, mission) -aircraft.performance.c_d0 = 0; +% Iterate over aircraft components +aircraft.performance.cd_0 = zeros(length(mission.segments), 1); +for i = 1 : length(mission.segments) + for j = 1 : length(aircraft.components) + aircraft.components{j}.cd_0 = zeros(length(mission.segments), 1); + if is_fuselage(aircraft.components{j}) + aircraft.components{j}.cd_0(i) = friction_coeff(aircraft.components{j}.length, mission.segments{i}.velocity, mission.segments{i}.speed_sound, mission.segments{i}.density, air_viscosity(mission.segments{i}.temperature)) *... + fuselage_form_factor(aircraft.components{j}.length / aircraft.components{j}.diameter) *... + aircraft.components{j}.interf_factor *... + aircraft.components{j}.area_wet / aircraft.components{j}.area_ref; + elseif is_wing(aircraft.components{j}) + m = mean(mission.segments{i}.velocity) / mean(mission.segments{i}.speed_sound); + bb = sqrt(1 - m^2); + aircraft.components{j}.cl_aa = aircraft.components{j}.airfoil.cl_aa * aircraft.components{j}.aspect_ratio /... + (2 + sqrt(4 + aircraft.components{j}.aspect_ratio^2 * bb^2 * (1 + tand(aircraft.components{j}.sweep_tc_max)^2 / bb^2))); -% Iterate over aircraft lifting surfaces -for i = 1 : length(aircraft.lifting_surfaces) - aircraft.performance.c_d0 = aircraft.performance.c_d0 + aircraft.lifting_surfaces{i}.c_d0 * aircraft.lifting_surfaces{i}.area_ref; -end + aircraft.components{j}.cd_0 = friction_coeff(aircraft.components{j}.mean_chord, mission.segments{i}.velocity, mission.segments{i}.speed_sound, mission.segments{i}.density, air_viscosity(mission.segments{i}.temperature)) *... + wing_form_factor(aircraft.components{j}.airfoil.xc_max, aircraft.components{j}.airfoil.tc_max, aircraft.components{j}.sweep_tc_max, m) *... + aircraft.components{j}.interf_factor *... + aircraft.components{j}.area_wet / aircraft.components{j}.area_ref; + end -% Iterate over aircraft fuselages -for i = 1 : length(aircraft.fuselages) - aircraft.performance.c_d0 = aircraft.performance.c_d0 + aircraft.fuselages{i}.c_d0 * aircraft.fuselages{i}.area_ref; + aircraft.performance.cd_0 = aircraft.performance.cd_0 + aircraft.components{j}.cd_0 .* aircraft.components{j}.area_ref; + end end +aircraft.performance.cd_0 = aircraft.performance.cd_0 ./ aircraft.performance.area_ref; + +function f = fuselage_form_factor(ld) +f = 1 + 60 / ld^3 + ld / 400; + +function f = wing_form_factor(xc_max, tc_max, sweep_tc_max, m) +f = (1 + 0.6 / xc_max * tc_max + 100 * tc_max^4) * 1.34 * m^0.18 * cosd(sweep_tc_max)^0.28; + +function test = is_fuselage(component) +test = strcmp(component.type, 'fuselage'); -aircraft.performance.c_d0 = aircraft.performance.c_d0 / aircraft.performance.area_ref; +function test = is_wing(component) +test = strcmp(component.type, 'wing'); diff --git a/estimate_mf_prop.m b/estimate_mf_prop.m deleted file mode 100644 index a25c58f..0000000 --- a/estimate_mf_prop.m +++ /dev/null @@ -1,19 +0,0 @@ -% Aircraft design tool -% -% Mario Bras (mbras@uvic.ca) and Ricardo Marques (ricardoemarques@uvic.ca) 2019 -% -% This file is subject to the license terms in the LICENSE file included in this distribution - -function mf_prop = estimate_mf_prop(f_install, pw_to, pw_ice, pw_gen, pw_em, pw_prop, config) -% f_install: Factor that accounts for all auxiliary systems for propulsion (> 1) -% pw_to: Shaft power to MTOW ratio -% pw_ice: Maximum shaft power to ICE weight ratio -% pw_gen: Maximum shaft power to generator weight ratio -% pw_em: Maximum shaft power to electric motor weight ratio -% pw_prop: Maximum shaft power to propeller weight ratio - -if strcmp(config, 'series') - mf_prop = f_install * pw_to / (pw_ice + pw_gen + pw_em + pw_prop); -elseif strcmp(config, 'parallel') - mf_prop = f_install * pw_to / (pw_ice + pw_em + pw_prop); -end diff --git a/estimate_mf_subs.m b/estimate_mf_subs.m deleted file mode 100644 index b5d8a34..0000000 --- a/estimate_mf_subs.m +++ /dev/null @@ -1,9 +0,0 @@ -% Aircraft design tool -% -% Mario Bras (mbras@uvic.ca) and Ricardo Marques (ricardoemarques@uvic.ca) 2019 -% -% This file is subject to the license terms in the LICENSE file included in this distribution - -function mf_subs = estimate_mf_subs(mf_prop, mf_struct) - -mf_subs = 3 / 7 * (mf_prop + mf_struct); diff --git a/examples/project.json b/examples/project.json index b898b5a..ee9c315 100644 --- a/examples/project.json +++ b/examples/project.json @@ -2,7 +2,7 @@ "$schema": "https://raw.githubusercontent.com/marioarbras/aircraft-design-tool/master/schemas/project.json", "concept": { "categories": { - "name": "Awesome New Concept", + "name": "New Concept", "pairs": [ [1, 0.14286], [7, 1] @@ -56,162 +56,205 @@ ] }, "mission": { - "mf_struct": 0.24, - "mf_subs": 0.0, "segments": [ { "type": "taxi", - "energy": { - "type": "electric" - }, - "time": 600.0, + "energy_source": "Main Battery", + "time": 120.0, "altitude": 0.0 }, { "type": "vertical_climb", - "energy": { - "type": "electric" - }, - "velocity": 20.0, - "altitude": [ 0.0, 1000.0] + "propulsion_type": "Prop Engine #2", + "energy_source": "Main Battery", + "velocity": 8.0, + "altitude": [0.0, 500.0] + }, + { + "type": "hover", + "propulsion_type": "Prop Engine #2", + "energy_source": "Main Battery", + "altitude": 500.0, + "time": 10.0 + }, + { + "type": "transition", + "propulsion_type": "Prop Engine #2", + "energy_source": "Main Battery", + "altitude": 500.0, + "transition_angle": 40.0, + "time": 120.0, + "velocity": [0.0, 40.0] + }, + { + "type": "climb", + "propulsion_type": "Prop Engine #1", + "energy_source": "Main Fuel Tank", + "velocity": 40.0, + "altitude": [500.0, 2500.0], + "angle": 7.2 }, { "type": "cruise", - "propulsion": { - "type": "prop" - }, - "energy": { - "type": "electric" - }, - "velocity": 100.0, - "range": 5000.0, - "altitude": 1000.0 + "propulsion_type": "Prop Engine #1", + "energy_source": "Main Battery", + "velocity": 50.0, + "range": 80000.0, + "altitude": 2500.0 }, { "type": "descent", - "energy": { - "type": "electric" - }, - "velocity": -60.0, - "altitude": [ 1000.0, 100.0], - "angle": -20.0 + "energy_source": "Main Battery", + "velocity": -40.0, + "altitude": [2500.0, 1000.0], + "angle": -5.0 }, { "type": "drop", "mass_drop": 1000.0, "time": 0.0, - "altitude": 100.0 + "altitude": 1000.0 }, { "type": "hold", - "propulsion": { - "type": "prop" - }, - "energy": { - "type": "fuel", - "sfc": 4.25e-5 - }, + "propulsion_type": "High Efficiency Prop Engine", + "energy_source": "Fuel Tank #2", "velocity": 40.0, "time": 300.0, - "altitude": 100.0 + "altitude": 1000.0 }, { "type": "load", "mass_reload": 1000.0, "time": 0.0, - "altitude": 100.0 + "altitude": 1000.0 }, { "type": "climb", - "propulsion": { - "type": "prop" - }, - "energy": { - "type": "fuel" - }, - "velocity": 50.0, - "altitude": [100.0, 1000.0], - "angle": 20.0 + "propulsion_type": "Prop Engine #1", + "energy_source": "Main Fuel Tank", + "velocity": 40.0, + "altitude": [1000.0, 2500.0], + "angle": 7.2 }, { "type": "cruise", - "propulsion": { - "type": "jet" - }, - "energy": { - "type": "fuel", - "sfc": 4.25e-5 - }, - "velocity": 100.0, - "range": 5000.0, - "altitude": 1000.0 + "propulsion_type": "Jet Engine", + "energy_source": "Main Fuel Tank", + "velocity": 50.0, + "range": 80000.0, + "altitude": 2500.0 + }, + { + "type": "descent", + "energy_source": "Prop Engine #1", + "velocity": -6.0, + "altitude": [2500.0, 400.0], + "angle": -5.0 }, { "type": "vertical_descent", - "energy": { - "type": "electric" - }, - "velocity": -50, - "altitude": [1000.0, 0.0] + "propulsion_type": "Prop Engine #2", + "energy_source": "Battery #2", + "velocity": -6.0, + "altitude": [400.0, 0.0] }, { "type": "landing", - "energy": { - "type": "electric" - }, - "time": 600.0, + "energy_source": "Main Battery", + "time": 120.0, "altitude": 0.0 } ] }, "aircraft": { - "mass_to": 4000.0, - "mass_payload": 680.0, - "propulsion": { - "gear_efficiency": 0.6, - "em_efficiency": 0.9, - "esc_efficiency": 0.8, - "config": "series", - "dist_efficiency": 0.6, - "prop_efficiency": 0.8, - "k_i": 1.15, - "c_d": 0.02, - "ss": 0.08, - "tip_velocity": 240.1, - "fm": 0.6, - "a": 10.0, - "n": 4.0, - "tt_tilt": 40.0, - "fwd_power": 500000.0, - "vert_power": 500000.0, - "mass": 500.0 + "payload": { + "mass": 680.0 }, - "energy": { - "e_spec_batt": 360000.0, - "batt_efficiency": 0.9, - "f_usable_batt": 0.9, - "reserve_batt": 1.2, - "reserve_fuel": 1.06 + "structure": { + "mass": 2000.0 }, + "subsystems": { + "mass": 50.0 + }, + "propulsion_types": [ + { + "name": "High Efficiency Prop Engine", + "type": "prop", + "efficiency": 0.8, + "specific_fuel_consumption": 4.25e-5, + "mass": 200.0 + }, + { + "name": "Prop Engine #1", + "type": "prop", + "efficiency": 0.6, + "specific_fuel_consumption": 4.25e-5, + "mass": 200.0 + }, + { + "name": "Jet Engine", + "type": "jet", + "specific_fuel_consumption": 4.25e-5, + "mass": 300.0 + }, + { + "name": "Prop Engine #2", + "type": "prop", + "mass": 500.0, + "tip_velocity": 240.1, + "efficiency": 0.6, + "specific_fuel_consumption": 4.25e-5, + "ss": 0.08, + "cd": 0.02, + "k_i": 1.15, + "fm": 0.6 + } + ], + "energy_sources": [ + { + "name": "Main Battery", + "type": "electric", + "specific_energy": 360000.0, + "efficiency": 0.9, + "reserve": 0.2 + }, + { + "name": "Battery #2", + "type": "electric", + "specific_energy": 360000.0, + "efficiency": 0.9 + }, + { + "name": "Main Fuel Tank", + "type": "fuel", + "specific_energy": 0.123, + "reserve": 0.06 + }, + { + "name": "Fuel Tank #2", + "type": "fuel", + "specific_energy": 0.123 + } + ], "performance": { - "c_d0": 0.02, + "cd_0": 0.02, "aspect_ratio": 5.0, - "ld_max": 15.0, - "dl": 700.0, - "area_ref": 15.0, - "area_wet": 30.0 + "ld_max": 8.0, + "disc_loading": 900 }, - "fuselages": [ + "components": [ { + "name": "Fuselage", + "type": "fuselage", "interf_factor": 1.0, "diameter": 1.0, "length": 4.0, "area_ref": 1.0, "area_wet": 2.0 - } - ], - "lifting_surfaces": [ + }, { + "name": "Main Wing", "type": "wing", "interf_factor": 1.0, "aspect_ratio": 5.0, @@ -231,7 +274,8 @@ "area_wet": 30.0 }, { - "type": "h-tail", + "name": "Horizontal Tail", + "type": "wing", "interf_factor": 1.0, "aspect_ratio": 5.0, "span": 2.0, @@ -250,7 +294,8 @@ "area_wet": 2.0 }, { - "type": "v-tail", + "name": "Vertical Tail", + "type": "wing", "interf_factor": 1.0, "aspect_ratio": 5.0, "span": 2.0, diff --git a/friction_coeff.m b/friction_coeff.m index 025a69a..8f6deb3 100644 --- a/friction_coeff.m +++ b/friction_coeff.m @@ -6,7 +6,7 @@ function cf = friction_coeff(x, v, a, rr, mm) - re = reynolds(x, v, rr, mm); +re = reynolds(x, v, rr, mm); if is_laminar(re) cf = 1.328 / sqrt(re); else diff --git a/fuselage.m b/fuselage.m deleted file mode 100644 index 5ca7630..0000000 --- a/fuselage.m +++ /dev/null @@ -1,18 +0,0 @@ -% Aircraft design tool -% -% Mario Bras (mbras@uvic.ca) and Ricardo Marques (ricardoemarques@uvic.ca) 2019 -% -% This file is subject to the license terms in the LICENSE file included in this distribution - -function aircraft = fuselage(aircraft, segment) - -% Iterate over aircraft lifting surfaces -for i = 1 : length(aircraft.fuselages) - aircraft.fuselages{i}.c_d0 = friction_coeff(aircraft.fuselages{i}.length, segment.velocity, segment.speed_sound, segment.density, air_viscosity(segment.temperature)) *... - form_factor(aircraft.fuselages{i}.length / aircraft.fuselages{i}.diameter) *... - aircraft.fuselages{i}.interf_factor *... - aircraft.fuselages{i}.area_wet / aircraft.fuselages{i}.area_ref; -end - -function f = form_factor(ld) -f = 1 + 60 / ld^3 + ld / 400; diff --git a/lift_slope.m b/lift_slope.m new file mode 100644 index 0000000..ddc4fb4 --- /dev/null +++ b/lift_slope.m @@ -0,0 +1,29 @@ +% Aircraft design tool +% +% Mario Bras (mbras@uvic.ca) and Ricardo Marques (ricardoemarques@uvic.ca) 2019 +% +% This file is subject to the license terms in the LICENSE file included in this distribution + +function aircraft = lift_slope(aircraft, mission) + +% Iterate over aircraft lifting surfaces +for i = 1 : length(mission.segments) + for j = 1 : length(aircraft.components) + aircraft.components{j}.cl_aa = zeros(length(mission.segments), 1); + if is_wing(aircraft.components{j}) + m = mean(mission.segments{i}.velocity) / mean(mission.segments{i}.speed_sound); + bb = sqrt(1 - m^2); + aircraft.components{j}.cl_aa = aircraft.components{j}.airfoil.cl_aa * aircraft.components{j}.aspect_ratio /... + (2 + sqrt(4 + aircraft.components{j}.aspect_ratio^2 * bb^2 * (1 + tand(aircraft.components{j}.sweep_tc_max)^2 / bb^2))); + end + end +end + +function f = form_factor(xc_max, tc_max, sweep_tc_max, m) +f = (1 + 0.6 / xc_max * tc_max + 100 * tc_max^4) * 1.34 * m^0.18 * cosd(sweep_tc_max)^0.28; + +function test = is_fuselage(component) +test = strcmp(component.type, 'fuselage'); + +function test = is_wing(component) +test = strcmp(component.type, 'wing'); diff --git a/lifting_surface.m b/lifting_surface.m deleted file mode 100644 index 4ee5349..0000000 --- a/lifting_surface.m +++ /dev/null @@ -1,24 +0,0 @@ -% Aircraft design tool -% -% Mario Bras (mbras@uvic.ca) and Ricardo Marques (ricardoemarques@uvic.ca) 2019 -% -% This file is subject to the license terms in the LICENSE file included in this distribution - -function aircraft = lifting_surface(aircraft, segment) - -m = segment.velocity / segment.speed_sound; - -% Iterate over aircraft lifting surfaces -for i = 1 : length(aircraft.lifting_surfaces) - bb = sqrt(1 - m^2); - aircraft.lifting_surfaces{i}.cl_aa = aircraft.lifting_surfaces{i}.airfoil.cl_aa * aircraft.lifting_surfaces{i}.aspect_ratio /... - (2 + sqrt(4 + aircraft.lifting_surfaces{i}.aspect_ratio^2 * bb^2 * (1 + tand(aircraft.lifting_surfaces{i}.sweep_tc_max)^2 / bb^2))); - - aircraft.lifting_surfaces{i}.c_d0 = friction_coeff(aircraft.lifting_surfaces{i}.mean_chord, segment.velocity, segment.speed_sound, segment.density, air_viscosity(segment.temperature)) *... - form_factor(aircraft.lifting_surfaces{i}.airfoil.xc_max, aircraft.lifting_surfaces{i}.airfoil.tc_max, aircraft.lifting_surfaces{i}.sweep_tc_max, m) *... - aircraft.lifting_surfaces{i}.interf_factor *... - aircraft.lifting_surfaces{i}.area_wet / aircraft.lifting_surfaces{i}.area_ref; -end - -function f = form_factor(xc_max, tc_max, sweep_tc_max, m) -f = (1 + 0.6 / xc_max * tc_max + 100 * tc_max^4) * 1.34 * m^0.18 * cosd(sweep_tc_max)^0.28; diff --git a/load_project.m b/load_project.m index 63d9bcd..6364e10 100644 --- a/load_project.m +++ b/load_project.m @@ -29,11 +29,15 @@ % Aircraft if isfield(data, 'aircraft') - if ~iscell(data.aircraft.fuselages) - data.aircraft.fuselages = num2cell(data.aircraft.fuselages); + if ~iscell(data.aircraft.propulsion_types) + data.aircraft.propulsion_types = num2cell(data.aircraft.propulsion_types); + end + + if ~iscell(data.aircraft.energy_sources) + data.aircraft.energy_sources = num2cell(data.aircraft.energy_sources); end - if ~iscell(data.aircraft.lifting_surfaces) - data.aircraft.lifting_surfaces = num2cell(data.aircraft.lifting_surfaces); + if ~iscell(data.aircraft.components) + data.aircraft.components = num2cell(data.aircraft.components); end end diff --git a/mtow.m b/mtow.m index b540e51..5ee3b9b 100644 --- a/mtow.m +++ b/mtow.m @@ -5,203 +5,222 @@ % This file is subject to the license terms in the LICENSE file included in this distribution function [mission, aircraft] = mtow(mission, aircraft, constants) +mass_to = fsolve(@(x)mtow_error(x, mission, aircraft, constants), 0, optimoptions('fsolve', 'Display','none')); +[~, mission, aircraft] = mtow_error(mass_to, mission, aircraft, constants); -% Initialize mission parameters -mission.mf_batt = 0; -mission.mf_fuel = 0; +function i = find_energy_source_index(aircraft, energy_source) +for i = 1 : length(aircraft.energy_sources) + if strcmp(aircraft.energy_sources{i}.name, energy_source) + break; + end +end + +function i = find_propulsion_type_index(aircraft, propulsion_type) +for i = 1 : length(aircraft.propulsion_types) + if strcmp(aircraft.propulsion_types{i}.name, propulsion_type) + break; + end +end + +function test = is_fuel(energy_source) +test = strcmp(energy_source.type, 'fuel'); + +function test = is_electric(energy_source) +test = strcmp(energy_source.type, 'electric'); + +function test = is_jet(propulsion) +test = strcmp(propulsion.type, 'jet'); + +function test = is_prop(propulsion) +test = strcmp(propulsion.type, 'prop'); + +function ld = get_ld(performance, segment, propulsion_type) +if strcmp(segment.type, 'cruise') + if is_jet(propulsion_type) + ld = 0.886 * performance.ld_max; + elseif is_prop(propulsion_type) + ld = performance.ld_max; + end +elseif strcmp(segment.type, 'hold') + if is_jet(propulsion_type) + ld = performance.ld_max; + elseif is_prop(propulsion_type) + ld = 0.886 * performance.ld_max; + end +end + +function [err, mission, aircraft] = mtow_error(mass_to, mission, aircraft, constants) + +% Initialize mission and aircraft parameters mission.time = 0; mission.range = 0; -mf_to_minus_fuel = 1; +mf_batt = 0; +mf_fuel = 0; +mass_ac = mass_to; + +for i = 1 : length(aircraft.energy_sources) + aircraft.energy_sources{i}.mass = 0; +end % Iterate over mission segments for i = 1 : length(mission.segments) - mission.segments{i}.mf_batt = 0; - mission.segments{i}.mf_fuel = 0; - [mission.segments{i}.temperature, mission.segments{i}.speed_sound, mission.segments{i}.pressure, mission.segments{i}.density] = atmosisa(mission.segments{i}.altitude); if strcmp(mission.segments{i}.type, 'taxi') % Taxi segment - if is_fuel(mission.segments{i}.energy) - mission.segments{i}.mf_fuel = 1 - 0.9725; - elseif is_electric(mission.segments{i}.energy) - mission.segments{i}.mf_batt = 0; + ei = find_energy_source_index(aircraft, mission.segments{i}.energy_source); + if is_fuel(aircraft.energy_sources{ei}) + mf_fuel = 1 - 0.9725; + elseif is_electric(aircraft.energy_sources{ei}) + mf_batt = 0; end - mission.segments{i}.range = 0; elseif strcmp(mission.segments{i}.type, 'hover') % Hover segment - if is_fuel(mission.segments{i}.energy) + ei = find_energy_source_index(aircraft, mission.segments{i}.energy_source); + pti = find_propulsion_type_index(aircraft, mission.segments{i}.propulsion_type); + if is_fuel(aircraft.energy_sources{ei}) errordlg('Not available'); % NOT AVAILABLE break; - elseif is_electric(mission.segments{i}.energy) - pl = aircraft.propulsion.fm * sqrt(2 * mission.segments{i}.density / aircraft.performance.dl); - mission.segments{i}.mf_batt = mission.segments{i}.time * constants.g / aircraft.energy.e_spec_batt / aircraft.propulsion.total_efficiency / aircraft.energy.batt_efficiency / aircraft.energy.f_usable_batt / pl; + elseif is_electric(aircraft.energy_sources{ei}) + pl = aircraft.propulsion_types{pti}.fm * sqrt(2 * mission.segments{i}.density / aircraft.performance.disc_loading); + mf_batt = mission.segments{i}.time * constants.g / aircraft.energy_sources{ei}.specific_energy / aircraft.propulsion_types{pti}.efficiency / aircraft.energy_sources{ei}.efficiency / pl; end - mission.segments{i}.range = 0; elseif strcmp(mission.segments{i}.type, 'climb') % Climb segment - if is_fuel(mission.segments{i}.energy) + ei = find_energy_source_index(aircraft, mission.segments{i}.energy_source); + if is_fuel(aircraft.energy_sources{ei}) mach = mission.segments{i}.velocity / mission.segments{i}.speed_sound(1); if mach < 1 - mission.segments{i}.mf_fuel = 1 - (1 - 0.04 * mach); + mf_fuel = 1 - (1 - 0.04 * mach); else - mission.segments{i}.mf_fuel = 1 - (0.96 - 0.03 * (mach - 1)); + mf_fuel = 1 - (0.96 - 0.03 * (mach - 1)); end - elseif is_electric(mission.segments{i}.energy) + elseif is_electric(aircraft.energy_sources{ei}) errordlg('Not available'); % NOT AVAILABLE break; end - mission.segments{i}.time = (mission.segments{i}.altitude(2) - mission.segments{i}.altitude(1)) / mission.segments{i}.velocity / sind(mission.segments{i}.angle); - mission.segments{i}.range = mission.segments{i}.velocity * mission.segments{i}.time * cosd(mission.segments{i}.angle); elseif strcmp(mission.segments{i}.type, 'vertical_climb') % Vertical climb segment - altitude_range = (mission.segments{i}.altitude(2) - mission.segments{i}.altitude(1)); - if is_fuel(mission.segments{i}.energy) + altitude_range = mission.segments{i}.altitude(2) - mission.segments{i}.altitude(1); + ei = find_energy_source_index(aircraft, mission.segments{i}.energy_source); + pti = find_propulsion_type_index(aircraft, mission.segments{i}.propulsion_type); + if is_fuel(aircraft.energy_sources{ei}) errordlg('Not available'); % NOT AVAILABLE break; - elseif is_electric(mission.segments{i}.energy) - pl = 1 / (mission.segments{i}.velocity - aircraft.propulsion.k_i / 2 * mission.segments{i}.velocity + aircraft.propulsion.k_i / 2 * sqrt(mission.segments{i}.velocity^2 + 2 * aircraft.performance.dl / mission.segments{i}.density(1)) + mission.segments{i}.density(1) * aircraft.propulsion.tip_velocity^3 / aircraft.performance.dl * aircraft.propulsion.ss * aircraft.propulsion.c_d / 8); % Power loading - mission.segments{i}.mf_batt = altitude_range * constants.g / aircraft.energy.e_spec_batt / aircraft.propulsion.total_efficiency / aircraft.energy.batt_efficiency / aircraft.energy.f_usable_batt / pl / mission.segments{i}.velocity; % Mass fraction for this segment + elseif is_electric(aircraft.energy_sources{ei}) + pl = 1 / (mission.segments{i}.velocity - aircraft.propulsion_types{pti}.k_i / 2 * mission.segments{i}.velocity + aircraft.propulsion_types{pti}.k_i / 2 * sqrt(mission.segments{i}.velocity^2 + 2 * aircraft.performance.disc_loading / mission.segments{i}.density(1)) + mission.segments{i}.density(1) * aircraft.propulsion_types{pti}.tip_velocity^3 / aircraft.performance.disc_loading * aircraft.propulsion_types{pti}.ss * aircraft.propulsion_types{pti}.cd / 8); % Power loading + mf_batt = altitude_range * constants.g / aircraft.energy_sources{ei}.specific_energy / aircraft.propulsion_types{pti}.efficiency / aircraft.energy_sources{ei}.efficiency / pl / mission.segments{i}.velocity; % Mass fraction for this segment end - mission.segments{i}.time = altitude_range * mission.segments{i}.velocity; - mission.segments{i}.range = 0; elseif strcmp(mission.segments{i}.type, 'acceleration') % Acceleration segment mach = mission.segments{i}.velocity / mission.segments{i}.a; - if is_fuel(mission.segments{i}.energy) + ei = find_energy_source_index(aircraft, mission.segments{i}.energy_source); + if is_fuel(aircraft.energy_sources{ei}) if i > 1 && mach == mission.segment{i - 1}.velocity / mission.segments{i - 1}.a - mission.segments{i}.mf_fuel = 0; + mf_fuel = 0; elseif mach < 1 - mission.segments{i}.mf_fuel = 1 - (1 - 0.04 * mach); + mf_fuel = 1 - (1 - 0.04 * mach); else - mission.segments{i}.mf_fuel = 1 - (0.96 - 0.03 * (mach - 1)); + mf_fuel = 1 - (0.96 - 0.03 * (mach - 1)); end - elseif is_electric(mission.segments{i}.energy) + elseif is_electric(aircraft.energy_sources{ei}) errordlg('Not available'); % NOT AVAILABLE break; end - mission.segments{i}.range = mission.segments{i}.velocity * mission.segments{i}.time; elseif strcmp(mission.segments{i}.type, 'cruise') % Cruise segment - ld = get_ld(aircraft.performance,mission.segments{i}); - - if is_fuel(mission.segments{i}.energy) - mission.segments{i}.mf_fuel = 1 - exp(-mission.segments{i}.range * mission.segments{i}.energy.sfc / mission.segments{i}.velocity / ld); - elseif is_electric(mission.segments{i}.energy) - mission.segments{i}.mf_batt = mission.segments{i}.range * constants.g / aircraft.energy.e_spec_batt / aircraft.propulsion.total_efficiency / aircraft.energy.batt_efficiency / aircraft.energy.f_usable_batt / ld; + pti = find_propulsion_type_index(aircraft, mission.segments{i}.propulsion_type); + ld = get_ld(aircraft.performance, mission.segments{i}, aircraft.propulsion_types{pti}); + ei = find_energy_source_index(aircraft, mission.segments{i}.energy_source); + + if is_fuel(aircraft.energy_sources{ei}) + mf_fuel = 1 - exp(-mission.segments{i}.range * aircraft.propulsion_types{pti}.specific_fuel_consumption / mission.segments{i}.velocity / ld); + elseif is_electric(aircraft.energy_sources{ei}) + mf_batt = mission.segments{i}.range * constants.g / aircraft.energy_sources{ei}.specific_energy / aircraft.propulsion_types{pti}.efficiency / aircraft.energy_sources{ei}.efficiency / ld; end - mission.segments{i}.time = mission.segments{i}.range * mission.segments{i}.velocity; elseif strcmp(mission.segments{i}.type, 'hold') % Hold segment - ld = get_ld(aircraft.performance, mission.segments{i}); - - if is_fuel(mission.segments{i}.energy) - mission.segments{i}.mf_fuel = 1 - exp(-mission.segments{i}.time * mission.segments{i}.energy.sfc / ld); - elseif is_electric(mission.segments{i}.energy) - mission.segments{i}.mf_batt = mission.segments{i}.time * mission.segments{i}.velocity * constants.g / aircraft.energy.e_spec_batt / aircraft.propulsion.total_efficiency / aircraft.energy.batt_efficiency / aircraft.energy.f_usable_batt / ld; + pti = find_propulsion_type_index(aircraft, mission.segments{i}.propulsion_type); + ld = get_ld(aircraft.performance, mission.segments{i}, aircraft.propulsion_types{pti}); + ei = find_energy_source_index(aircraft, mission.segments{i}.energy_source); + if is_fuel(aircraft.energy_sources{ei}) + mf_fuel = 1 - exp(-mission.segments{i}.time * aircraft.propulsion_types{pti}.specific_fuel_consumption / ld); + elseif is_electric(aircraft.energy_sources{ei}) + mf_batt = mission.segments{i}.time * mission.segments{i}.velocity * constants.g / aircraft.energy_sources{ei}.specific_energy / aircraft.propulsion_types{pti}.efficiency / aircraft.energy_sources{ei}.efficiency / ld; end - mission.segments{i}.range = mission.segments{i}.velocity * mission.segments{i}.time; - % elseif strcmp(mission.segments{i}.type, 'combat') % Combat segment - % if is_fuel(mission.segments{i}.energy) - % % mission.segments{i}.mf_fuel = 1 - % TODO - % elseif is_electric(mission.segments{i}.energy) + elseif strcmp(mission.segments{i}.type, 'combat') % Combat segment + % ei = find_energy_source_index(aircraft, mission.segments{i}.energy_source); + % if is_fuel(aircraft.energy_sources{ei}) + % % mf_fuel = 1 - % TODO + % elseif is_electric(aircraft.energy_sources{ei}) % errordlg('Not available'); % NOT AVAILABLE % break; % end elseif strcmp(mission.segments{i}.type, 'descent') % Descent segment - if is_fuel(mission.segments{i}.energy) - mission.segments{i}.mf_fuel = 0; - elseif is_electric(mission.segments{i}.energy) - mission.segments{i}.mf_batt = 0; + ei = find_energy_source_index(aircraft, mission.segments{i}.energy_source); + if is_fuel(aircraft.energy_sources{ei}) + mf_fuel = 0; + elseif is_electric(aircraft.energy_sources{ei}) + mf_batt = 0; end - mission.segments{i}.time = abs(mission.segments{i}.altitude(2) - mission.segments{i}.altitude(1)) / abs(mission.segments{i}.velocity) / abs(sind(mission.segments{i}.angle)); - mission.segments{i}.range = abs(mission.segments{i}.velocity) * mission.segments{i}.time * abs(cosd(mission.segments{i}.angle)); elseif strcmp(mission.segments{i}.type, 'vertical_descent') % Vertical descent segment altitude_range = abs(mission.segments{i}.altitude(2) - mission.segments{i}.altitude(1)); - if is_fuel(mission.segments{i}.energy) + ei = find_energy_source_index(aircraft, mission.segments{i}.energy_source); + pti = find_propulsion_type_index(aircraft, mission.segments{i}.propulsion_type); + if is_fuel(aircraft.energy_sources{ei}) errordlg('Not available'); % NOT AVAILABLE break; - elseif is_electric(mission.segments{i}.energy) - v_i = sqrt(aircraft.performance.dl / 2 / mission.segments{i}.density(2)); % Induced velocity in hover + elseif is_electric(aircraft.energy_sources{ei}) + v_i = sqrt(aircraft.performance.disc_loading / 2 / mission.segments{i}.density(2)); % Induced velocity in hover if mission.segments{i}.velocity / v_i <= -2 % If this condition is met, the vertical climb equation is used for descent, else, an empirical equation is employed - pl = 1 / (mission.segments{i}.velocity - aircraft.propulsion.k_i / 2 * (mission.segments{i}.velocity + sqrt(mission.segments{i}.velocity^2 - 2 * aircraft.performance.dl / mission.segments{i}.density(2))) + mission.segments{i}.density(2) * aircraft.propulsion.tip_velocity^3 / aircraft.performance.dl * aircraft.propulsion.ss * aircraft.propulsion.c_d / 8); + pl = 1 / (mission.segments{i}.velocity - aircraft.propulsion_types{pti}.k_i / 2 * (mission.segments{i}.velocity + sqrt(mission.segments{i}.velocity^2 - 2 * aircraft.performance.disc_loading / mission.segments{i}.density(2))) + mission.segments{i}.density(2) * aircraft.propulsion_types{pti}.tip_velocity^3 / aircraft.performance.disc_loading * aircraft.propulsion_types{pti}.ss * aircraft.propulsion_types{pti}.cd / 8); else - v_d = v_i * (aircraft.propulsion.k_i - 1.125 * mission.segments{i}.velocity / v_i - 1.372 * (mission.segments{i}.velocity / v_i)^2 - 1.718 * (mission.segments{i}.velocity / v_i)^3 - 0.655 * (mission.segments{i}.velocity / v_i)^4); % Induced velocity in descent according to an empirical relation (see lecture slides) - pl = 1 / (mission.segments{i}.velocity + aircraft.propulsion.k_i * v_d + mission.segments{i}.density(2) * aircraft.propulsion.tip_velocity^3 / aircraft.performance.dl * aircraft.propulsion.ss * aircraft.propulsion.c_d / 8); + v_d = v_i * (aircraft.propulsion_types{pti}.k_i - 1.125 * mission.segments{i}.velocity / v_i - 1.372 * (mission.segments{i}.velocity / v_i)^2 - 1.718 * (mission.segments{i}.velocity / v_i)^3 - 0.655 * (mission.segments{i}.velocity / v_i)^4); % Induced velocity in descent according to an empirical relation (see lecture slides) + pl = 1 / (mission.segments{i}.velocity + aircraft.propulsion_types{pti}.k_i * v_d + mission.segments{i}.density(2) * aircraft.propulsion_types{pti}.tip_velocity^3 / aircraft.performance.disc_loading * aircraft.propulsion_types{pti}.ss * aircraft.propulsion_types{pti}.cd / 8); end if pl > 0 - mission.segments{i}.mf_batt = altitude_range * constants.g / aircraft.energy.e_spec_batt / aircraft.propulsion.total_efficiency / aircraft.energy.batt_efficiency / aircraft.energy.f_usable_batt / pl / abs(mission.segments{i}.velocity); + mf_batt = altitude_range * constants.g / aircraft.energy_sources{ei}.specific_energy / aircraft.propulsion_types{pti}.efficiency / aircraft.energy_sources{ei}.efficiency / pl / abs(mission.segments{i}.velocity); else - mission.segments{i}.mf_batt = 0; + mf_batt = 0; end end - mission.segments{i}.time = altitude_range * abs(mission.segments{i}.velocity); - mission.segments{i}.range = 0; elseif strcmp(mission.segments{i}.type, 'landing') % Landing segment - if is_fuel(mission.segments{i}.energy) - mission.segments{i}.mf_fuel = 1 - 0.9725; - elseif is_electric(mission.segments{i}.energy) - mission.segments{i}.mf_batt = 0; + ei = find_energy_source_index(aircraft, mission.segments{i}.energy_source); + if is_fuel(aircraft.energy_sources{ei}) + mf_fuel = 1 - 0.9725; + elseif is_electric(aircraft.energy_sources{ei}) + mf_batt = 0; end - mission.segments{i}.range = 0; elseif strcmp(mission.segments{i}.type, 'drop') % Drop segment - mission.segments{i}.range = 0; + mass_ac = mass_ac - mission.segments{i}.mass_drop; elseif strcmp(mission.segments{i}.type, 'load') % Load segment - mission.segments{i}.range = 0; + mass_ac = mass_ac + mission.segments{i}.mass_reload; end % Accumulate mission parameters - mission.mf_batt = mission.mf_batt + mission.segments{i}.mf_batt; - mf_to_minus_fuel = mf_to_minus_fuel * (1 - mission.segments{i}.mf_fuel); + if isfield(mission.segments{i}, "energy_source") + ei = find_energy_source_index(aircraft, mission.segments{i}.energy_source); + if is_fuel(aircraft.energy_sources{ei}) + aircraft.energy_sources{ei}.mass = aircraft.energy_sources{ei}.mass + mf_fuel * mass_ac; + mass_ac = mass_ac - aircraft.energy_sources{ei}.mass; % Aircraft mass at end of segment + elseif is_electric(aircraft.energy_sources{ei}) + aircraft.energy_sources{ei}.mass = aircraft.energy_sources{ei}.mass + mf_batt * mass_to; + end + end mission.time = mission.time + mission.segments{i}.time; mission.range = mission.range + mission.segments{i}.range; end -mission.mf_fuel = 1 - mf_to_minus_fuel; - -% Iterate over drop and reload mission segments -for i = 1 : length(mission.segments) - if strcmp(mission.segments{i}.type, 'drop') % Drop segment - mission.mf_fuel = mission.mf_fuel - accumulate_mf_fuel(mission, i, length(mission.segments)) * mission.segments{i}.mass_drop / aircraft.mass_to; - end - if strcmp(mission.segments{i}.type, 'load') % Load segment - mission.mf_fuel = mission.mf_fuel + accumulate_mf_fuel(mission, i, length(mission.segments)) * mission.segments{i}.mass_reload / aircraft.mass_to; +% Accumulate energy mass +energy_mass = 0; +for i = 1 : length(aircraft.energy_sources) + energy_mass = energy_mass + aircraft.energy_sources{i}.mass; + if isfield(aircraft.energy_sources{i}, 'reserve') + energy_mass = energy_mass + aircraft.energy_sources{i}.mass * aircraft.energy_sources{i}.reserve; end end -mission.mf_batt = aircraft.energy.reserve_batt * mission.mf_batt; -mission.mf_fuel = aircraft.energy.reserve_fuel * mission.mf_fuel; -aircraft.mass_to = fsolve(@(x)mtow_error(x, mission, aircraft), aircraft.mass_to, optimoptions('fsolve', 'Display','none')); % aircraft.mass_to = aircraft.mass_payload / (1 - mission.mf_struct - mission.mf_subs - mission.mf_prop - mission.mf_fuel - mission.mf_batt); -aircraft.weight_to = aircraft.mass_to * constants.g; - -function test = is_fuel(energy) -test = strcmp(energy.type, 'fuel'); - -function test = is_electric(energy) -test = strcmp(energy.type, 'electric'); - -function test = is_jet(propulsion) -test = strcmp(propulsion.type, 'jet'); - -function test = is_prop(propulsion) -test = strcmp(propulsion.type, 'prop'); - -function mf_fuel = accumulate_mf_fuel(mission, first, last) -mf_fuel = 1; -for i = first : last - mf_fuel = mf_fuel * (1 - mission.segments{i}.mf_fuel); +% Accumulate propulsion mass +propulsion_mass = 0; +for i = 1 : length(aircraft.propulsion_types) + propulsion_mass = propulsion_mass + aircraft.propulsion_types{i}.mass; end -mf_fuel = 1 - mf_fuel; -function ld = get_ld(performance, segment) -if strcmp(segment.type, 'cruise') - if is_jet(segment.propulsion) - ld = 0.886 * performance.ld_max; - elseif is_prop(segment.propulsion) - ld = performance.ld_max; - end -elseif strcmp(segment.type, 'hold') - if is_jet(segment.propulsion) - ld = performance.ld_max; - elseif is_prop(segment.propulsion) - ld = 0.886 * performance.ld_max; - end -end +aircraft.mass = mass_to; -function err = mtow_error(mass_to, mission, aircraft) -err = aircraft.mass_payload / (1 - mission.mf_struct - mission.mf_subs - mission.mf_prop - mission.mf_fuel - mission.mf_batt) - mass_to; \ No newline at end of file +err = aircraft.mass - energy_mass - propulsion_mass - aircraft.payload.mass - aircraft.structure.mass - aircraft.subsystems.mass; diff --git a/plot_mission.m b/plot_mission.m index 8d5f9e1..f11d2c8 100644 --- a/plot_mission.m +++ b/plot_mission.m @@ -6,39 +6,64 @@ function plot_mission(mission) -x = 0; - figure; legend; hold on; +a = gca; +a.XLabel.String = 'Range'; +a.YLabel.String = "Altitude"; % Iterate over mission segments +x = 0; for i = 1 : length(mission.segments) if strcmp(mission.segments{i}.type, 'taxi') % Taxi segment - + plot([x, x + 500], [mission.segments{i}.altitude, mission.segments{i}.altitude], 'DisplayName', ['Taxi (', num2str(i-1), '-', num2str(i), ')']); + text([x, x + 500], [mission.segments{i}.altitude, mission.segments{i}.altitude], {num2str(i-1), num2str(i)}); + x = x + 500; elseif strcmp(mission.segments{i}.type, 'hover') % Hover segment - + plot([x, x + 500], [mission.segments{i}.altitude, mission.segments{i}.altitude], 'DisplayName', ['Hover (', num2str(i-1), '-', num2str(i), ')']); + text(x + 500, mission.segments{i}.altitude, num2str(i)); + x = x + 500; + elseif strcmp(mission.segments{i}.type, 'transition') % Transition segment + plot([x, x + 500], [mission.segments{i}.altitude, mission.segments{i}.altitude], 'DisplayName', ['Transition (', num2str(i-1), '-', num2str(i), ')']); + text(x + 500, mission.segments{i}.altitude, num2str(i)); + x = x + 500; elseif strcmp(mission.segments{i}.type, 'climb') % Climb segment - plot([x, x + mission.segments{i}.range], [mission.segments{i}.altitude(1), mission.segments{i}.altitude(2)], 'DisplayName', 'Climb'); + plot([x, x + mission.segments{i}.range], [mission.segments{i}.altitude(1), mission.segments{i}.altitude(2)], 'DisplayName', ['Climb (', num2str(i-1), '-', num2str(i), ')']); + text(x + mission.segments{i}.range, mission.segments{i}.altitude(2), num2str(i)); x = x + mission.segments{i}.range; elseif strcmp(mission.segments{i}.type, 'vertical_climb') % Vertical climb segment - plot([x, x + mission.segments{i}.range], [mission.segments{i}.altitude(1), mission.segments{i}.altitude(2)], 'DisplayName', 'Vertical descent'); + plot([x, x + mission.segments{i}.range], [mission.segments{i}.altitude(1), mission.segments{i}.altitude(2)], 'DisplayName', ['Vertical Climb (', num2str(i-1), '-', num2str(i), ')']); + text(x + mission.segments{i}.range, mission.segments{i}.altitude(2), num2str(i)); x = x + mission.segments{i}.range; elseif strcmp(mission.segments{i}.type, 'acceleration') % Acceleration segment - elseif strcmp(mission.segments{i}.type, 'cruise') % Cruise segment - plot([x, x + mission.segments{i}.range], [mission.segments{i}.altitude, mission.segments{i}.altitude], 'DisplayName', 'Cruise'); + plot([x, x + mission.segments{i}.range], [mission.segments{i}.altitude, mission.segments{i}.altitude], 'DisplayName', ['Cruise (', num2str(i-1), '-', num2str(i), ')']); + text(x + mission.segments{i}.range, mission.segments{i}.altitude, num2str(i)); x = x + mission.segments{i}.range; elseif strcmp(mission.segments{i}.type, 'hold') % Hold segment - plot([x, x + mission.segments{i}.range], [mission.segments{i}.altitude, mission.segments{i}.altitude], 'DisplayName', 'Hold'); + plot([x, x + mission.segments{i}.range], [mission.segments{i}.altitude, mission.segments{i}.altitude], 'DisplayName', ['Hold (', num2str(i-1), '-', num2str(i), ')']); + text(x + mission.segments{i}.range, mission.segments{i}.altitude, num2str(i)); x = x + mission.segments{i}.range; elseif strcmp(mission.segments{i}.type, 'descent') % Descent segment - plot([x, x + mission.segments{i}.range], [mission.segments{i}.altitude(1), mission.segments{i}.altitude(2)], 'DisplayName', 'Descent'); + plot([x, x + mission.segments{i}.range], [mission.segments{i}.altitude(1), mission.segments{i}.altitude(2)], 'DisplayName', ['Descent (', num2str(i-1), '-', num2str(i), ')']); + text(x + mission.segments{i}.range, mission.segments{i}.altitude(2), num2str(i)); x = x + mission.segments{i}.range; elseif strcmp(mission.segments{i}.type, 'vertical_descent') % Vertical descent segment - plot([x, x + mission.segments{i}.range], [mission.segments{i}.altitude(1), mission.segments{i}.altitude(2)], 'DisplayName', 'Vertical descent'); + plot([x, x + mission.segments{i}.range], [mission.segments{i}.altitude(1), mission.segments{i}.altitude(2)], 'DisplayName', ['Vertical Descent (', num2str(i-1), '-', num2str(i), ')']); + text(x + mission.segments{i}.range, mission.segments{i}.altitude(2), num2str(i)); x = x + mission.segments{i}.range; + elseif strcmp(mission.segments{i}.type, 'drop') % Drop segment + plot([x, x + 500], [mission.segments{i}.altitude, mission.segments{i}.altitude], 'DisplayName', ['Drop (', num2str(i-1), '-', num2str(i), ')']); + text([x + 250, x + 500], [mission.segments{i}.altitude - 10, mission.segments{i}.altitude], {'\downarrow', num2str(i)}); + x = x + 500; + elseif strcmp(mission.segments{i}.type, 'load') % Load segment + plot([x, x + 500], [mission.segments{i}.altitude, mission.segments{i}.altitude], 'DisplayName', ['Load (', num2str(i-1), '-', num2str(i), ')']); + text([x + 250, x + 500], [mission.segments{i}.altitude - 10, mission.segments{i}.altitude], {'\uparrow', num2str(i)}); + x = x + 500; elseif strcmp(mission.segments{i}.type, 'landing') % Landing segment - + plot([x, x + 500], [mission.segments{i}.altitude, mission.segments{i}.altitude], 'DisplayName', ['Landing (', num2str(i-1), '-', num2str(i), ')']); + text(x + 500, mission.segments{i}.altitude, num2str(i)); + x = x + 500; end end diff --git a/prettyjson.m b/prettyjson.m new file mode 100644 index 0000000..e0791a5 --- /dev/null +++ b/prettyjson.m @@ -0,0 +1,71 @@ +function [less_ugly] = prettyjson(ugly) +% Makes JSON strings (relatively) pretty +% Probably inefficient + +% Mostly meant for structures with simple strings and arrays; +% gets confused and !!mangles!! JSON when strings contain [ ] { or }. + + MAX_ARRAY_WIDTH = 80; + TAB = ' '; + + ugly = strrep(ugly, '{', sprintf('{\n')); + ugly = strrep(ugly, '}', sprintf('\n}')); + ugly = strrep(ugly, ',"', sprintf(', \n"')); + ugly = strrep(ugly, ',{', sprintf(', \n{')); + + indent = 0; + lines = splitlines(ugly); + + for i = 1:length(lines) + line = lines{i}; + next_indent = 0; + + % Count brackets + open_brackets = length(strfind(line, '[')); + close_brackets = length(strfind(line, ']')); + + open_braces = length(strfind(line, '{')); + close_braces = length(strfind(line, '}')); + + if close_brackets > open_brackets || close_braces > open_braces + indent = indent - 1; + end + + if open_brackets > close_brackets + line = strrep(line, '[', sprintf('[\n')); + next_indent = 1; + elseif open_brackets < close_brackets + line = strrep(line, ']', sprintf('\n]')); + next_indent = -1; + elseif open_brackets == close_brackets && length(line) > MAX_ARRAY_WIDTH + first_close_bracket = strfind(line, ']'); + if first_close_bracket > MAX_ARRAY_WIDTH % Just a long array -> each element on a new line + line = strrep(line, '[', sprintf('[\n%s', TAB)); + line = strrep(line, ']', sprintf('\n]')); + line = strrep(line, ',', sprintf(', \n%s', TAB)); % Add indents! + else % Nested array, probably 2d, first level is not too wide -> each sub-array on a new line + line = strrep(line, '[[', sprintf('[\n%s[', TAB)); + line = strrep(line, '],', sprintf('], \n%s', TAB)); % Add indents! + line = strrep(line, ']]', sprintf(']\n]')); + end + end + + sublines = splitlines(line); + for j = 1:length(sublines) + if j > 1 % todo: dumb to do this check at every line... + sublines{j} = sprintf('%s%s', repmat(TAB, 1, indent+next_indent), sublines{j}); + else + sublines{j} = sprintf('%s%s', repmat(TAB, 1, indent), sublines{j}); + end + end + + if open_brackets > close_brackets || open_braces > close_braces + indent = indent + 1; + end + indent = indent + next_indent; + lines{i} = strjoin(sublines, newline); + + end + + less_ugly = strjoin(lines, newline); +end \ No newline at end of file diff --git a/save_project.m b/save_project.m index ac38e00..fd58890 100644 --- a/save_project.m +++ b/save_project.m @@ -6,10 +6,7 @@ function save_project(data, filename) -str = jsonencode(data); -str = strrep(str, ',', sprintf(',\r')); -str = strrep(str, '[{', sprintf('[\r{\r')); -str = strrep(str, '}]', sprintf('\r}\r]')); +str = prettyjson(jsonencode(data)); fid = fopen(filename,'wt'); fprintf(fid, str);