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 01/10] 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 02/10] 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 03/10] 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 04/10] 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
From a7758272d1be6f3a598a4a81491cc39181b4bb2a Mon Sep 17 00:00:00 2001
From: Meric Taneri <meric.taneri@tum.de>
Date: Mon, 24 Feb 2025 15:30:07 +0100
Subject: [PATCH 05/10] bug-fixes to solve issues caused by disabling the gust
and field length constraints
---
.../constraint_analysis/ca_minfinder.h | 60 +++++++++++--------
.../src/energy_based/energy_based.cpp | 18 ++++--
.../src/energy_based/energy_based.h | 3 +-
3 files changed, 50 insertions(+), 31 deletions(-)
diff --git a/constraint_analysis/include/constraint_analysis/ca_minfinder.h b/constraint_analysis/include/constraint_analysis/ca_minfinder.h
index 67412817..ea0e3777 100644
--- a/constraint_analysis/include/constraint_analysis/ca_minfinder.h
+++ b/constraint_analysis/include/constraint_analysis/ca_minfinder.h
@@ -58,24 +58,27 @@ class Simple_Analysis : Analysis
{
public:
std::vector<double> dominant_thrust_to_weight = {};
- std::vector<double> wing_loading_boundaries = {};
+ double lower_wing_loading_boundary = {};
+ double upper_wing_loading_boundary = {};
std::vector<double> wing_loading_vec = {};
double design_thrust_to_weight = 1000.0;
double design_wing_loading = 0.0;
- Simple_Analysis(std::vector<Constraint> constraint_list, std::vector<double> boundaries) : Analysis(constraint_list)
+ Simple_Analysis(std::vector<Constraint> constraint_list, std::vector<double> lower_boundaries, std::vector<double> upper_boundaries) : Analysis(constraint_list)
{
dominant_thrust_to_weight = constraint_list[0].thrust_to_weight;
wing_loading_vec = constraint_list[0].wing_loading;
- wing_loading_boundaries = boundaries;
+ lower_wing_loading_boundary = *std::max_element(lower_boundaries.begin(), lower_boundaries.end());
+ upper_wing_loading_boundary = *std::min_element(upper_boundaries.begin(), upper_boundaries.end());
};
- Simple_Analysis(std::vector<Constraint> constraint_list, std::shared_ptr<node> acXML, std::vector<double> boundaries) :
+ Simple_Analysis(std::vector<Constraint> constraint_list, std::shared_ptr<node> acXML, std::vector<double> lower_boundaries, std::vector<double> upper_boundaries) :
Analysis(constraint_list, acXML)
{
dominant_thrust_to_weight = constraint_list[0].thrust_to_weight;
wing_loading_vec = constraint_list[0].wing_loading;
- wing_loading_boundaries = boundaries;
+ lower_wing_loading_boundary = *std::max_element(lower_boundaries.begin(), lower_boundaries.end());
+ upper_wing_loading_boundary = *std::min_element(upper_boundaries.begin(), upper_boundaries.end());
};
void find_dominant_curve() {
@@ -90,7 +93,7 @@ public:
void find_design_point(const double& safety_factor)
{
- std::sort(wing_loading_boundaries.begin(), wing_loading_boundaries.end());
+ std::vector<double> wing_loading_boundaries = { lower_wing_loading_boundary, upper_wing_loading_boundary };
std::vector<double> xlims = {};
for (auto x_limit : wing_loading_boundaries) {
auto it = std::distance(wing_loading_vec.begin(), std::lower_bound(wing_loading_vec.begin(), wing_loading_vec.end(), x_limit));
@@ -112,29 +115,38 @@ public:
void find_converged_design_point(const double& previous_TW, const double& previous_WS, const double& safety_factor, const double& buffer_factor)
{
- std::sort(wing_loading_boundaries.begin(), wing_loading_boundaries.end());
- std::vector<double> xlims = {};
- for (auto x_limit : wing_loading_boundaries) {
- auto it = std::distance(wing_loading_vec.begin(), std::lower_bound(wing_loading_vec.begin(), wing_loading_vec.end(), x_limit));
- xlims.push_back(it);
+ std::vector<double> wing_loading_boundaries = { lower_wing_loading_boundary, upper_wing_loading_boundary };
+ if (upper_wing_loading_boundary < lower_wing_loading_boundary)
+ {
+ design_thrust_to_weight = previous_TW;
+ design_wing_loading = previous_WS;
+ std::cout << "infeasable design space, loadings are not adjusted" << std::endl;
}
- std::vector<double> limited_TW(dominant_thrust_to_weight.begin() + xlims[0], dominant_thrust_to_weight.begin() + xlims[1]);
- int i = xlims[0];
- for (double point : limited_TW)
+ else
{
- if (point <= design_thrust_to_weight)
+ std::vector<double> xlims = {};
+ for (auto x_limit : wing_loading_boundaries) {
+ auto it = std::distance(wing_loading_vec.begin(), std::lower_bound(wing_loading_vec.begin(), wing_loading_vec.end(), x_limit));
+ xlims.push_back(it);
+ }
+ std::vector<double> limited_TW(dominant_thrust_to_weight.begin() + xlims[0], dominant_thrust_to_weight.begin() + xlims[1]);
+ int i = xlims[0];
+ for (double point : limited_TW)
{
- design_thrust_to_weight = point;
- design_wing_loading = wing_loading_vec[i];
+ if (point <= design_thrust_to_weight)
+ {
+ design_thrust_to_weight = point;
+ design_wing_loading = wing_loading_vec[i];
+ }
+ i++;
}
- i++;
- }
- design_thrust_to_weight = design_thrust_to_weight * (1.0 + safety_factor);
+ design_thrust_to_weight = design_thrust_to_weight * (1.0 + safety_factor);
- double relative_change_TW = sqrt(pow(1 - previous_TW / design_thrust_to_weight, 2));
- if (relative_change_TW <= buffer_factor/100.0) { design_thrust_to_weight = previous_TW; };
- double relative_change_WS = sqrt(pow(1 - previous_WS / design_wing_loading, 2));
- if (relative_change_WS <= buffer_factor/100.0) { design_wing_loading = previous_WS; };
+ double relative_change_TW = sqrt(pow(1 - previous_TW / design_thrust_to_weight, 2));
+ if (relative_change_TW <= buffer_factor / 100.0) { design_thrust_to_weight = previous_TW; };
+ double relative_change_WS = sqrt(pow(1 - previous_WS / design_wing_loading, 2));
+ if (relative_change_WS <= buffer_factor / 100.0) { design_wing_loading = previous_WS; };
+ };
};
template<typename T>
diff --git a/constraint_analysis/src/energy_based/energy_based.cpp b/constraint_analysis/src/energy_based/energy_based.cpp
index f82e6b59..a72b0a5a 100644
--- a/constraint_analysis/src/energy_based/energy_based.cpp
+++ b/constraint_analysis/src/energy_based/energy_based.cpp
@@ -83,6 +83,9 @@ void EnergyBased::operator()()
W_over_S_data[i] = (WS_start + i) * 9.81;
}
+ this->lower_boundaries.push_back(50.0 * 9.81);
+ this->upper_boundaries.push_back(820.0 * 9.81);
+
Engine engine = Engine(this->engine_directory());
Mattingly::constraint_analysis constraint_analysis_tool;
@@ -113,7 +116,7 @@ void EnergyBased::operator()()
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);
+ Simple_Analysis min_finder = Simple_Analysis(this->constraint_list, this->aircraft_xml(), this->lower_boundaries, this->upper_boundaries);
/* Find dominant curve of the design chart*/
min_finder.find_dominant_curve();
@@ -155,23 +158,26 @@ void EnergyBased::operator()()
}
this->constraint_list = {};
- this->boundaries = {};
+ this->lower_boundaries = { 50. * 9.81 };
+ this->upper_boundaries = { 820. * 9.81 };
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);
+ Simple_Analysis min_finder_2 = Simple_Analysis(this->constraint_list, this->aircraft_xml(), this->lower_boundaries, this->upper_boundaries);
min_finder_2.find_dominant_curve();
ConstraintAnalysis_Plot ca_plot;
+ std::vector<double> boundaries = { min_finder_2.lower_wing_loading_boundary , min_finder_2.upper_wing_loading_boundary };
+
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)
+ for (auto boundary : boundaries)
{
ca_plot.add_curve(boundary, "W/S limit");
}
@@ -347,7 +353,7 @@ void EnergyBased::assess_constraints(std::vector<double>& W_over_S_data, Mattin
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);
+ this->upper_boundaries.push_back(constraint_border);
}
else if (constraint->name == "gust")
{
@@ -363,7 +369,7 @@ void EnergyBased::assess_constraints(std::vector<double>& W_over_S_data, Mattin
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);
+ this->lower_boundaries.push_back(constraint_border);
}
}
};
diff --git a/constraint_analysis/src/energy_based/energy_based.h b/constraint_analysis/src/energy_based/energy_based.h
index b8c56b24..3e48c381 100644
--- a/constraint_analysis/src/energy_based/energy_based.h
+++ b/constraint_analysis/src/energy_based/energy_based.h
@@ -70,7 +70,8 @@ public:
void assess_constraints(std::vector<double>& W_over_S_data, Mattingly::constraint_analysis& constraint_analysis_tool, readMission& mission_data, Engine& engine);
std::vector<Constraint> constraint_list = {};
- std::vector<double> boundaries = {};
+ std::vector<double> lower_boundaries = {};
+ std::vector<double> upper_boundaries = {};
}; // namespace design
#endif // SRC_CONSTRAINT_ANALYSIS_ENERGY_BASED_H_
\ No newline at end of file
--
GitLab
From 814b1d846ee9985e227b5cccda819f40b77f0f8f Mon Sep 17 00:00:00 2001
From: Meric Taneri <meric.taneri@tum.de>
Date: Thu, 27 Feb 2025 08:01:01 +0100
Subject: [PATCH 06/10] approach speed constraint added - will be improved
---
.../constraint_analysis_conf.xml | 16 ++++++++++++--
.../constraint_analysis/ca_functions.h | 1 +
constraint_analysis/src/ca_functions.cpp | 22 +++++++++++++++++++
.../src/energy_based/energy_based.cpp | 17 +++++++++++++-
4 files changed, 53 insertions(+), 3 deletions(-)
diff --git a/constraint_analysis/constraint_analysis_conf.xml b/constraint_analysis/constraint_analysis_conf.xml
index b48efabd..24d27639 100644
--- a/constraint_analysis/constraint_analysis_conf.xml
+++ b/constraint_analysis/constraint_analysis_conf.xml
@@ -5,7 +5,7 @@
<value>CSMR-2020.xml</value>
</aircraft_exchange_file_name>
<aircraft_exchange_file_directory description="Specify the direction in which the aircraft exchange file can be found">
- <value>../projects/CSMR/CSMR-2020/</value>
+ <value>../projects/</value>
</aircraft_exchange_file_directory>
<own_tool_level description="Specify the tool level of this tool">
<value>3</value>
@@ -54,7 +54,7 @@
<value>Energy_Based</value>
</method>
<mode_selection description="The method of sizign, mode_0 changes both T/W and W/S, mode_1 changes only the T/W">
- <value>mode_0</value>
+ <value>mode_1</value>
</mode_selection>
<buffer_factor description="The buffer in percentage, if the rate of change of the design point parameters is less than this percentage, it remains constant">
<value>0.0</value>
@@ -113,6 +113,18 @@
<active>
<value>true</value>
</active>
+ <Mach>
+ <value>0.2</value>
+ </Mach>
+ <free_roll_time>
+ <value>5.0</value>
+ </free_roll_time>
+ <k_TD>
+ <value>1.15</value>
+ </k_TD>
+ <obstacle_height>
+ <value>10</value>
+ </obstacle_height>
</approach_speed>
<landing_field_length>
<active>
diff --git a/constraint_analysis/include/constraint_analysis/ca_functions.h b/constraint_analysis/include/constraint_analysis/ca_functions.h
index a5126853..3c1e611c 100644
--- a/constraint_analysis/include/constraint_analysis/ca_functions.h
+++ b/constraint_analysis/include/constraint_analysis/ca_functions.h
@@ -184,6 +184,7 @@ namespace Mattingly {
*/
auto min_speed(double CL_max, double altitude, double Mach, double weight_fraction) -> const double;
+ auto constraint_analysis::approach_speed(double V_app, double CL_max, double weight_fraction, double s_A, double C_D, double h_obs, double k_TD, double altitude) -> const double;
private:
// Membervariables
diff --git a/constraint_analysis/src/ca_functions.cpp b/constraint_analysis/src/ca_functions.cpp
index 738e8186..43496ff3 100644
--- a/constraint_analysis/src/ca_functions.cpp
+++ b/constraint_analysis/src/ca_functions.cpp
@@ -195,5 +195,27 @@ namespace Mattingly {
return W_over_S;
}
+ auto constraint_analysis::approach_speed(double V_app, double CL_max, double weight_fraction, double s_A, double C_D, double h_obs, double k_TD, double altitude) -> const double {
+
+ double W_over_S;
+
+ double V_stall = 0.0;
+
+ // this is a 4th order polynomial with coefficients (a*x^4 + b*x^3 + c*x^2 + d*x + e = 0)
+ double a = -k_TD * CL_max;
+ double b = 0.;
+ double c = pow(V_app, 2.) * CL_max + 2. * h_obs * CL_max - s_A * k_TD;
+ double d = 0.;
+ double e = -s_A * V_app;
+
+ double del = pow(c, 2.) - 4. * a * e;
+
+ V_stall = std::sqrt((-c + sqrt(del)) / (2. * a));
+ V_stall = std::sqrt((-c - sqrt(del)) / (2. * a));
+
+
+ return pow(V_stall,2.)/(2.*weight_fraction/(atm.getDensity(altitude)*CL_max));
+ }
+
}
diff --git a/constraint_analysis/src/energy_based/energy_based.cpp b/constraint_analysis/src/energy_based/energy_based.cpp
index a72b0a5a..671ebf37 100644
--- a/constraint_analysis/src/energy_based/energy_based.cpp
+++ b/constraint_analysis/src/energy_based/energy_based.cpp
@@ -337,7 +337,22 @@ void EnergyBased::assess_constraints(std::vector<double>& W_over_S_data, Mattin
}
else if (constraint->name == "approach_speed")
{
-
+ double t_FR = constraint->at("free_roll_time/value");
+ double k_TD = constraint->at("k_TD/value");
+ double h_obs = constraint->at("obstacle_height/value");
+ double M = constraint->at("Mach/value");
+ double V_app = this->aircraft_xml()->at(
+ "aircraft_exchange_file/requirements_and_specifications/requirements/top_level_aircraft_requirements/flight_envelope/maximum_approach_speed/value");
+ double CL_max_L = this->aircraft_xml()->at("aircraft_exchange_file/analysis/aerodynamics/lift_coefficients/C_LmaxLanding/value");
+ double s_G_L = this->aircraft_xml()->at("aircraft_exchange_file/requirements_and_specifications/requirements/top_level_aircraft_requirements/takeoff_distance/value");
+ double s_A = s_G_L * 0.4; // find a way to get the touchdown speed
+ 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 approach = constraint_analysis_tool.approach_speed(V_app, CL_max_L, weight_fraction, s_A, CD_max_L, h_obs, k_TD, 0.0);
+ this->upper_boundaries.push_back(approach);
}
else if (constraint->name == "landing_field_length")
{
--
GitLab
From 45fd8d1b26ba25ba1f901b8f29c9bb4ca5b2b5e4 Mon Sep 17 00:00:00 2001
From: Meric Taneri <meric.taneri@tum.de>
Date: Thu, 27 Feb 2025 08:06:11 +0100
Subject: [PATCH 07/10] landing field length adjustment
---
constraint_analysis/src/energy_based/energy_based.cpp | 4 ++--
1 file changed, 2 insertions(+), 2 deletions(-)
diff --git a/constraint_analysis/src/energy_based/energy_based.cpp b/constraint_analysis/src/energy_based/energy_based.cpp
index 671ebf37..1e261450 100644
--- a/constraint_analysis/src/energy_based/energy_based.cpp
+++ b/constraint_analysis/src/energy_based/energy_based.cpp
@@ -344,7 +344,7 @@ void EnergyBased::assess_constraints(std::vector<double>& W_over_S_data, Mattin
double V_app = this->aircraft_xml()->at(
"aircraft_exchange_file/requirements_and_specifications/requirements/top_level_aircraft_requirements/flight_envelope/maximum_approach_speed/value");
double CL_max_L = this->aircraft_xml()->at("aircraft_exchange_file/analysis/aerodynamics/lift_coefficients/C_LmaxLanding/value");
- double s_G_L = this->aircraft_xml()->at("aircraft_exchange_file/requirements_and_specifications/requirements/top_level_aircraft_requirements/takeoff_distance/value");
+ double s_G_L = this->aircraft_xml()->at("aircraft_exchange_file/requirements_and_specifications/requirements/top_level_aircraft_requirements/landing_field_length/value");
double s_A = s_G_L * 0.4; // find a way to get the touchdown speed
readPolar landing_config(this->polar(), "landing", "linear", M);
double CD_max_L = landing_config.interpCD(CL_max_L);
@@ -361,7 +361,7 @@ void EnergyBased::assess_constraints(std::vector<double>& W_over_S_data, Mattin
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");
+ double s_G_L = 0.6*this->aircraft_xml()->at("aircraft_exchange_file/requirements_and_specifications/requirements/top_level_aircraft_requirements/landing_field_length/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");
--
GitLab
From b5d5e95ffd58371d989528bf4f215d69c41d84aa Mon Sep 17 00:00:00 2001
From: Meric Taneri <meric.taneri@tum.de>
Date: Thu, 27 Feb 2025 08:10:46 +0100
Subject: [PATCH 08/10] small fix
---
constraint_analysis/src/ca_functions.cpp | 4 +---
1 file changed, 1 insertion(+), 3 deletions(-)
diff --git a/constraint_analysis/src/ca_functions.cpp b/constraint_analysis/src/ca_functions.cpp
index 43496ff3..8d6e6230 100644
--- a/constraint_analysis/src/ca_functions.cpp
+++ b/constraint_analysis/src/ca_functions.cpp
@@ -210,9 +210,7 @@ namespace Mattingly {
double del = pow(c, 2.) - 4. * a * e;
- V_stall = std::sqrt((-c + sqrt(del)) / (2. * a));
- V_stall = std::sqrt((-c - sqrt(del)) / (2. * a));
-
+ V_stall = std::max(std::sqrt((-c + sqrt(del)) / (2. * a)), std::sqrt((-c - sqrt(del)) / (2. * a)));
return pow(V_stall,2.)/(2.*weight_fraction/(atm.getDensity(altitude)*CL_max));
}
--
GitLab
From 2c3b72f7a1c0c498e6dfa2a90634eebec90cdd4d Mon Sep 17 00:00:00 2001
From: Meric Taneri <meric.taneri@tum.de>
Date: Thu, 27 Feb 2025 10:04:22 +0100
Subject: [PATCH 09/10] fix
---
.../constraint_analysis_conf.xml | 16 ++++------------
constraint_analysis/src/ca_functions.cpp | 6 ++----
.../src/energy_based/energy_based.cpp | 19 +++----------------
3 files changed, 9 insertions(+), 32 deletions(-)
diff --git a/constraint_analysis/constraint_analysis_conf.xml b/constraint_analysis/constraint_analysis_conf.xml
index 24d27639..87827f45 100644
--- a/constraint_analysis/constraint_analysis_conf.xml
+++ b/constraint_analysis/constraint_analysis_conf.xml
@@ -54,7 +54,7 @@
<value>Energy_Based</value>
</method>
<mode_selection description="The method of sizign, mode_0 changes both T/W and W/S, mode_1 changes only the T/W">
- <value>mode_1</value>
+ <value>mode_0</value>
</mode_selection>
<buffer_factor description="The buffer in percentage, if the rate of change of the design point parameters is less than this percentage, it remains constant">
<value>0.0</value>
@@ -109,27 +109,19 @@
<value>0.508</value>
</minimum_climb_rate>
</SEP>
- <approach_speed>
+ <landing_field_length>
<active>
<value>true</value>
</active>
<Mach>
<value>0.2</value>
</Mach>
- <free_roll_time>
- <value>5.0</value>
- </free_roll_time>
<k_TD>
- <value>1.15</value>
+ <value>1.10</value>
</k_TD>
<obstacle_height>
- <value>10</value>
+ <value>15.4</value>
</obstacle_height>
- </approach_speed>
- <landing_field_length>
- <active>
- <value>true</value>
- </active>
<Mach>
<value>0.2</value>
</Mach>
diff --git a/constraint_analysis/src/ca_functions.cpp b/constraint_analysis/src/ca_functions.cpp
index 8d6e6230..a4e9375f 100644
--- a/constraint_analysis/src/ca_functions.cpp
+++ b/constraint_analysis/src/ca_functions.cpp
@@ -195,10 +195,8 @@ namespace Mattingly {
return W_over_S;
}
- auto constraint_analysis::approach_speed(double V_app, double CL_max, double weight_fraction, double s_A, double C_D, double h_obs, double k_TD, double altitude) -> const double {
-
- double W_over_S;
-
+ auto constraint_analysis::approach_speed(double V_app, double CL_max, double weight_fraction, double s_A, double C_D, double h_obs, double k_TD, double altitude) -> const double
+ {
double V_stall = 0.0;
// this is a 4th order polynomial with coefficients (a*x^4 + b*x^3 + c*x^2 + d*x + e = 0)
diff --git a/constraint_analysis/src/energy_based/energy_based.cpp b/constraint_analysis/src/energy_based/energy_based.cpp
index 803ba58e..cbb35926 100644
--- a/constraint_analysis/src/energy_based/energy_based.cpp
+++ b/constraint_analysis/src/energy_based/energy_based.cpp
@@ -335,17 +335,15 @@ void EnergyBased::assess_constraints(std::vector<double>& W_over_S_data, Mattin
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")
{
- double t_FR = constraint->at("free_roll_time/value");
double k_TD = constraint->at("k_TD/value");
double h_obs = constraint->at("obstacle_height/value");
double M = constraint->at("Mach/value");
double V_app = this->aircraft_xml()->at(
"aircraft_exchange_file/requirements_and_specifications/requirements/top_level_aircraft_requirements/flight_envelope/maximum_approach_speed/value");
double CL_max_L = this->aircraft_xml()->at("aircraft_exchange_file/analysis/aerodynamics/lift_coefficients/C_LmaxLanding/value");
- double s_G_L = this->aircraft_xml()->at("aircraft_exchange_file/requirements_and_specifications/requirements/top_level_aircraft_requirements/landing_field_length/value");
- double s_A = s_G_L * 0.4; // find a way to get the touchdown speed
+ double s_A = h_obs / std::tan(3. * std::numbers::pi / 180.);
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");
@@ -353,21 +351,10 @@ void EnergyBased::assess_constraints(std::vector<double>& W_over_S_data, Mattin
double weight_fraction = max_landing_weight / MTOM;
const auto approach = constraint_analysis_tool.approach_speed(V_app, CL_max_L, weight_fraction, s_A, CD_max_L, h_obs, k_TD, 0.0);
this->upper_boundaries.push_back(approach);
- }
- 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 = 0.6 * this->aircraft_xml()->at("aircraft_exchange_file/requirements_and_specifications/requirements/top_level_aircraft_requirements/landing_field_length/value");
- // s_G_L as braking distance = 60% of landing field length
- 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;
+ double s_G_L = 0.6 * this->aircraft_xml()->at("aircraft_exchange_file/requirements_and_specifications/requirements/top_level_aircraft_requirements/landing_field_length/value") - s_A;
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->upper_boundaries.push_back(constraint_border);
}
--
GitLab
From 06ce5c2aa56a1c6ddf606a53f560042874905b3e Mon Sep 17 00:00:00 2001
From: Meric Taneri <meric.taneri@tum.de>
Date: Thu, 27 Feb 2025 12:21:58 +0100
Subject: [PATCH 10/10] free roll added
---
.../constraint_analysis_conf.xml | 7 +++++--
.../include/constraint_analysis/ca_functions.h | 5 ++++-
constraint_analysis/src/ca_functions.cpp | 18 ++++++++++++++++--
.../src/energy_based/energy_based.cpp | 6 +++++-
4 files changed, 30 insertions(+), 6 deletions(-)
diff --git a/constraint_analysis/constraint_analysis_conf.xml b/constraint_analysis/constraint_analysis_conf.xml
index 87827f45..5d913201 100644
--- a/constraint_analysis/constraint_analysis_conf.xml
+++ b/constraint_analysis/constraint_analysis_conf.xml
@@ -5,7 +5,7 @@
<value>CSMR-2020.xml</value>
</aircraft_exchange_file_name>
<aircraft_exchange_file_directory description="Specify the direction in which the aircraft exchange file can be found">
- <value>../projects/</value>
+ <value>../projects2/</value>
</aircraft_exchange_file_directory>
<own_tool_level description="Specify the tool level of this tool">
<value>3</value>
@@ -117,8 +117,11 @@
<value>0.2</value>
</Mach>
<k_TD>
- <value>1.10</value>
+ <value>1.15</value>
</k_TD>
+ <t_FR>
+ <value>3</value>
+ </t_FR>
<obstacle_height>
<value>15.4</value>
</obstacle_height>
diff --git a/constraint_analysis/include/constraint_analysis/ca_functions.h b/constraint_analysis/include/constraint_analysis/ca_functions.h
index 3c1e611c..0dab1748 100644
--- a/constraint_analysis/include/constraint_analysis/ca_functions.h
+++ b/constraint_analysis/include/constraint_analysis/ca_functions.h
@@ -184,7 +184,10 @@ namespace Mattingly {
*/
auto min_speed(double CL_max, double altitude, double Mach, double weight_fraction) -> const double;
- auto constraint_analysis::approach_speed(double V_app, double CL_max, double weight_fraction, double s_A, double C_D, double h_obs, double k_TD, double altitude) -> const double;
+ auto approach_speed(double V_app, double CL_max, double weight_fraction, double s_A, double C_D, double h_obs, double k_TD, double altitude) -> const double;
+
+ auto get_V_stall(double V_app, double k_TD, double CL_max, double CD, double h_obs, double s_A) -> const double;
+
private:
// Membervariables
diff --git a/constraint_analysis/src/ca_functions.cpp b/constraint_analysis/src/ca_functions.cpp
index a4e9375f..f176ebd0 100644
--- a/constraint_analysis/src/ca_functions.cpp
+++ b/constraint_analysis/src/ca_functions.cpp
@@ -204,14 +204,28 @@ namespace Mattingly {
double b = 0.;
double c = pow(V_app, 2.) * CL_max + 2. * h_obs * CL_max - s_A * k_TD;
double d = 0.;
- double e = -s_A * V_app;
+ double e = -s_A * V_app * C_D;
double del = pow(c, 2.) - 4. * a * e;
V_stall = std::max(std::sqrt((-c + sqrt(del)) / (2. * a)), std::sqrt((-c - sqrt(del)) / (2. * a)));
-
+
return pow(V_stall,2.)/(2.*weight_fraction/(atm.getDensity(altitude)*CL_max));
}
+ auto constraint_analysis::get_V_stall(double V_app, double k_TD, double CL_max, double CD, double h_obs, double s_A) -> const double
+ {
+ // this is a 4th order polynomial with coefficients (a*x^4 + b*x^3 + c*x^2 + d*x + e = 0)
+ double a = -k_TD * CL_max;
+ double b = 0.;
+ double c = pow(V_app, 2.) * CL_max + 2. * h_obs * CL_max - s_A * k_TD;
+ double d = 0.;
+ double e = -s_A * V_app * CD;
+
+ double del = pow(c, 2.) - 4. * a * e;
+
+ return std::max(std::sqrt((-c + sqrt(del)) / (2. * a)), std::sqrt((-c - sqrt(del)) / (2. * a)));
+ }
+
}
diff --git a/constraint_analysis/src/energy_based/energy_based.cpp b/constraint_analysis/src/energy_based/energy_based.cpp
index cbb35926..eb6c05e7 100644
--- a/constraint_analysis/src/energy_based/energy_based.cpp
+++ b/constraint_analysis/src/energy_based/energy_based.cpp
@@ -340,6 +340,7 @@ void EnergyBased::assess_constraints(std::vector<double>& W_over_S_data, Mattin
double k_TD = constraint->at("k_TD/value");
double h_obs = constraint->at("obstacle_height/value");
double M = constraint->at("Mach/value");
+ double t_FR = constraint->at("t_FR/value");
double V_app = this->aircraft_xml()->at(
"aircraft_exchange_file/requirements_and_specifications/requirements/top_level_aircraft_requirements/flight_envelope/maximum_approach_speed/value");
double CL_max_L = this->aircraft_xml()->at("aircraft_exchange_file/analysis/aerodynamics/lift_coefficients/C_LmaxLanding/value");
@@ -351,10 +352,13 @@ void EnergyBased::assess_constraints(std::vector<double>& W_over_S_data, Mattin
double weight_fraction = max_landing_weight / MTOM;
const auto approach = constraint_analysis_tool.approach_speed(V_app, CL_max_L, weight_fraction, s_A, CD_max_L, h_obs, k_TD, 0.0);
this->upper_boundaries.push_back(approach);
+ double V_stall = constraint_analysis_tool.get_V_stall(V_app, k_TD, CL_max_L, CD_max_L, h_obs, s_A);
+ double s_FR = V_stall * k_TD * t_FR;
+
/* Landing */
double thrust_lapse = calculate_thrust_lapse(M, 0.0);
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 = 0.6 * this->aircraft_xml()->at("aircraft_exchange_file/requirements_and_specifications/requirements/top_level_aircraft_requirements/landing_field_length/value") - s_A;
+ double s_G_L = 0.6 * this->aircraft_xml()->at("aircraft_exchange_file/requirements_and_specifications/requirements/top_level_aircraft_requirements/landing_field_length/value") - s_A - s_FR;
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->upper_boundaries.push_back(constraint_border);
}
--
GitLab