From 6e79ec2831140dfe3d221d0101ce7a9ef4ad3184 Mon Sep 17 00:00:00 2001
From: Meric Taneri <meric.taneri@tum.de>
Date: Fri, 21 Feb 2025 09:18:22 +0100
Subject: [PATCH 1/4] plotting for mode_1 is improved
---
.../src/energy_based/energy_based_refactored.cpp | 12 ++++++++++++
1 file changed, 12 insertions(+)
diff --git a/constraint_analysis/src/energy_based/energy_based_refactored.cpp b/constraint_analysis/src/energy_based/energy_based_refactored.cpp
index d7e9e093..59233c69 100644
--- a/constraint_analysis/src/energy_based/energy_based_refactored.cpp
+++ b/constraint_analysis/src/energy_based/energy_based_refactored.cpp
@@ -140,7 +140,19 @@ void EnergyBased::operator()()
/* Plot the complete design chart with constraints, boundaries and design point*/
std::string plotting = this->configuration()->at("module_configuration_file/control_settings/plot_output/enable/value");
+
if (plotting == "true") {
+ /*Define relevant Wing_loading range. Selected values should be valid for all possible configurations*/
+ WS_start = 50;
+ WS_end = 820;
+ WS_size = WS_end - WS_start + 1; // Range from 50 to 820, inclusive
+ std::vector<double> W_over_S_data(WS_size);
+
+ /*Transform from kg/m^2 into N/m^2*/
+ for (int i = 0; i < WS_size; ++i) {
+ W_over_S_data[i] = (WS_start + i) * 9.81;
+ }
+ assess_constraints(W_over_S_data, constraint_analysis_tool, mission_data, engine);
ConstraintAnalysis_Plot ca_plot;
ca_plot.fill_infeasible_area(W_over_S_data, min_finder.dominant_thrust_to_weight, boundaries);
--
GitLab
From 621477c15889f1616da30b011c60dcd84b357146 Mon Sep 17 00:00:00 2001
From: Meric Taneri <meric.taneri@tum.de>
Date: Sun, 23 Feb 2025 22:21:58 +0100
Subject: [PATCH 2/4] plotting fix
---
.../energy_based/energy_based_refactored.cpp | 26 ++++++++++++-------
1 file changed, 16 insertions(+), 10 deletions(-)
diff --git a/constraint_analysis/src/energy_based/energy_based_refactored.cpp b/constraint_analysis/src/energy_based/energy_based_refactored.cpp
index 59233c69..4b4d16b6 100644
--- a/constraint_analysis/src/energy_based/energy_based_refactored.cpp
+++ b/constraint_analysis/src/energy_based/energy_based_refactored.cpp
@@ -105,7 +105,9 @@ void EnergyBased::operator()()
break;
}
}
-
+
+ Mattingly::constraint_analysis constraint_analysis_plotting = constraint_analysis_tool;
+
readMission mission_data(missionCSV);
assess_constraints(W_over_S_data, constraint_analysis_tool, mission_data, engine);
@@ -140,22 +142,26 @@ void EnergyBased::operator()()
/* Plot the complete design chart with constraints, boundaries and design point*/
std::string plotting = this->configuration()->at("module_configuration_file/control_settings/plot_output/enable/value");
-
if (plotting == "true") {
- /*Define relevant Wing_loading range. Selected values should be valid for all possible configurations*/
- WS_start = 50;
- WS_end = 820;
- WS_size = WS_end - WS_start + 1; // Range from 50 to 820, inclusive
- std::vector<double> W_over_S_data(WS_size);
- /*Transform from kg/m^2 into N/m^2*/
+ WS_start = 50;
for (int i = 0; i < WS_size; ++i) {
W_over_S_data[i] = (WS_start + i) * 9.81;
}
- assess_constraints(W_over_S_data, constraint_analysis_tool, mission_data, engine);
+
+ this->constraint_list = {};
+ this->boundaries = {};
+
+ assess_constraints(W_over_S_data, constraint_analysis_plotting, mission_data, engine);
+
+ Simple_Analysis min_finder_2 = Simple_Analysis(this->constraint_list, this->aircraft_xml(), this->boundaries);
+
+ min_finder_2.find_dominant_curve();
+
ConstraintAnalysis_Plot ca_plot;
- ca_plot.fill_infeasible_area(W_over_S_data, min_finder.dominant_thrust_to_weight, boundaries);
+ ca_plot.fill_infeasible_area(W_over_S_data, min_finder_2.dominant_thrust_to_weight, boundaries);
+
for (auto ca_case : this->constraint_list)
{
ca_plot.add_curve(ca_case.wing_loading, ca_case.thrust_to_weight, ca_case.constraint_case);
--
GitLab
From 6314d7953daa1b7d4b7b621bbb3b2f0e3a379e9a Mon Sep 17 00:00:00 2001
From: Meric Taneri <meric.taneri@tum.de>
Date: Mon, 24 Feb 2025 08:38:01 +0100
Subject: [PATCH 3/4] refactoring completed
---
constraint_analysis/CMakeLists.txt | 2 +-
.../src/energy_based/energy_based.cpp | 406 ++++++++++--------
.../energy_based/energy_based_refactored.cpp | 365 ----------------
3 files changed, 226 insertions(+), 547 deletions(-)
delete mode 100644 constraint_analysis/src/energy_based/energy_based_refactored.cpp
diff --git a/constraint_analysis/CMakeLists.txt b/constraint_analysis/CMakeLists.txt
index 20918bf4..564eb7d2 100644
--- a/constraint_analysis/CMakeLists.txt
+++ b/constraint_analysis/CMakeLists.txt
@@ -7,7 +7,7 @@ set(MODULE_SOURCES
src/constraint_analysis.cpp
src/ca_parser.cpp
src/io/aircraft_xml.cpp
- src/energy_based/energy_based_refactored.cpp
+ src/energy_based/energy_based.cpp
)
add_executable(${MODULE_NAME}
diff --git a/constraint_analysis/src/energy_based/energy_based.cpp b/constraint_analysis/src/energy_based/energy_based.cpp
index 7cc47455..4b4d16b6 100644
--- a/constraint_analysis/src/energy_based/energy_based.cpp
+++ b/constraint_analysis/src/energy_based/energy_based.cpp
@@ -62,6 +62,7 @@ void EnergyBased::operator()()
int WS_start = 0;
int WS_end = 0;
size_t WS_size = 1;
+
if (mode_selection == "mode_0")
{
/*Define relevant Wing_loading range. Selected values should be valid for all possible configurations*/
@@ -71,8 +72,8 @@ void EnergyBased::operator()()
}
else if (mode_selection == "mode_1")
{
- WS_start = this->aircraft_xml()->at("aircraft_exchange_file/sizing_point/wing_loading/value");
- WS_size = 1; // Constant W/S
+ WS_start = this->aircraft_xml()->at("aircraft_exchange_file/sizing_point/wing_loading/value");
+ WS_size = 1; // Constant W/S
}
std::vector<double> W_over_S_data(WS_size);
@@ -82,6 +83,8 @@ void EnergyBased::operator()()
W_over_S_data[i] = (WS_start + i) * 9.81;
}
+ Engine engine = Engine(this->engine_directory());
+
Mattingly::constraint_analysis constraint_analysis_tool;
atmosphere atm;
@@ -92,31 +95,8 @@ void EnergyBased::operator()()
constraint_analysis_tool.atm = atm;
- double CD = 0.0;
- double CL = 0.0;
- std::vector<double> CD_vector = {};
- std::vector<double> CL_vector = {};
-
- Engine engine = Engine(this->engine_directory());
- auto calculate_thrust_lapse = [&atm, &engine](double M, double altitude) {
- const double thrust_lapse = engine.get_thrust_lapse("cruise", atm, M, altitude);
- return thrust_lapse;
- };
-
- auto get_num_engines = [&]() {
- double num_engines = 0;
-
- // Get the node vector for all defined propulsion systems
- auto propulsors = this->aircraft_xml()->getVector(
- "aircraft_exchange_file/requirements_and_specifications/design_specification/propulsion/propulsor");
-
- // Increment num_engines for each propulsor
- for (const auto& propulsor : propulsors) {
- num_engines += 1;
- }
- return num_engines;
- };
std::filesystem::path missionCSV = this->mission_directory();
+
for (const auto& entry : fs::directory_iterator(this->mission_directory())) {
std::string filename = entry.path().filename().string();
@@ -125,157 +105,15 @@ void EnergyBased::operator()()
break;
}
}
- readMission mission_data(missionCSV);
-
- /* Take-off Segment */
- double CL_max_TO = this->aircraft_xml()->at("aircraft_exchange_file/analysis/aerodynamics/lift_coefficients/C_LmaxT-O/value");
- double s_G_TO = this->aircraft_xml()->at("aircraft_exchange_file/requirements_and_specifications/requirements/top_level_aircraft_requirements/takeoff_distance/value");
- double takeoff_climb_angle = this->configuration()->at("module_configuration_file/program_settings/constraint_selection/standard_set/takeoff_climb/takeoff_climb_angle/value");
- double M_TO = this->configuration()->at("module_configuration_file/program_settings/constraint_selection/standard_set/takeoff_ground_roll/Mach_TO/value");
- double thrust_lapse_TO = calculate_thrust_lapse(M_TO, 0.0);
- double weight_fraction_TO = mission_data.get_beta(" takeoff", 1.0); //from mission analysis
- const auto to_ground_roll = constraint_analysis_tool.takeoff_ground_roll(
- W_over_S_data,
- weight_fraction_TO,
- thrust_lapse_TO,
- M_TO,
- 0.0,
- s_G_TO,
- CL_max_TO);
-
- readPolar takeoff_config(this->polar(), "takeoff", "linear", M_TO);
- std::vector<double> coeffs = takeoff_config.quadratic_coefficients();
- double K_1 = coeffs[0];
- double K_2 = coeffs[1];
- double C_D0 = coeffs[2];
-
- for (int i = 0; i < W_over_S_data.size(); ++i) {
- CD = takeoff_config.getCD(W_over_S_data[i], 1.0, weight_fraction_TO, M_TO, 0.0);
- CD_vector.push_back(CD);
- };
- const auto to_climb_angle = constraint_analysis_tool.takeoff_climb_angle(
- W_over_S_data,
- CD_vector,
- weight_fraction_TO,
- thrust_lapse_TO,
- M_TO,
- 0.0,
- CL_max_TO,
- takeoff_climb_angle);
-
- /* OEI Climb */
- double climb_speed = this->aircraft_xml()->at(
- "aircraft_exchange_file/requirements_and_specifications/requirements/top_level_aircraft_requirements/design_mission/climb_speed_schedule/climb_speed_below_FL100/value");
- double climb_gradient_OEI = this->configuration()->at("module_configuration_file/program_settings/constraint_selection/standard_set/OEI_climb/climb_gradient_OEI/value");
- double M_climb = climb_speed / atm.getSpeedOfSound(0.0);
-
- double num_engines = get_num_engines();
-
- readPolar climb_config(this->polar(), "takeoff", "linear", M_TO);
- for (int i = 0; i < W_over_S_data.size(); ++i) {
- CD = climb_config.getCD(W_over_S_data[i], 1.0, weight_fraction_TO, (climb_speed / atm.getSpeedOfSound(0.0)), 0.0);
- CD_vector.push_back(CD);
- };
+ Mattingly::constraint_analysis constraint_analysis_plotting = constraint_analysis_tool;
- double thrust_factor_OEI = thrust_lapse_TO *((num_engines - 1)/ num_engines);
- double climb_rate_OEI = climb_gradient_OEI * climb_speed;
- double weight_fraction_climb_OEI = mission_data.get_beta(" takeoff");
- /*minimal climb rate of 2.4% defined in CS25 for second segment climb*/
- const auto climb_OEI = constraint_analysis_tool.constant_speed_climb(
- W_over_S_data,
- CD_vector,
- weight_fraction_climb_OEI,
- thrust_factor_OEI,
- M_climb,
- 0.0,
- climb_rate_OEI); /* DEFINE CLIMB RATE CONSTRAINT */
-
- /* Service Ceiling */
- double altitude_max = this->aircraft_xml()->at(
- "aircraft_exchange_file/requirements_and_specifications/requirements/top_level_aircraft_requirements/flight_envelope/maximum_operating_altitude/value");
- double M_cruise = this->aircraft_xml()->at(
- "aircraft_exchange_file/requirements_and_specifications/requirements/top_level_aircraft_requirements/study_mission/initial_cruise_mach_number/value");
-
- double M_polar = M_cruise;
- double weight_fraction_climb_service_ceiling = mission_data.get_beta(" climb_to_cruise");
-
- bool polar_read = false;
- readPolar clean_config_M_max(this->polar(), "clean", "linear", M_TO);
- while (!polar_read) {
- try {
- clean_config_M_max = readPolar(this->polar(), "clean", "linear", M_polar);
- if (clean_config_M_max.CD_polar.size() >= 1)
- {
- polar_read = true;
- }
- else {
- M_polar += 0.01;
- };
- }
- catch (...) {
- M_polar += 0.01;
- }
- }
- double thrust_lapse_segment = calculate_thrust_lapse(M_cruise, altitude_max);
- CD_vector = {};
- CL_vector = {};
- coeffs = clean_config_M_max.quadratic_coefficients();
- K_1 = coeffs[0];
- K_2 = coeffs[1];
- C_D0 = coeffs[2];
-
- for (int i = 0; i < W_over_S_data.size(); ++i) {
- CD = clean_config_M_max.getCD(W_over_S_data[i], 1.0, weight_fraction_climb_service_ceiling, M_polar, altitude_max);
- CD_vector.push_back(CD);
- };
+ readMission mission_data(missionCSV);
- double minimum_climb_rate = this->configuration()->at("module_configuration_file/program_settings/constraint_selection/standard_set/SEP/minimum_climb_rate/value");
-
- /* Service ceiling with the minimum climb rate requirement */
- const auto service_ceiling = constraint_analysis_tool.constant_speed_climb(
- W_over_S_data,
- CD_vector,
- weight_fraction_climb_service_ceiling,
- thrust_lapse_segment,
- M_cruise,
- altitude_max,
- minimum_climb_rate);
-
- /* Landing */
- double thrust_lapse_landing = calculate_thrust_lapse(M_TO, 0.0);
- double CL_max_L = this->aircraft_xml()->at("aircraft_exchange_file/analysis/aerodynamics/lift_coefficients/C_LmaxLanding/value");
- double my_B = this->aircraft_xml()->at("aircraft_exchange_file/requirements_and_specifications/requirements/additional_requirements/landing_gear/braking_coefficient/value");
- double s_G_L = this->aircraft_xml()->at("aircraft_exchange_file/requirements_and_specifications/requirements/top_level_aircraft_requirements/takeoff_distance/value");
- readPolar landing_config(this->polar(), "landing", "linear", M_TO);
- double CD_max_L = landing_config.interpCD(CL_max_L);
- double max_landing_weight = this->aircraft_xml()->at("aircraft_exchange_file/analysis/masses_cg_inertia/maximum_landing_mass/mass_properties/mass/value");
- double MTOM = this->aircraft_xml()->at("aircraft_exchange_file/analysis/masses_cg_inertia/maximum_takeoff_mass/mass_properties/mass/value");
- double weight_fraction_landing = max_landing_weight / MTOM;
- const auto landing = constraint_analysis_tool.braking_roll(CD_max_L, CL_max_L, weight_fraction_landing, thrust_lapse_landing, M_TO, 0.0, my_B, s_G_L);
-
- /* Gust */
- double AoA_1 = landing_config.interpalpha(CL_max_L);
- double AoA_2 = landing_config.interpalpha(CL_max_L - CL_max_L * 0.9);
- double CL_alpha = (CL_max_L - (CL_max_L - CL_max_L * 0.9)) / ((AoA_1 - AoA_2) * PI / 180.0);
- double V_TO_L = M_TO * atm.getSpeedOfSound(0.0);
- double w_g = this->configuration()->at("module_configuration_file/program_settings/constraint_selection/standard_set/gust/gust_speed/value");
- double dn_G = this->configuration()->at("module_configuration_file/program_settings/constraint_selection/standard_set/gust/gust_load_factor/value");
- const auto gust = constraint_analysis_tool.gust(CL_alpha, 0.0, V_TO_L, w_g, dn_G, weight_fraction_TO);
-
- /*Define Constraint out of each non-vertical constraint case*/
- Constraint Case1(W_over_S_data, to_ground_roll, "Takeoff Ground Roll");
- Constraint Case2(W_over_S_data, to_climb_angle, "Takeoff Climb Angle");
- Constraint Case3(W_over_S_data, climb_OEI, "OEI Climb");
- Constraint Case4(W_over_S_data, service_ceiling, "Service Ceiling");
-
- std::vector<Constraint> constraint_list = { Case1, Case2, Case3, Case4 };
-
- /* Define boundaries out of vertical constraint cases*/
- std::vector<double> boundaries = { gust, landing };
+ assess_constraints(W_over_S_data, constraint_analysis_tool, mission_data, engine);
/* Define Simple_Analysis*/
- Simple_Analysis min_finder = Simple_Analysis(constraint_list, this->aircraft_xml(), boundaries);
+ Simple_Analysis min_finder = Simple_Analysis(this->constraint_list, this->aircraft_xml(), this->boundaries);
/* Find dominant curve of the design chart*/
min_finder.find_dominant_curve();
@@ -294,28 +132,234 @@ void EnergyBased::operator()()
/*Update design point in ACxml*/
min_finder.update_design_point();
- std::string plot_dir_string=this->plot_dir().string();
+ std::string plot_dir_string = this->plot_dir().string();
// Add output of design point in logfile
- myRuntimeInfo->out << "New wing loading: " << min_finder.design_wing_loading * 9.81 << " N/m^2. New thrust to weight ratio: " << min_finder.design_thrust_to_weight << std::endl;
+ myRuntimeInfo->out << "New wing loading: " << min_finder.design_wing_loading * 9.81 << " N/m^2. New thrust to weight ratio: " << min_finder.design_thrust_to_weight << std::endl;
min_finder.write_design_point(plot_dir_string + "/csv_files/design_point.csv");
/* Plot the complete design chart with constraints, boundaries and design point*/
std::string plotting = this->configuration()->at("module_configuration_file/control_settings/plot_output/enable/value");
if (plotting == "true") {
+
+ WS_start = 50;
+ for (int i = 0; i < WS_size; ++i) {
+ W_over_S_data[i] = (WS_start + i) * 9.81;
+ }
+
+ this->constraint_list = {};
+ this->boundaries = {};
+
+ assess_constraints(W_over_S_data, constraint_analysis_plotting, mission_data, engine);
+
+ Simple_Analysis min_finder_2 = Simple_Analysis(this->constraint_list, this->aircraft_xml(), this->boundaries);
+
+ min_finder_2.find_dominant_curve();
+
ConstraintAnalysis_Plot ca_plot;
- ca_plot.fill_infeasible_area(W_over_S_data, min_finder.dominant_thrust_to_weight, boundaries);
- ca_plot.add_curve(W_over_S_data, to_ground_roll, "Takeoff Ground Roll");
- ca_plot.add_curve(W_over_S_data, to_climb_angle, "Takeoff Climb Angle");
- ca_plot.add_curve(W_over_S_data, climb_OEI, "OEI Climb");
- ca_plot.add_curve(W_over_S_data, service_ceiling, "Service Ceiling");
- ca_plot.add_curve(landing, "Landing");
- ca_plot.add_curve(gust, "Gust");
+ ca_plot.fill_infeasible_area(W_over_S_data, min_finder_2.dominant_thrust_to_weight, boundaries);
+
+ for (auto ca_case : this->constraint_list)
+ {
+ ca_plot.add_curve(ca_case.wing_loading, ca_case.thrust_to_weight, ca_case.constraint_case);
+ }
+ for (auto boundary : this->boundaries)
+ {
+ ca_plot.add_curve(boundary, "W/S limit");
+ }
ca_plot.add_curve({ min_finder.design_wing_loading * 9.81 }, { min_finder.design_thrust_to_weight }, "Design Point");
ca_plot.save_figure(plot_dir_string + "/constraint_plot.svg");
};
}
+
+void EnergyBased::assess_constraints(std::vector<double>& W_over_S_data, Mattingly::constraint_analysis& constraint_analysis_tool, readMission& mission_data, Engine& engine)
+{
+ auto calculate_thrust_lapse = [&constraint_analysis_tool, &engine](double M, double altitude) {
+ const double thrust_lapse = engine.get_thrust_lapse("cruise", constraint_analysis_tool.atm, M, altitude);
+ return thrust_lapse;
+ };
+
+ auto get_num_engines = [&]() {
+ double num_engines = 0;
+
+ // Get the node vector for all defined propulsion systems
+ auto propulsors = this->aircraft_xml()->getVector(
+ "aircraft_exchange_file/requirements_and_specifications/design_specification/propulsion/propulsor");
+
+ // Increment num_engines for each propulsor
+ for (const auto& propulsor : propulsors) {
+ num_engines += 1;
+ }
+ return num_engines;
+ };
+
+ auto std_constraints = this->configuration()->at(
+ "module_configuration_file/program_settings/constraint_selection/standard_set").getChildren();
+
+ for (auto constraint : std_constraints)
+ {
+ std::string activated = constraint->at("active/value");
+ if (activated == "true")
+ {
+ double CD = 0.0;
+ double CL = 0.0;
+ std::vector<double> CD_vector = {};
+ std::vector<double> CL_vector = {};
+
+ if (constraint->name == "takeoff_ground_roll")
+ {
+ double CL_max_TO = this->aircraft_xml()->at("aircraft_exchange_file/analysis/aerodynamics/lift_coefficients/C_LmaxT-O/value");
+ double s_G_TO = this->aircraft_xml()->at("aircraft_exchange_file/requirements_and_specifications/requirements/top_level_aircraft_requirements/takeoff_distance/value");
+ double M = constraint->at("Mach_TO/value");
+ double thrust_lapse = calculate_thrust_lapse(M, 0.0);
+ double weight_fraction = mission_data.get_beta(" takeoff", 1.0); //from mission analysis
+ const auto constraint_vector = constraint_analysis_tool.takeoff_ground_roll(
+ W_over_S_data,
+ weight_fraction,
+ thrust_lapse,
+ M,
+ 0.0,
+ s_G_TO,
+ CL_max_TO);
+ Constraint Case(W_over_S_data, constraint_vector, constraint->name);
+ this->constraint_list.push_back(Case);
+ }
+ else if (constraint->name == "takeoff_climb")
+ {
+ double M = constraint->at("Mach_TO/value");
+ double takeoff_climb_angle = constraint->at("takeoff_climb_angle/value");
+ double thrust_lapse = calculate_thrust_lapse(M, 0.0);
+ double weight_fraction = mission_data.get_beta(" takeoff", 1.0); //from mission analysis
+ double CL_max_TO = this->aircraft_xml()->at("aircraft_exchange_file/analysis/aerodynamics/lift_coefficients/C_LmaxT-O/value");
+ readPolar takeoff_config(this->polar(), "takeoff", "linear", M);
+ CD = takeoff_config.interpCD(CL_max_TO);
+ for (int i = 0; i < W_over_S_data.size(); ++i) {
+ CD_vector.push_back(CD);
+ };
+ const auto constraint_vector = constraint_analysis_tool.takeoff_climb_angle(
+ W_over_S_data,
+ CD_vector,
+ weight_fraction,
+ thrust_lapse,
+ M,
+ 0.0,
+ CL_max_TO,
+ takeoff_climb_angle);
+ Constraint Case(W_over_S_data, constraint_vector, constraint->name);
+ this->constraint_list.push_back(Case);
+ }
+ else if (constraint->name == "OEI_climb")
+ {
+ /* OEI Climb */
+ double climb_speed = this->aircraft_xml()->at(
+ "aircraft_exchange_file/requirements_and_specifications/requirements/top_level_aircraft_requirements/design_mission/climb_speed_schedule/climb_speed_below_FL100/value");
+ double climb_gradient_OEI = constraint->at("climb_gradient_OEI/value");
+ double M = climb_speed / constraint_analysis_tool.atm.getSpeedOfSound(0.0);
+ double weight_fraction = mission_data.get_beta(" takeoff"); //from mission analysis
+ double thrust_lapse = calculate_thrust_lapse(M, 0.0);
+ double num_engines = get_num_engines();
+
+ readPolar climb_config(this->polar(), "takeoff", "linear", M);
+ for (int i = 0; i < W_over_S_data.size(); ++i) {
+ CD = climb_config.getCD(W_over_S_data[i], 1.0, weight_fraction, (climb_speed / constraint_analysis_tool.atm.getSpeedOfSound(0.0)), 0.0);
+ CD_vector.push_back(CD);
+ };
+
+ double thrust_factor_OEI = thrust_lapse * ((num_engines - 1) / num_engines);
+ double climb_rate_OEI = climb_gradient_OEI * climb_speed;
+ /*minimal climb rate of 2.4% defined in CS25 for second segment climb*/
+ const auto constraint_vector = constraint_analysis_tool.constant_speed_climb(
+ W_over_S_data,
+ CD_vector,
+ weight_fraction,
+ thrust_factor_OEI,
+ M,
+ 0.0,
+ climb_rate_OEI); /* DEFINE CLIMB RATE CONSTRAINT */
+ Constraint Case(W_over_S_data, constraint_vector, constraint->name);
+ this->constraint_list.push_back(Case);
+ }
+ else if (constraint->name == "top_of_climb")
+ {
+
+ }
+ else if (constraint->name == "cruise")
+ {
+
+ }
+ else if (constraint->name == "SEP")
+ {
+ /* Service Ceiling */
+ double altitude = this->aircraft_xml()->at(
+ "aircraft_exchange_file/requirements_and_specifications/requirements/top_level_aircraft_requirements/flight_envelope/maximum_operating_altitude/value");
+ double M = this->aircraft_xml()->at(
+ "aircraft_exchange_file/requirements_and_specifications/requirements/top_level_aircraft_requirements/study_mission/initial_cruise_mach_number/value");
+
+ double weight_fraction = mission_data.get_beta(" climb_to_cruise");
+ readPolar clean_config_M_cruise(this->polar(), "clean", "linear", M);
+ double thrust_lapse = calculate_thrust_lapse(M, altitude);
+ CD_vector = {};
+ CL_vector = {};
+
+ for (int i = 0; i < W_over_S_data.size(); ++i) {
+ CD = clean_config_M_cruise.getCD(W_over_S_data[i], 1.0, weight_fraction, M, altitude);
+ CD_vector.push_back(CD);
+ }
+
+ double minimum_climb_rate = constraint->at("minimum_climb_rate/value");
+
+ /* Service ceiling with the minimum climb rate requirement */
+ const auto constraint_vector = constraint_analysis_tool.constant_speed_climb(
+ W_over_S_data,
+ CD_vector,
+ weight_fraction,
+ thrust_lapse,
+ M,
+ altitude,
+ minimum_climb_rate);
+ Constraint Case(W_over_S_data, constraint_vector, constraint->name);
+ this->constraint_list.push_back(Case);
+ }
+ else if (constraint->name == "approach_speed")
+ {
+
+ }
+ else if (constraint->name == "landing_field_length")
+ {
+ /* Landing */
+ double M = constraint->at("Mach/value");
+ double thrust_lapse = calculate_thrust_lapse(M, 0.0);
+ double CL_max_L = this->aircraft_xml()->at("aircraft_exchange_file/analysis/aerodynamics/lift_coefficients/C_LmaxLanding/value");
+ double my_B = this->aircraft_xml()->at("aircraft_exchange_file/requirements_and_specifications/requirements/additional_requirements/landing_gear/braking_coefficient/value");
+ double s_G_L = this->aircraft_xml()->at("aircraft_exchange_file/requirements_and_specifications/requirements/top_level_aircraft_requirements/takeoff_distance/value");
+ readPolar landing_config(this->polar(), "landing", "linear", M);
+ double CD_max_L = landing_config.interpCD(CL_max_L);
+ double max_landing_weight = this->aircraft_xml()->at("aircraft_exchange_file/analysis/masses_cg_inertia/maximum_landing_mass/mass_properties/mass/value");
+ double MTOM = this->aircraft_xml()->at("aircraft_exchange_file/analysis/masses_cg_inertia/maximum_takeoff_mass/mass_properties/mass/value");
+ double weight_fraction = max_landing_weight / MTOM;
+ const auto constraint_border = constraint_analysis_tool.braking_roll(CD_max_L, CL_max_L, weight_fraction, thrust_lapse, M, 0.0, my_B, s_G_L);
+ this->boundaries.push_back(constraint_border);
+ }
+ else if (constraint->name == "gust")
+ {
+ /* Gust */
+ double M = constraint->at("Mach/value");
+ readPolar landing_config(this->polar(), "landing", "linear", M);
+ double weight_fraction = mission_data.get_beta(" takeoff", 1.0); //from mission analysis
+ double CL_max_L = this->aircraft_xml()->at("aircraft_exchange_file/analysis/aerodynamics/lift_coefficients/C_LmaxLanding/value");
+ double AoA_1 = landing_config.interpalpha(CL_max_L);
+ double AoA_2 = landing_config.interpalpha(CL_max_L - CL_max_L * 0.9);
+ double CL_alpha = (CL_max_L - (CL_max_L - CL_max_L * 0.9)) / ((AoA_1 - AoA_2) * PI / 180.0);
+ double V_TO_L = this->aircraft_xml()->at("aircraft_exchange_file/requirements_and_specifications/requirements/top_level_aircraft_requirements/flight_envelope/maximum_approach_speed/value");
+ double w_g = this->configuration()->at("module_configuration_file/program_settings/constraint_selection/standard_set/gust/gust_speed/value");
+ double dn_G = this->configuration()->at("module_configuration_file/program_settings/constraint_selection/standard_set/gust/gust_load_factor/value");
+ const auto constraint_border = constraint_analysis_tool.gust(CL_alpha, 0.0, V_TO_L, w_g, dn_G, weight_fraction);
+ this->boundaries.push_back(constraint_border);
+ }
+ }
+ };
+};
\ No newline at end of file
diff --git a/constraint_analysis/src/energy_based/energy_based_refactored.cpp b/constraint_analysis/src/energy_based/energy_based_refactored.cpp
deleted file mode 100644
index 4b4d16b6..00000000
--- a/constraint_analysis/src/energy_based/energy_based_refactored.cpp
+++ /dev/null
@@ -1,365 +0,0 @@
-/*
- * UNICADO - UNIversity Conceptual Aircraft Design and Optimization
- *
- * Copyright (C) 2024 UNICADO consortium
- *
- * 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 <https://www.gnu.org/licenses/>.
- *
- * Description:
- * This file is part of UNICADO.
- */
-
-#include "energy_based.h"
-#include <algorithm>
-#include <filesystem>
-#include <format>
-#include <list>
-#include <ranges>
-#include <cmath>
-#include <vector>
-#include <numbers>
-#include <iostream>
-#include <string>
-#include <moduleBasics/module.h>
-#include <engine/engine.h>
-#include "../../include/constraint_analysis/ca_functions.h"
-#include "../../include/constraint_analysis/ca_parser.h"
-#include "../constraint_analysis.h"
-#include "../../include/constraint_analysis/ca_plotting.h"
-#include "../io/aircraft_xml.h"
-#include "../../include/constraint_analysis/ca_minfinder.h"
-
-
-void EnergyBased::initialize()
-{
- if (!this->aircraft_xml()) {
- throw std::runtime_error("[EnergyBased::initialize()] aircraft_xml is null!");
- };
- if (!this->configuration()) {
- throw std::runtime_error("[EnergyBased::initialize()] config_xml is null!");
- };
- if (!this->polar()) {
- throw std::runtime_error("[EnergyBased::initialize()] polar_xml is null!");
- };
-}
-
-void EnergyBased::operator()()
-{
-
- const std::string mode_selection(this->configuration()->at("module_configuration_file/program_settings/mode_selection/value"));
-
- int WS_start = 0;
- int WS_end = 0;
- size_t WS_size = 1;
-
- if (mode_selection == "mode_0")
- {
- /*Define relevant Wing_loading range. Selected values should be valid for all possible configurations*/
- WS_start = 50;
- WS_end = 820;
- WS_size = WS_end - WS_start + 1; // Range from 50 to 820, inclusive
- }
- else if (mode_selection == "mode_1")
- {
- WS_start = this->aircraft_xml()->at("aircraft_exchange_file/sizing_point/wing_loading/value");
- WS_size = 1; // Constant W/S
- }
-
- std::vector<double> W_over_S_data(WS_size);
-
- /*Transform from kg/m^2 into N/m^2*/
- for (int i = 0; i < WS_size; ++i) {
- W_over_S_data[i] = (WS_start + i) * 9.81;
- }
-
- Engine engine = Engine(this->engine_directory());
-
- Mattingly::constraint_analysis constraint_analysis_tool;
-
- atmosphere atm;
-
- // Set atmosphere according to the delta ISA requirement for the sizing process
- double temperature_ISA = ISA_TEMPERATURE + this->aircraft_xml()->at("aircraft_exchange_file/requirements_and_specifications/requirements/top_level_aircraft_requirements/design_mission/delta_ISA/value");
- atm.setAtmosphere(0.0, temperature_ISA, ISA_PRESSURE);
-
- constraint_analysis_tool.atm = atm;
-
- std::filesystem::path missionCSV = this->mission_directory();
-
- for (const auto& entry : fs::directory_iterator(this->mission_directory())) {
- std::string filename = entry.path().filename().string();
-
- if (filename.find("design_mission") != std::string::npos && filename.find("_out.csv") != std::string::npos) {
- missionCSV /= filename;
- break;
- }
- }
-
- Mattingly::constraint_analysis constraint_analysis_plotting = constraint_analysis_tool;
-
- readMission mission_data(missionCSV);
-
- assess_constraints(W_over_S_data, constraint_analysis_tool, mission_data, engine);
-
- /* Define Simple_Analysis*/
- Simple_Analysis min_finder = Simple_Analysis(this->constraint_list, this->aircraft_xml(), this->boundaries);
-
- /* Find dominant curve of the design chart*/
- min_finder.find_dominant_curve();
-
- double safety_factor = this->configuration()->at("module_configuration_file/program_settings/safety_factor/value");
-
- const double buffer_factor(this->configuration()->at("module_configuration_file/program_settings/buffer_factor/value"));
-
- /*Get the previous design point*/
- double previous_WS = this->aircraft_xml()->at("aircraft_exchange_file/sizing_point/wing_loading/value");
- double previous_TW = this->aircraft_xml()->at("aircraft_exchange_file/sizing_point/thrust_to_weight/value");
-
- /*Find design point with safety factor*/
- min_finder.find_converged_design_point(previous_TW, previous_WS * 9.81, safety_factor, buffer_factor);
-
- /*Update design point in ACxml*/
- min_finder.update_design_point();
-
- std::string plot_dir_string = this->plot_dir().string();
-
- // Add output of design point in logfile
-
- myRuntimeInfo->out << "New wing loading: " << min_finder.design_wing_loading * 9.81 << " N/m^2. New thrust to weight ratio: " << min_finder.design_thrust_to_weight << std::endl;
-
- min_finder.write_design_point(plot_dir_string + "/csv_files/design_point.csv");
-
- /* Plot the complete design chart with constraints, boundaries and design point*/
- std::string plotting = this->configuration()->at("module_configuration_file/control_settings/plot_output/enable/value");
- if (plotting == "true") {
-
- WS_start = 50;
- for (int i = 0; i < WS_size; ++i) {
- W_over_S_data[i] = (WS_start + i) * 9.81;
- }
-
- this->constraint_list = {};
- this->boundaries = {};
-
- assess_constraints(W_over_S_data, constraint_analysis_plotting, mission_data, engine);
-
- Simple_Analysis min_finder_2 = Simple_Analysis(this->constraint_list, this->aircraft_xml(), this->boundaries);
-
- min_finder_2.find_dominant_curve();
-
- ConstraintAnalysis_Plot ca_plot;
-
- ca_plot.fill_infeasible_area(W_over_S_data, min_finder_2.dominant_thrust_to_weight, boundaries);
-
- for (auto ca_case : this->constraint_list)
- {
- ca_plot.add_curve(ca_case.wing_loading, ca_case.thrust_to_weight, ca_case.constraint_case);
- }
- for (auto boundary : this->boundaries)
- {
- ca_plot.add_curve(boundary, "W/S limit");
- }
- ca_plot.add_curve({ min_finder.design_wing_loading * 9.81 }, { min_finder.design_thrust_to_weight }, "Design Point");
-
- ca_plot.save_figure(plot_dir_string + "/constraint_plot.svg");
- };
-}
-
-void EnergyBased::assess_constraints(std::vector<double>& W_over_S_data, Mattingly::constraint_analysis& constraint_analysis_tool, readMission& mission_data, Engine& engine)
-{
- auto calculate_thrust_lapse = [&constraint_analysis_tool, &engine](double M, double altitude) {
- const double thrust_lapse = engine.get_thrust_lapse("cruise", constraint_analysis_tool.atm, M, altitude);
- return thrust_lapse;
- };
-
- auto get_num_engines = [&]() {
- double num_engines = 0;
-
- // Get the node vector for all defined propulsion systems
- auto propulsors = this->aircraft_xml()->getVector(
- "aircraft_exchange_file/requirements_and_specifications/design_specification/propulsion/propulsor");
-
- // Increment num_engines for each propulsor
- for (const auto& propulsor : propulsors) {
- num_engines += 1;
- }
- return num_engines;
- };
-
- auto std_constraints = this->configuration()->at(
- "module_configuration_file/program_settings/constraint_selection/standard_set").getChildren();
-
- for (auto constraint : std_constraints)
- {
- std::string activated = constraint->at("active/value");
- if (activated == "true")
- {
- double CD = 0.0;
- double CL = 0.0;
- std::vector<double> CD_vector = {};
- std::vector<double> CL_vector = {};
-
- if (constraint->name == "takeoff_ground_roll")
- {
- double CL_max_TO = this->aircraft_xml()->at("aircraft_exchange_file/analysis/aerodynamics/lift_coefficients/C_LmaxT-O/value");
- double s_G_TO = this->aircraft_xml()->at("aircraft_exchange_file/requirements_and_specifications/requirements/top_level_aircraft_requirements/takeoff_distance/value");
- double M = constraint->at("Mach_TO/value");
- double thrust_lapse = calculate_thrust_lapse(M, 0.0);
- double weight_fraction = mission_data.get_beta(" takeoff", 1.0); //from mission analysis
- const auto constraint_vector = constraint_analysis_tool.takeoff_ground_roll(
- W_over_S_data,
- weight_fraction,
- thrust_lapse,
- M,
- 0.0,
- s_G_TO,
- CL_max_TO);
- Constraint Case(W_over_S_data, constraint_vector, constraint->name);
- this->constraint_list.push_back(Case);
- }
- else if (constraint->name == "takeoff_climb")
- {
- double M = constraint->at("Mach_TO/value");
- double takeoff_climb_angle = constraint->at("takeoff_climb_angle/value");
- double thrust_lapse = calculate_thrust_lapse(M, 0.0);
- double weight_fraction = mission_data.get_beta(" takeoff", 1.0); //from mission analysis
- double CL_max_TO = this->aircraft_xml()->at("aircraft_exchange_file/analysis/aerodynamics/lift_coefficients/C_LmaxT-O/value");
- readPolar takeoff_config(this->polar(), "takeoff", "linear", M);
- CD = takeoff_config.interpCD(CL_max_TO);
- for (int i = 0; i < W_over_S_data.size(); ++i) {
- CD_vector.push_back(CD);
- };
- const auto constraint_vector = constraint_analysis_tool.takeoff_climb_angle(
- W_over_S_data,
- CD_vector,
- weight_fraction,
- thrust_lapse,
- M,
- 0.0,
- CL_max_TO,
- takeoff_climb_angle);
- Constraint Case(W_over_S_data, constraint_vector, constraint->name);
- this->constraint_list.push_back(Case);
- }
- else if (constraint->name == "OEI_climb")
- {
- /* OEI Climb */
- double climb_speed = this->aircraft_xml()->at(
- "aircraft_exchange_file/requirements_and_specifications/requirements/top_level_aircraft_requirements/design_mission/climb_speed_schedule/climb_speed_below_FL100/value");
- double climb_gradient_OEI = constraint->at("climb_gradient_OEI/value");
- double M = climb_speed / constraint_analysis_tool.atm.getSpeedOfSound(0.0);
- double weight_fraction = mission_data.get_beta(" takeoff"); //from mission analysis
- double thrust_lapse = calculate_thrust_lapse(M, 0.0);
- double num_engines = get_num_engines();
-
- readPolar climb_config(this->polar(), "takeoff", "linear", M);
- for (int i = 0; i < W_over_S_data.size(); ++i) {
- CD = climb_config.getCD(W_over_S_data[i], 1.0, weight_fraction, (climb_speed / constraint_analysis_tool.atm.getSpeedOfSound(0.0)), 0.0);
- CD_vector.push_back(CD);
- };
-
- double thrust_factor_OEI = thrust_lapse * ((num_engines - 1) / num_engines);
- double climb_rate_OEI = climb_gradient_OEI * climb_speed;
- /*minimal climb rate of 2.4% defined in CS25 for second segment climb*/
- const auto constraint_vector = constraint_analysis_tool.constant_speed_climb(
- W_over_S_data,
- CD_vector,
- weight_fraction,
- thrust_factor_OEI,
- M,
- 0.0,
- climb_rate_OEI); /* DEFINE CLIMB RATE CONSTRAINT */
- Constraint Case(W_over_S_data, constraint_vector, constraint->name);
- this->constraint_list.push_back(Case);
- }
- else if (constraint->name == "top_of_climb")
- {
-
- }
- else if (constraint->name == "cruise")
- {
-
- }
- else if (constraint->name == "SEP")
- {
- /* Service Ceiling */
- double altitude = this->aircraft_xml()->at(
- "aircraft_exchange_file/requirements_and_specifications/requirements/top_level_aircraft_requirements/flight_envelope/maximum_operating_altitude/value");
- double M = this->aircraft_xml()->at(
- "aircraft_exchange_file/requirements_and_specifications/requirements/top_level_aircraft_requirements/study_mission/initial_cruise_mach_number/value");
-
- double weight_fraction = mission_data.get_beta(" climb_to_cruise");
- readPolar clean_config_M_cruise(this->polar(), "clean", "linear", M);
- double thrust_lapse = calculate_thrust_lapse(M, altitude);
- CD_vector = {};
- CL_vector = {};
-
- for (int i = 0; i < W_over_S_data.size(); ++i) {
- CD = clean_config_M_cruise.getCD(W_over_S_data[i], 1.0, weight_fraction, M, altitude);
- CD_vector.push_back(CD);
- }
-
- double minimum_climb_rate = constraint->at("minimum_climb_rate/value");
-
- /* Service ceiling with the minimum climb rate requirement */
- const auto constraint_vector = constraint_analysis_tool.constant_speed_climb(
- W_over_S_data,
- CD_vector,
- weight_fraction,
- thrust_lapse,
- M,
- altitude,
- minimum_climb_rate);
- Constraint Case(W_over_S_data, constraint_vector, constraint->name);
- this->constraint_list.push_back(Case);
- }
- else if (constraint->name == "approach_speed")
- {
-
- }
- else if (constraint->name == "landing_field_length")
- {
- /* Landing */
- double M = constraint->at("Mach/value");
- double thrust_lapse = calculate_thrust_lapse(M, 0.0);
- double CL_max_L = this->aircraft_xml()->at("aircraft_exchange_file/analysis/aerodynamics/lift_coefficients/C_LmaxLanding/value");
- double my_B = this->aircraft_xml()->at("aircraft_exchange_file/requirements_and_specifications/requirements/additional_requirements/landing_gear/braking_coefficient/value");
- double s_G_L = this->aircraft_xml()->at("aircraft_exchange_file/requirements_and_specifications/requirements/top_level_aircraft_requirements/takeoff_distance/value");
- readPolar landing_config(this->polar(), "landing", "linear", M);
- double CD_max_L = landing_config.interpCD(CL_max_L);
- double max_landing_weight = this->aircraft_xml()->at("aircraft_exchange_file/analysis/masses_cg_inertia/maximum_landing_mass/mass_properties/mass/value");
- double MTOM = this->aircraft_xml()->at("aircraft_exchange_file/analysis/masses_cg_inertia/maximum_takeoff_mass/mass_properties/mass/value");
- double weight_fraction = max_landing_weight / MTOM;
- const auto constraint_border = constraint_analysis_tool.braking_roll(CD_max_L, CL_max_L, weight_fraction, thrust_lapse, M, 0.0, my_B, s_G_L);
- this->boundaries.push_back(constraint_border);
- }
- else if (constraint->name == "gust")
- {
- /* Gust */
- double M = constraint->at("Mach/value");
- readPolar landing_config(this->polar(), "landing", "linear", M);
- double weight_fraction = mission_data.get_beta(" takeoff", 1.0); //from mission analysis
- double CL_max_L = this->aircraft_xml()->at("aircraft_exchange_file/analysis/aerodynamics/lift_coefficients/C_LmaxLanding/value");
- double AoA_1 = landing_config.interpalpha(CL_max_L);
- double AoA_2 = landing_config.interpalpha(CL_max_L - CL_max_L * 0.9);
- double CL_alpha = (CL_max_L - (CL_max_L - CL_max_L * 0.9)) / ((AoA_1 - AoA_2) * PI / 180.0);
- double V_TO_L = this->aircraft_xml()->at("aircraft_exchange_file/requirements_and_specifications/requirements/top_level_aircraft_requirements/flight_envelope/maximum_approach_speed/value");
- double w_g = this->configuration()->at("module_configuration_file/program_settings/constraint_selection/standard_set/gust/gust_speed/value");
- double dn_G = this->configuration()->at("module_configuration_file/program_settings/constraint_selection/standard_set/gust/gust_load_factor/value");
- const auto constraint_border = constraint_analysis_tool.gust(CL_alpha, 0.0, V_TO_L, w_g, dn_G, weight_fraction);
- this->boundaries.push_back(constraint_border);
- }
- }
- };
-};
\ No newline at end of file
--
GitLab
From f79aaaf41b4c6e96afe5947dfb81cdd1b385238b Mon Sep 17 00:00:00 2001
From: Meric Taneri <meric.taneri@tum.de>
Date: Mon, 24 Feb 2025 11:14:08 +0100
Subject: [PATCH 4/4] plotting should be fixed
---
constraint_analysis/src/energy_based/energy_based.cpp | 5 +++++
1 file changed, 5 insertions(+)
diff --git a/constraint_analysis/src/energy_based/energy_based.cpp b/constraint_analysis/src/energy_based/energy_based.cpp
index 9a3ea3d7..f82e6b59 100644
--- a/constraint_analysis/src/energy_based/energy_based.cpp
+++ b/constraint_analysis/src/energy_based/energy_based.cpp
@@ -145,6 +145,11 @@ void EnergyBased::operator()()
if (plotting == "true") {
WS_start = 50;
+ WS_end = 820;
+ WS_size = WS_end - WS_start + 1; // Range from 50 to 820, inclusive
+
+ std::vector<double> W_over_S_data(WS_size);
+
for (int i = 0; i < WS_size; ++i) {
W_over_S_data[i] = (WS_start + i) * 9.81;
}
--
GitLab