From b2227eb2d6924453cd00641d80b0fb14fcbaaf50 Mon Sep 17 00:00:00 2001 From: tshibuk Date: Fri, 2 Jun 2023 18:47:43 +0900 Subject: [PATCH 1/8] add assertions, change solar calculations --- src/dynamics/thermal/heater.cpp | 7 +++++ src/dynamics/thermal/heater.hpp | 6 ++++ src/dynamics/thermal/heater_controller.cpp | 7 ++++- src/dynamics/thermal/heater_controller.hpp | 7 +++++ src/dynamics/thermal/heatload.cpp | 13 +++++++-- src/dynamics/thermal/heatload.hpp | 6 ++++ src/dynamics/thermal/initialize_heater.cpp | 19 +++++++++++- src/dynamics/thermal/initialize_heatload.cpp | 27 +++++++++++++++-- src/dynamics/thermal/initialize_node.cpp | 24 ++++++++------- .../thermal/initialize_temperature.cpp | 10 +++++++ src/dynamics/thermal/node.cpp | 29 ++++++++++++------- src/dynamics/thermal/node.hpp | 10 +++++-- src/dynamics/thermal/temperature.cpp | 25 ++++++++-------- src/dynamics/thermal/temperature.hpp | 12 ++++---- src/environment/global/physical_constants.hpp | 1 + 15 files changed, 155 insertions(+), 48 deletions(-) diff --git a/src/dynamics/thermal/heater.cpp b/src/dynamics/thermal/heater.cpp index 2da6e158f..37b956121 100644 --- a/src/dynamics/thermal/heater.cpp +++ b/src/dynamics/thermal/heater.cpp @@ -1,10 +1,12 @@ #include "heater.hpp" +#include #include using namespace std; Heater::Heater(const unsigned int heater_id, const double power_rating_W) : heater_id_(heater_id), power_rating_W_(power_rating_W) { + AssertHeaterParams(); heater_status_ = HeaterStatus::kOff; power_output_W_ = 0; } @@ -26,3 +28,8 @@ void Heater::PrintParam(void) { cout << " status : " << static_cast(heater_status_) << endl; cout << " power output : " << power_output_W_ << endl; } + +void Heater::AssertHeaterParams(void) { + assert(heater_id_ >= 1); // Heater ID must be larger than 1 + assert(power_rating_W_ >= 0); // Power rating must be larger than 0 +} diff --git a/src/dynamics/thermal/heater.hpp b/src/dynamics/thermal/heater.hpp index a923fa52d..6472a6615 100644 --- a/src/dynamics/thermal/heater.hpp +++ b/src/dynamics/thermal/heater.hpp @@ -83,6 +83,12 @@ class Heater { * @brief Print power_rating_W, heater_status_, power_output_W_ */ void PrintParam(void); + + /** + * @fn AssertHeaterParams + * @brief Check Heater Parameters + */ + void AssertHeaterParams(void); }; #endif // S2E_DYNAMICS_THERMAL_HEATER_HPP_ diff --git a/src/dynamics/thermal/heater_controller.cpp b/src/dynamics/thermal/heater_controller.cpp index 8bb9e6db0..10280528d 100644 --- a/src/dynamics/thermal/heater_controller.cpp +++ b/src/dynamics/thermal/heater_controller.cpp @@ -1,11 +1,14 @@ #include "heater_controller.hpp" +#include #include using namespace std; HeaterController::HeaterController(const double lower_threshold_degC, const double upper_threshold_degC) - : lower_threshold_degC_(lower_threshold_degC), upper_threshold_degC_(upper_threshold_degC) {} + : lower_threshold_degC_(lower_threshold_degC), upper_threshold_degC_(upper_threshold_degC) { + AssertHeaterControllerParams(); +} HeaterController::~HeaterController() {} @@ -21,3 +24,5 @@ void HeaterController::ControlHeater(Heater* p_heater, double temperature_degC) } } } + +void HeaterController::AssertHeaterControllerParams() { assert(upper_threshold_degC_ > lower_threshold_degC_); } diff --git a/src/dynamics/thermal/heater_controller.hpp b/src/dynamics/thermal/heater_controller.hpp index b678d6da1..db271ec84 100644 --- a/src/dynamics/thermal/heater_controller.hpp +++ b/src/dynamics/thermal/heater_controller.hpp @@ -10,6 +10,7 @@ #include #include #include + #include "heater.hpp" class HeaterController { @@ -72,6 +73,12 @@ class HeaterController { * @param[in] temperature_degC */ void ControlHeater(Heater* p_heater, double temperature_degC); + + /** + * @fn AssertHeaterControllerParams + * @brief Check Parameters of HeaterController + */ + void AssertHeaterControllerParams(void); }; #endif // S2E_DYNAMICS_THERMAL_HEATER_CONTROLLER_HPP_ diff --git a/src/dynamics/thermal/heatload.cpp b/src/dynamics/thermal/heatload.cpp index 6aa33b67f..b0ac5e106 100644 --- a/src/dynamics/thermal/heatload.cpp +++ b/src/dynamics/thermal/heatload.cpp @@ -1,5 +1,6 @@ #include "heatload.hpp" +#include #include #include @@ -15,8 +16,8 @@ Heatload::Heatload(int node_id, std::vector time_table_s, std::vector= 1); + assert(internal_heatload_table_W_.size() >= 1); + // size of time_table_s_ and internal_heatload_table_W_ must be the same + assert(time_table_s_.size() == internal_heatload_table_W_.size()); +} diff --git a/src/dynamics/thermal/heatload.hpp b/src/dynamics/thermal/heatload.hpp index ad8272d70..805b1456b 100644 --- a/src/dynamics/thermal/heatload.hpp +++ b/src/dynamics/thermal/heatload.hpp @@ -99,6 +99,12 @@ class Heatload { * @param[in] heater_heatload_W */ inline void SetHeaterHeatload_W(double heater_heatload_W) { heater_heatload_W_ = heater_heatload_W; } + + /** + * @fn AssertHeatloadParams + * @brief Check if Heatload Parameters are Correct + */ + void AssertHeatloadParams(); }; #endif // S2E_DYNAMICS_THERMAL_HEATLOAD_HPP_ diff --git a/src/dynamics/thermal/initialize_heater.cpp b/src/dynamics/thermal/initialize_heater.cpp index c00c7fc19..fb488ad54 100644 --- a/src/dynamics/thermal/initialize_heater.cpp +++ b/src/dynamics/thermal/initialize_heater.cpp @@ -5,12 +5,26 @@ #include "initialize_heater.hpp" -// [FIXME] Write discription of csv file for heater initialization +#include + +/* Import heater properties by reading CSV File (heaters.csv) + +[File Formats of heater.csv] +column 1: Heater_id(int, Use values larger than or equal to 1) +column 2: Power Rating (double, [W]) +column 3: Lower threshold of control (double, [degC]) +column 4: Upper threshold of control (double, [degC]) + +First row is for Header, data begins from the second row +*/ Heater InitHeater(const std::vector& heater_str) { using std::stod; using std::stoi; + size_t heater_str_size_defined = 4; // Correct size of heater_str + assert(heater_str.size() == heater_str_size_defined); // Check if size of heater_str is correct + int heater_id = 0; double power_rating_W = 0; // [W] @@ -28,6 +42,9 @@ Heater InitHeater(const std::vector& heater_str) { HeaterController InitHeaterController(const std::vector& heater_str) { using std::stod; + size_t heater_str_size_defined = 4; // Correct size of heater_str + assert(heater_str.size() == heater_str_size_defined); // Check if size of heater_str is correct + // Index to read from heater_str for each parameter int index_lower_threshold = 2; int index_upper_threshold = 3; diff --git a/src/dynamics/thermal/initialize_heatload.cpp b/src/dynamics/thermal/initialize_heatload.cpp index c76b79be5..888353593 100644 --- a/src/dynamics/thermal/initialize_heatload.cpp +++ b/src/dynamics/thermal/initialize_heatload.cpp @@ -5,13 +5,36 @@ #include "initialize_heatload.hpp" -// [FIXME] Write discription of csv file for heatload initialization +#include + +/* Import heatload properties by reading CSV File (heatload.csv) + +[File Formats of heatload.csv] +First Row: Array of Simulation Times [s] +After Second Row: + First Column: Node ID + After Second Column: Array of Heatload Values [W] + +Heatload object calculates heatload value of certain time by interpolating values of the heatload array + +Data Example: +NodeID/Time, 0, 99, 100, 300, 500 +0, 10, 10, 20, 10, 20 + +In this example, heatload at time 200s is calculated as 15W from interpolation + +The number of nodes and heatloads must be the same + +*/ Heatload InitHeatload(const std::vector& time_str, const std::vector& internal_heatload_str) { using std::stod; using std::stoi; - // [FIXME] Include assertion for size of time_str and internal_heatload_str to be larger than 2 + assert(time_str.size() >= 2); // Number of columns must be larger than 2 + assert(internal_heatload_str.size() >= 2); // Number of columns must be larger than 2 + assert(time_str.size() == internal_heatload_str.size()); // Number of columns must be same + std::vector time_table_s(time_str.size() - 1); // exclude index std::vector internal_heatload_table_W(internal_heatload_str.size() - 1); // exclude index diff --git a/src/dynamics/thermal/initialize_node.cpp b/src/dynamics/thermal/initialize_node.cpp index 920a10225..fce137243 100644 --- a/src/dynamics/thermal/initialize_node.cpp +++ b/src/dynamics/thermal/initialize_node.cpp @@ -5,15 +5,16 @@ #include "initialize_node.hpp" +#include #include #include #include #include #include -/* Import node properties by reading CSV File (Node.csv) +/* Import node properties by reading CSV File (node.csv) -[File Formats of Node.csv] +[File Formats of node.csv] column 1: Node_id(int) column 2: Node_label(string) column 3: Node_type(int, 0: diffusive, 1: boundary), Arithmetic node to be implemented as future work @@ -38,14 +39,17 @@ Node InitNode(const std::vector& node_str) { using std::stod; using std::stoi; - int node_id = 0; // node number - std::string node_label = "temp"; // node name - int node_type_int = 0; // node type - int heater_id = 0; // heater node index - double temperature_K = 0; // [K] - double capacity_J_K = 0; // [J/K] - double alpha = 0; // [] - double area_m2 = 0; // [m^2] + size_t node_str_size_defined = 11; // Correct size of node_str + assert(node_str.size() == node_str_size_defined); // Check if size of node_str is correct + + int node_id = 0; // node number + std::string node_label = "temp"; // node name + int node_type_int = 0; // node type + int heater_id = 0; // heater node index + double temperature_K = 0; // [K] + double capacity_J_K = 0; // [J/K] + double alpha = 0; // [] + double area_m2 = 0; // [m^2] // Index to read from node_str for each parameter int index_node_id = 0; diff --git a/src/dynamics/thermal/initialize_temperature.cpp b/src/dynamics/thermal/initialize_temperature.cpp index 3508961c2..79c34b01c 100644 --- a/src/dynamics/thermal/initialize_temperature.cpp +++ b/src/dynamics/thermal/initialize_temperature.cpp @@ -5,6 +5,7 @@ #include "initialize_temperature.hpp" +#include #include #include #include @@ -70,6 +71,7 @@ Temperature* InitTemperature(const std::string file_name, const double rk_prop_s vector> heater_str_list; // string vector of heater property data vector> heatload_str_list; // string vector of heatload property data unsigned int node_num = 1; + unsigned int heatload_num = 1; unsigned int heater_num = 1; bool is_calc_enabled = mainIni.ReadEnable("THERMAL", "calculation"); @@ -100,6 +102,7 @@ Temperature* InitTemperature(const std::string file_name, const double rk_prop_s /*since we don't know the number of node_list yet, set node_num=100 temporary. Recall that Nodes_num are given to this function only to reseve memory*/ + heatload_num = heatload_str_list.size() - 1; auto times_itr = heatload_str_list.begin(); // First Row is Time Data for (auto itr = heatload_str_list.begin() + 1; itr != heatload_str_list.end(); ++itr) { heatload_list.push_back(InitHeatload(*times_itr, *itr)); @@ -118,6 +121,8 @@ Temperature* InitTemperature(const std::string file_name, const double rk_prop_s node_list.push_back(InitNode(*itr)); } + assert(node_num == heatload_num); // Number of nodes and heatload lists must be the same + // Read Heater Properties from CSV File string filepath_heater = file_path + "heaters.csv"; IniAccess conf_heater(filepath_heater); @@ -141,6 +146,11 @@ Temperature* InitTemperature(const std::string file_name, const double rk_prop_s conf_cij.ReadCsvDoubleWithHeader(conductance_matrix, node_num, 1, 1); conf_rij.ReadCsvDoubleWithHeader(radiation_matrix, node_num, 1, 1); + assert(conductance_matrix.size() == node_num); // Dimension must be same as node_num + assert(radiation_matrix.size() == node_num); // Dimension must be same as node_num + assert(conductance_matrix.size() == conductance_matrix[0].size()); // Must be square matrix + assert(radiation_matrix.size() == radiation_matrix[0].size()); // Must be square matrix + Temperature* temperature; temperature = new Temperature(conductance_matrix, radiation_matrix, node_list, heatload_list, heater_list, heater_controller_list, node_num, rk_prop_step_s, is_calc_enabled, solar_calc_setting, debug); diff --git a/src/dynamics/thermal/node.cpp b/src/dynamics/thermal/node.cpp index 6eef4c865..f865bbdb2 100644 --- a/src/dynamics/thermal/node.cpp +++ b/src/dynamics/thermal/node.cpp @@ -5,6 +5,7 @@ #include "node.hpp" +#include #include #include @@ -20,25 +21,22 @@ Node::Node(const int node_id, const string node_name, const NodeType node_type, capacity_J_K_(capacity_J_K), alpha_(alpha), area_m2_(area_m2), - normal_vector_b_(normal_vector_b), - node_type_(node_type) { + node_type_(node_type), + normal_vector_b_(normal_vector_b) { + AssertNodeParams(); solar_radiation_W_ = 0; } Node::~Node() {} -double Node::CalcSolarRadiation_W(libra::Vector<3> sun_direction_b) { - // FIXME: constants - double R = 6.96E+8; // Distance from sun - double T = 5778; // sun surface temperature [K] - double sigma = 5.67E-8; // stephan-boltzman constant - double a = sun_direction_b.CalcNorm(); // distance from sun - double S = pow((R / a), 2) * sigma * pow(T, 4); // Solar Constant at s/c position S[W/m^2] - double cos_theta = InnerProduct(sun_direction_b, normal_vector_b_) / a / normal_vector_b_.CalcNorm(); +double Node::CalcSolarRadiation_W(libra::Vector<3> sun_position_b_m) { + double sun_distance_m = sun_position_b_m.CalcNorm(); + double solar_flux_W_m2 = environment::solar_constant_W_m2 * pow((environment::astronomical_unit_m / sun_distance_m), 2); + double cos_theta = InnerProduct(sun_position_b_m, normal_vector_b_) / sun_distance_m / normal_vector_b_.CalcNorm(); // calculate sun_power if (cos_theta > 0) - solar_radiation_W_ = S * area_m2_ * alpha_ * cos_theta; + solar_radiation_W_ = solar_flux_W_m2 * area_m2_ * alpha_ * cos_theta; else solar_radiation_W_ = 0; return solar_radiation_W_; @@ -66,3 +64,12 @@ void Node::PrintParam(void) { } cout << endl; } + +void Node::AssertNodeParams(void) { + assert(node_id_ >= 0); // Node ID should be larger than 0 + assert(heater_id_ >= 0); // Heater ID should be larger than 0 + assert(temperature_K_ >= environment::absolute_zero_degC); // Temperature must be larger than zero kelvin + assert(capacity_J_K_ >= 0); // Capacity must be larger than 0, use 0 when node is boundary or arithmetic + assert(0 <= alpha_ && alpha_ <= 1); // alpha must be between 0 and 1 + assert(area_m2_ >= 0); // Area must be larger than 0 +} diff --git a/src/dynamics/thermal/node.hpp b/src/dynamics/thermal/node.hpp index 0b96e506e..84d251c1b 100644 --- a/src/dynamics/thermal/node.hpp +++ b/src/dynamics/thermal/node.hpp @@ -67,10 +67,10 @@ class Node { * @fn CalcSolarRadiation_W * @brief Calculate solar radiation [W] from sun direction, alpha, area, and normal vector * - * @param sun_direction_b: Sun direction in body frame + * @param sun_position_b_m: Sun position in body frame [m] * @return double: Solar Radiation [W] */ - double CalcSolarRadiation_W(libra::Vector<3> sun_direction_b); // 太陽入射熱を計算 + double CalcSolarRadiation_W(libra::Vector<3> sun_position_b_m); // Getter /** @@ -137,6 +137,12 @@ class Node { * @brief Print parameters of node in debug output */ void PrintParam(void); + + /** + * @fn AssertNodeParams + * @brief Check if Node Parameters are Correct + */ + void AssertNodeParams(void); }; #endif // S2E_DYNAMICS_THERMAL_NODE_HPP_ diff --git a/src/dynamics/thermal/temperature.cpp b/src/dynamics/thermal/temperature.cpp index 74cccba2c..86c7cdab1 100644 --- a/src/dynamics/thermal/temperature.cpp +++ b/src/dynamics/thermal/temperature.cpp @@ -44,13 +44,13 @@ Temperature::Temperature() { Temperature::~Temperature() {} -void Temperature::Propagate(libra::Vector<3> sun_direction_b, const double time_end_s) { +void Temperature::Propagate(libra::Vector<3> sun_position_b_m, const double time_end_s) { if (!is_calc_enabled_) return; while (time_end_s - propagation_time_s_ - propagation_step_s_ > 1.0e-6) { - CalcRungeOneStep(propagation_time_s_, propagation_step_s_, sun_direction_b, node_num_); + CalcRungeOneStep(propagation_time_s_, propagation_step_s_, sun_position_b_m, node_num_); propagation_time_s_ += propagation_step_s_; } - CalcRungeOneStep(propagation_time_s_, time_end_s - propagation_time_s_, sun_direction_b, node_num_); + CalcRungeOneStep(propagation_time_s_, time_end_s - propagation_time_s_, sun_position_b_m, node_num_); propagation_time_s_ = time_end_s; UpdateHeaterStatus(); @@ -65,9 +65,9 @@ void Temperature::Propagate(libra::Vector<3> sun_direction_b, const double time_ cout << setprecision(4) << itr->GetSolarRadiation_W() << " "; } cout << "SunDir: "; - double norm_sun_direction_b = sun_direction_b.CalcNorm(); + double sun_distance_m = sun_position_b_m.CalcNorm(); for (int i = 0; i < 3; i++) { - cout << setprecision(3) << sun_direction_b[i] / norm_sun_direction_b << " "; + cout << setprecision(3) << sun_position_b_m[i] / sun_distance_m << " "; } cout << "Heatload: "; for (auto itr = heatloads_.begin(); itr != heatloads_.end(); ++itr) { @@ -77,7 +77,7 @@ void Temperature::Propagate(libra::Vector<3> sun_direction_b, const double time_ } } -void Temperature::CalcRungeOneStep(double time_now_s, double time_step_s, libra::Vector<3> sun_direction_b, int node_num) { +void Temperature::CalcRungeOneStep(double time_now_s, double time_step_s, libra::Vector<3> sun_position_b_m, int node_num) { vector temperatures_now_K(node_num); for (int i = 0; i < node_num; i++) { temperatures_now_K[i] = nodes_[i].GetTemperature_K(); @@ -86,22 +86,22 @@ void Temperature::CalcRungeOneStep(double time_now_s, double time_step_s, libra: vector k1(node_num), k2(node_num), k3(node_num), k4(node_num); vector xk2(node_num), xk3(node_num), xk4(node_num); - k1 = CalcTemperatureDifferentials(temperatures_now_K, time_now_s, sun_direction_b, node_num); + k1 = CalcTemperatureDifferentials(temperatures_now_K, time_now_s, sun_position_b_m, node_num); for (int i = 0; i < node_num; i++) { xk2[i] = temperatures_now_K[i] + (time_step_s / 2.0) * k1[i]; } - k2 = CalcTemperatureDifferentials(xk2, (time_now_s + time_step_s / 2.0), sun_direction_b, node_num); + k2 = CalcTemperatureDifferentials(xk2, (time_now_s + time_step_s / 2.0), sun_position_b_m, node_num); for (int i = 0; i < node_num; i++) { xk3[i] = temperatures_now_K[i] + (time_step_s / 2.0) * k2[i]; } - k3 = CalcTemperatureDifferentials(xk3, (time_now_s + time_step_s / 2.0), sun_direction_b, node_num); + k3 = CalcTemperatureDifferentials(xk3, (time_now_s + time_step_s / 2.0), sun_position_b_m, node_num); for (int i = 0; i < node_num; i++) { xk4[i] = temperatures_now_K[i] + time_step_s * k3[i]; } - k4 = CalcTemperatureDifferentials(xk4, (time_now_s + time_step_s), sun_direction_b, node_num); + k4 = CalcTemperatureDifferentials(xk4, (time_now_s + time_step_s), sun_position_b_m, node_num); vector temperatures_next_K(node_num); for (int i = 0; i < node_num; i++) { @@ -113,7 +113,7 @@ void Temperature::CalcRungeOneStep(double time_now_s, double time_step_s, libra: } } -vector Temperature::CalcTemperatureDifferentials(vector temperatures_K, double t, libra::Vector<3> sun_direction, int node_num) { +vector Temperature::CalcTemperatureDifferentials(vector temperatures_K, double t, libra::Vector<3> sun_position_b_m, int node_num) { // TODO: consider the following unused arguments are really needed UNUSED(temperatures_K); @@ -122,7 +122,7 @@ vector Temperature::CalcTemperatureDifferentials(vector temperat heatloads_[i].SetElapsedTime_s(t); if (nodes_[i].GetNodeType() == NodeType::kDiffusive) { if (solar_calc_setting_ == SolarCalcSetting::kEnable) { - double solar_radiation_W = nodes_[i].CalcSolarRadiation_W(sun_direction); + double solar_radiation_W = nodes_[i].CalcSolarRadiation_W(sun_position_b_m); heatloads_[i].SetSolarHeatload_W(solar_radiation_W); } double heater_power_W = GetHeaterPower_W(i); @@ -157,7 +157,6 @@ double Temperature::GetHeaterPower_W(int node_id) { } void Temperature::UpdateHeaterStatus(void) { - // [FIXME] Heater status doesn't get updated... for (auto itr = nodes_.begin(); itr != nodes_.end(); ++itr) { int heater_id = itr->GetHeaterId(); if (heater_id > 0) { diff --git a/src/dynamics/thermal/temperature.hpp b/src/dynamics/thermal/temperature.hpp index 72962a41b..cb6a36133 100644 --- a/src/dynamics/thermal/temperature.hpp +++ b/src/dynamics/thermal/temperature.hpp @@ -49,21 +49,21 @@ class Temperature : public ILoggable { * * @param[in] time_now_s: Current elapsed time [s] * @param[in] time_step_s: Time step of RK4 [s] - * @param[in] sun_direction_b: Sun direction in body frame + * @param[in] sun_position_b_m: Sun position in body frame [m] * @param[in] node_num: Number of nodes */ - void CalcRungeOneStep(double time_now_s, double time_step_s, libra::Vector<3> sun_direction_b, int node_num); + void CalcRungeOneStep(double time_now_s, double time_step_s, libra::Vector<3> sun_position_b_m, int node_num); /** * @fn CalcTemperatureDifferentials * @brief Calculate differential of thermal equilibrium equation * * @param temperatures_K: [UNUSED] Temperatures of each node [K] * @param time_now_s: Current elapsed time [s] - * @param sun_direction_b: Sun direction in body frame + * @param[in] sun_position_b_m: Sun position in body frame [m] * @param node_num: Number of nodes * @return std::vector: Differential of thermal equilibrium equation at time now */ - std::vector CalcTemperatureDifferentials(std::vector temperatures_K, double time_now_s, const libra::Vector<3> sun_direction_b, + std::vector CalcTemperatureDifferentials(std::vector temperatures_K, double time_now_s, const libra::Vector<3> sun_position_b_m, int node_num); public: @@ -103,10 +103,10 @@ class Temperature : public ILoggable { * @fn Propagate * @brief Propagate thermal calculation until time_end_s * - * @param sun_direction_b: Sun direction in body frame + * @param[in] sun_position_b_m: Sun position in body frame [m] * @param time_end_s: Time to finish propagation [s] */ - void Propagate(libra::Vector<3> sun_direction_b, const double time_end_s); + void Propagate(libra::Vector<3> sun_position_b_m, const double time_end_s); // Getter /** diff --git a/src/environment/global/physical_constants.hpp b/src/environment/global/physical_constants.hpp index 40a79daea..c2acd7de0 100644 --- a/src/environment/global/physical_constants.hpp +++ b/src/environment/global/physical_constants.hpp @@ -23,6 +23,7 @@ DEFINE_PHYSICAL_CONSTANT(speed_of_light_m_s, 299792458.0L) //! DEFINE_PHYSICAL_CONSTANT(boltzmann_constant_J_K, 1.380649e-23L) //!< Boltzmann constant [J/K] DEFINE_PHYSICAL_CONSTANT(stefan_boltzmann_constant_W_m2K4, 5.670374419e-8L) //!< Stefan-Boltzmann constant [W/m2K4] DEFINE_PHYSICAL_CONSTANT(absolute_zero_degC, -273.15L) //!< Absolute zero [degC] +DEFINE_PHYSICAL_CONSTANT(solar_constant_W_m2, 1.366e3L) //!< Solar constant [W/m2] } // namespace physics inline namespace astronomy { From 8b8742b8cf8eff29cf55371731160eb106cbd96b Mon Sep 17 00:00:00 2001 From: tshibuk Date: Tue, 6 Jun 2023 14:03:20 +0900 Subject: [PATCH 2/8] refactor Node::CalcSolarRadiation_W based on SolarRadiationPressureEnvironment::GetPowerDensity_W_m2() --- src/dynamics/dynamics.cpp | 3 ++- src/dynamics/thermal/initialize_temperature.cpp | 4 ++-- src/dynamics/thermal/initialize_temperature.hpp | 2 +- src/dynamics/thermal/node.cpp | 3 +-- src/dynamics/thermal/node.hpp | 2 +- src/dynamics/thermal/temperature.cpp | 9 ++++++--- src/dynamics/thermal/temperature.hpp | 14 ++++++++------ .../solar_radiation_pressure_environment.cpp | 2 ++ .../solar_radiation_pressure_environment.hpp | 17 ++++++++++++----- 9 files changed, 35 insertions(+), 21 deletions(-) diff --git a/src/dynamics/dynamics.cpp b/src/dynamics/dynamics.cpp index 10c403e98..04d0886a7 100644 --- a/src/dynamics/dynamics.cpp +++ b/src/dynamics/dynamics.cpp @@ -29,7 +29,8 @@ void Dynamics::Initialize(const SimulationConfiguration* simulation_configuratio local_celestial_information.GetGlobalInformation().GetCenterBodyGravityConstant_m3_s2(), "ORBIT", relative_information); attitude_ = InitAttitude(simulation_configuration->spacecraft_file_list_[spacecraft_id], orbit_, &local_celestial_information, simulation_time->GetAttitudeRkStepTime_s(), structure->GetKinematicsParameters().GetInertiaTensor_b_kgm2(), spacecraft_id); - temperature_ = InitTemperature(simulation_configuration->spacecraft_file_list_[spacecraft_id], simulation_time->GetThermalRkStepTime_s()); + temperature_ = InitTemperature(simulation_configuration->spacecraft_file_list_[spacecraft_id], simulation_time->GetThermalRkStepTime_s(), + &(local_environment_->GetSolarRadiationPressure())); // To get initial value orbit_->UpdateByAttitude(attitude_->GetQuaternion_i2b()); diff --git a/src/dynamics/thermal/initialize_temperature.cpp b/src/dynamics/thermal/initialize_temperature.cpp index 79c34b01c..6e0a6ddc2 100644 --- a/src/dynamics/thermal/initialize_temperature.cpp +++ b/src/dynamics/thermal/initialize_temperature.cpp @@ -58,7 +58,7 @@ First row is time data using std::string; using std::vector; -Temperature* InitTemperature(const std::string file_name, const double rk_prop_step_s) { +Temperature* InitTemperature(const std::string file_name, const double rk_prop_step_s, const SolarRadiationPressureEnvironment* srp_environment) { auto mainIni = IniAccess(file_name); vector node_list; @@ -153,6 +153,6 @@ Temperature* InitTemperature(const std::string file_name, const double rk_prop_s Temperature* temperature; temperature = new Temperature(conductance_matrix, radiation_matrix, node_list, heatload_list, heater_list, heater_controller_list, node_num, - rk_prop_step_s, is_calc_enabled, solar_calc_setting, debug); + rk_prop_step_s, srp_environment, is_calc_enabled, solar_calc_setting, debug); return temperature; } diff --git a/src/dynamics/thermal/initialize_temperature.hpp b/src/dynamics/thermal/initialize_temperature.hpp index 1b248a1e7..ee841f0bd 100644 --- a/src/dynamics/thermal/initialize_temperature.hpp +++ b/src/dynamics/thermal/initialize_temperature.hpp @@ -15,6 +15,6 @@ * @param[in] rk_prop_step_s: time step interval for temperature propagation integration * @return Temperature* */ -Temperature* InitTemperature(const std::string file_name, const double rk_prop_step_s); +Temperature* InitTemperature(const std::string file_name, const double rk_prop_step_s, const SolarRadiationPressureEnvironment* srp_environment); #endif // S2E_DYNAMICS_THERMAL_INITIALIZE_TEMPERATURE_HPP_ diff --git a/src/dynamics/thermal/node.cpp b/src/dynamics/thermal/node.cpp index f865bbdb2..2a225610d 100644 --- a/src/dynamics/thermal/node.cpp +++ b/src/dynamics/thermal/node.cpp @@ -29,9 +29,8 @@ Node::Node(const int node_id, const string node_name, const NodeType node_type, Node::~Node() {} -double Node::CalcSolarRadiation_W(libra::Vector<3> sun_position_b_m) { +double Node::CalcSolarRadiation_W(libra::Vector<3> sun_position_b_m, double solar_flux_W_m2) { double sun_distance_m = sun_position_b_m.CalcNorm(); - double solar_flux_W_m2 = environment::solar_constant_W_m2 * pow((environment::astronomical_unit_m / sun_distance_m), 2); double cos_theta = InnerProduct(sun_position_b_m, normal_vector_b_) / sun_distance_m / normal_vector_b_.CalcNorm(); // calculate sun_power diff --git a/src/dynamics/thermal/node.hpp b/src/dynamics/thermal/node.hpp index 84d251c1b..8324bd8c4 100644 --- a/src/dynamics/thermal/node.hpp +++ b/src/dynamics/thermal/node.hpp @@ -70,7 +70,7 @@ class Node { * @param sun_position_b_m: Sun position in body frame [m] * @return double: Solar Radiation [W] */ - double CalcSolarRadiation_W(libra::Vector<3> sun_position_b_m); + double CalcSolarRadiation_W(libra::Vector<3> sun_position_b_m, double solar_flux_W_m2); // Getter /** diff --git a/src/dynamics/thermal/temperature.cpp b/src/dynamics/thermal/temperature.cpp index 86c7cdab1..51eac3bf6 100644 --- a/src/dynamics/thermal/temperature.cpp +++ b/src/dynamics/thermal/temperature.cpp @@ -15,7 +15,8 @@ using namespace std; Temperature::Temperature(const vector> conductance_matrix_W_K, const vector> radiation_matrix_m2, vector nodes, vector heatloads, vector heaters, vector heater_controllers, const int node_num, - const double propagation_step_s, const bool is_calc_enabled, const SolarCalcSetting solar_calc_setting, const bool debug) + const double propagation_step_s, const SolarRadiationPressureEnvironment* srp_environment, const bool is_calc_enabled, + const SolarCalcSetting solar_calc_setting, const bool debug) : conductance_matrix_W_K_(conductance_matrix_W_K), radiation_matrix_m2_(radiation_matrix_m2), nodes_(nodes), @@ -23,7 +24,8 @@ Temperature::Temperature(const vector> conductance_matrix_W_K, co heaters_(heaters), heater_controllers_(heater_controllers), node_num_(node_num), - propagation_step_s_(propagation_step_s), // ルンゲクッタ積分時間刻み幅 + propagation_step_s_(propagation_step_s), + srp_environment_(srp_environment), is_calc_enabled_(is_calc_enabled), solar_calc_setting_(solar_calc_setting), debug_(debug) { @@ -121,8 +123,9 @@ vector Temperature::CalcTemperatureDifferentials(vector temperat for (int i = 0; i < node_num; i++) { heatloads_[i].SetElapsedTime_s(t); if (nodes_[i].GetNodeType() == NodeType::kDiffusive) { + double solar_flux_W_m2 = srp_environment_->GetPowerDensity_W_m2(); if (solar_calc_setting_ == SolarCalcSetting::kEnable) { - double solar_radiation_W = nodes_[i].CalcSolarRadiation_W(sun_position_b_m); + double solar_radiation_W = nodes_[i].CalcSolarRadiation_W(sun_position_b_m, solar_flux_W_m2); heatloads_[i].SetSolarHeatload_W(solar_radiation_W); } double heater_power_W = GetHeaterPower_W(i); diff --git a/src/dynamics/thermal/temperature.hpp b/src/dynamics/thermal/temperature.hpp index cb6a36133..ba001272d 100644 --- a/src/dynamics/thermal/temperature.hpp +++ b/src/dynamics/thermal/temperature.hpp @@ -6,6 +6,7 @@ #ifndef S2E_DYNAMICS_THERMAL_TEMPERATURE_HPP_ #define S2E_DYNAMICS_THERMAL_TEMPERATURE_HPP_ +#include #include #include #include @@ -38,10 +39,11 @@ class Temperature : public ILoggable { std::vector heater_controllers_; // vector of heater controllers int node_num_; // number of nodes double propagation_step_s_; // propagation step [s] - double propagation_time_s_; // Incremented time inside class Temperature [s], finish propagation when reaching end_time - bool is_calc_enabled_; // Whether temperature calculation is enabled - SolarCalcSetting solar_calc_setting_; // setting for solar calculation - bool debug_; // Activate debug output or not + double propagation_time_s_; // Incremented time inside class Temperature [s], finish propagation when reaching end_time + const SolarRadiationPressureEnvironment* srp_environment_; // SolarRadiationPressureEnvironment for calculating solar flux + bool is_calc_enabled_; // Whether temperature calculation is enabled + SolarCalcSetting solar_calc_setting_; // setting for solar calculation + bool debug_; // Activate debug output or not /** * @fn CalcRungeOneStep @@ -85,8 +87,8 @@ class Temperature : public ILoggable { */ Temperature(const std::vector> conductance_matrix_W_K, const std::vector> radiation_matrix_m2, std::vector nodes, std::vector heatloads, std::vector heaters, std::vector heater_controllers, - const int node_num, const double propagation_step_s, const bool is_calc_enabled, const SolarCalcSetting solar_calc_setting, - const bool debug); + const int node_num, const double propagation_step_s, const SolarRadiationPressureEnvironment* srp_environment, const bool is_calc_enabled, + const SolarCalcSetting solar_calc_setting, const bool debug); /** * @fn Temperature * @brief Construct a new Temperature object, used when thermal calculation is disabled. diff --git a/src/environment/local/solar_radiation_pressure_environment.cpp b/src/environment/local/solar_radiation_pressure_environment.cpp index d056e7292..d7b204d93 100644 --- a/src/environment/local/solar_radiation_pressure_environment.cpp +++ b/src/environment/local/solar_radiation_pressure_environment.cpp @@ -19,6 +19,8 @@ SolarRadiationPressureEnvironment::SolarRadiationPressureEnvironment(LocalCelest sun_radius_m_ = local_celestial_information_->GetGlobalInformation().GetMeanRadiusFromName_m("SUN"); } +SolarRadiationPressureEnvironment::SolarRadiationPressureEnvironment() {} + void SolarRadiationPressureEnvironment::UpdateAllStates() { if (!IsCalcEnabled) return; diff --git a/src/environment/local/solar_radiation_pressure_environment.hpp b/src/environment/local/solar_radiation_pressure_environment.hpp index ecb3c5c29..086683629 100644 --- a/src/environment/local/solar_radiation_pressure_environment.hpp +++ b/src/environment/local/solar_radiation_pressure_environment.hpp @@ -23,6 +23,13 @@ class SolarRadiationPressureEnvironment : public ILoggable { * @param [in] local_celestial_information: Local celestial information */ SolarRadiationPressureEnvironment(LocalCelestialInformation* local_celestial_information); + + /** + * @fn SolarRadiationPressureEnvironment + * @brief Constructor (Use for constructing Temperature with no params) + */ + SolarRadiationPressureEnvironment(); + /** * @fn ~SolarRadiationPressureEnvironment * @brief Destructor @@ -80,11 +87,11 @@ class SolarRadiationPressureEnvironment : public ILoggable { virtual std::string GetLogValue() const; private: - double solar_radiation_pressure_N_m2_; //!< Solar radiation pressure [N/m^2] - double solar_constant_W_m2_ = 1366.0; //!< Solar constant [W/m^2] TODO: We need to change the value depends on sun activity. - double shadow_coefficient_ = 1.0; //!< Shadow function - double sun_radius_m_; //!< Sun radius [m] - std::string shadow_source_name_; //!< Shadow source name + double solar_radiation_pressure_N_m2_; //!< Solar radiation pressure [N/m^2] + double solar_constant_W_m2_ = 1366.0; //!< Solar constant [W/m^2] TODO: We need to change the value depends on sun activity. + double shadow_coefficient_ = 1.0; //!< Shadow function + double sun_radius_m_; //!< Sun radius [m] + std::string shadow_source_name_; //!< Shadow source name LocalCelestialInformation* local_celestial_information_; //!< Local celestial information From f6a50e49bde39f3931e5b33a6ea0662fd7029d14 Mon Sep 17 00:00:00 2001 From: tshibuk Date: Tue, 6 Jun 2023 14:14:01 +0900 Subject: [PATCH 3/8] fix format --- src/dynamics/thermal/temperature.hpp | 12 ++++++------ .../local/solar_radiation_pressure_environment.cpp | 4 ++-- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/dynamics/thermal/temperature.hpp b/src/dynamics/thermal/temperature.hpp index ba001272d..2e71d1618 100644 --- a/src/dynamics/thermal/temperature.hpp +++ b/src/dynamics/thermal/temperature.hpp @@ -39,11 +39,11 @@ class Temperature : public ILoggable { std::vector heater_controllers_; // vector of heater controllers int node_num_; // number of nodes double propagation_step_s_; // propagation step [s] - double propagation_time_s_; // Incremented time inside class Temperature [s], finish propagation when reaching end_time + double propagation_time_s_; // Incremented time inside class Temperature [s], finish propagation when reaching end_time const SolarRadiationPressureEnvironment* srp_environment_; // SolarRadiationPressureEnvironment for calculating solar flux - bool is_calc_enabled_; // Whether temperature calculation is enabled - SolarCalcSetting solar_calc_setting_; // setting for solar calculation - bool debug_; // Activate debug output or not + bool is_calc_enabled_; // Whether temperature calculation is enabled + SolarCalcSetting solar_calc_setting_; // setting for solar calculation + bool debug_; // Activate debug output or not /** * @fn CalcRungeOneStep @@ -87,8 +87,8 @@ class Temperature : public ILoggable { */ Temperature(const std::vector> conductance_matrix_W_K, const std::vector> radiation_matrix_m2, std::vector nodes, std::vector heatloads, std::vector heaters, std::vector heater_controllers, - const int node_num, const double propagation_step_s, const SolarRadiationPressureEnvironment* srp_environment, const bool is_calc_enabled, - const SolarCalcSetting solar_calc_setting, const bool debug); + const int node_num, const double propagation_step_s, const SolarRadiationPressureEnvironment* srp_environment, + const bool is_calc_enabled, const SolarCalcSetting solar_calc_setting, const bool debug); /** * @fn Temperature * @brief Construct a new Temperature object, used when thermal calculation is disabled. diff --git a/src/environment/local/solar_radiation_pressure_environment.cpp b/src/environment/local/solar_radiation_pressure_environment.cpp index d7b204d93..4a673a7da 100644 --- a/src/environment/local/solar_radiation_pressure_environment.cpp +++ b/src/environment/local/solar_radiation_pressure_environment.cpp @@ -86,11 +86,11 @@ void SolarRadiationPressureEnvironment::CalcShadowCoefficient(std::string shadow } else if (c < fabs(a - b) && a > b) // The occultation is partial but maximum { shadow_coefficient_ = 1.0 - (b * b) / (a * a); - } else if (fabs(a - b) <= c && c <= (a + b)) // spacecraft is in penumbra + } else if (fabs(a - b) <= c && c <= (a + b)) // spacecraft is in penumbra { double A = a * a * acos(x / a) + b * b * acos((c - x) / b) - c * y; // The area of the occulted segment of the apparent solar disk shadow_coefficient_ = 1.0 - A / (libra::pi * a * a); - } else { // no occultation takes place + } else { // no occultation takes place assert(c > (a + b)); shadow_coefficient_ = 1.0; } From 0d9014360485cda9bb98c37e8c551284a6044bb4 Mon Sep 17 00:00:00 2001 From: tshibuk Date: Fri, 9 Jun 2023 11:19:27 +0900 Subject: [PATCH 4/8] fix format --- .../local/solar_radiation_pressure_environment.cpp | 4 ++-- .../local/solar_radiation_pressure_environment.hpp | 10 +++++----- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/environment/local/solar_radiation_pressure_environment.cpp b/src/environment/local/solar_radiation_pressure_environment.cpp index 4a673a7da..d7b204d93 100644 --- a/src/environment/local/solar_radiation_pressure_environment.cpp +++ b/src/environment/local/solar_radiation_pressure_environment.cpp @@ -86,11 +86,11 @@ void SolarRadiationPressureEnvironment::CalcShadowCoefficient(std::string shadow } else if (c < fabs(a - b) && a > b) // The occultation is partial but maximum { shadow_coefficient_ = 1.0 - (b * b) / (a * a); - } else if (fabs(a - b) <= c && c <= (a + b)) // spacecraft is in penumbra + } else if (fabs(a - b) <= c && c <= (a + b)) // spacecraft is in penumbra { double A = a * a * acos(x / a) + b * b * acos((c - x) / b) - c * y; // The area of the occulted segment of the apparent solar disk shadow_coefficient_ = 1.0 - A / (libra::pi * a * a); - } else { // no occultation takes place + } else { // no occultation takes place assert(c > (a + b)); shadow_coefficient_ = 1.0; } diff --git a/src/environment/local/solar_radiation_pressure_environment.hpp b/src/environment/local/solar_radiation_pressure_environment.hpp index 086683629..77d31e3e1 100644 --- a/src/environment/local/solar_radiation_pressure_environment.hpp +++ b/src/environment/local/solar_radiation_pressure_environment.hpp @@ -87,11 +87,11 @@ class SolarRadiationPressureEnvironment : public ILoggable { virtual std::string GetLogValue() const; private: - double solar_radiation_pressure_N_m2_; //!< Solar radiation pressure [N/m^2] - double solar_constant_W_m2_ = 1366.0; //!< Solar constant [W/m^2] TODO: We need to change the value depends on sun activity. - double shadow_coefficient_ = 1.0; //!< Shadow function - double sun_radius_m_; //!< Sun radius [m] - std::string shadow_source_name_; //!< Shadow source name + double solar_radiation_pressure_N_m2_; //!< Solar radiation pressure [N/m^2] + double solar_constant_W_m2_ = 1366.0; //!< Solar constant [W/m^2] TODO: We need to change the value depends on sun activity. + double shadow_coefficient_ = 1.0; //!< Shadow function + double sun_radius_m_; //!< Sun radius [m] + std::string shadow_source_name_; //!< Shadow source name LocalCelestialInformation* local_celestial_information_; //!< Local celestial information From 6895d9a8cda61d83d49ea9097b3a3425b7a6fd5b Mon Sep 17 00:00:00 2001 From: tshibuk Date: Fri, 9 Jun 2023 11:26:34 +0900 Subject: [PATCH 5/8] fix format --- src/dynamics/thermal/initialize_heatload.cpp | 8 ++++---- src/dynamics/thermal/initialize_node.cpp | 16 ++++++++-------- 2 files changed, 12 insertions(+), 12 deletions(-) diff --git a/src/dynamics/thermal/initialize_heatload.cpp b/src/dynamics/thermal/initialize_heatload.cpp index 888353593..a0e16176e 100644 --- a/src/dynamics/thermal/initialize_heatload.cpp +++ b/src/dynamics/thermal/initialize_heatload.cpp @@ -31,14 +31,14 @@ Heatload InitHeatload(const std::vector& time_str, const std::vecto using std::stod; using std::stoi; - assert(time_str.size() >= 2); // Number of columns must be larger than 2 - assert(internal_heatload_str.size() >= 2); // Number of columns must be larger than 2 - assert(time_str.size() == internal_heatload_str.size()); // Number of columns must be same + assert(time_str.size() >= 2); // Number of columns must be larger than 2 + assert(internal_heatload_str.size() >= 2); // Number of columns must be larger than 2 + assert(time_str.size() == internal_heatload_str.size()); // Number of columns must be same std::vector time_table_s(time_str.size() - 1); // exclude index std::vector internal_heatload_table_W(internal_heatload_str.size() - 1); // exclude index - int node_id = stoi(internal_heatload_str[0]); // First data of internal_heatload_str is node id + int node_id = stoi(internal_heatload_str[0]); // First data of internal_heatload_str is node id // read table for (unsigned int i = 0; i < time_str.size() - 1; ++i) { diff --git a/src/dynamics/thermal/initialize_node.cpp b/src/dynamics/thermal/initialize_node.cpp index fce137243..c58b4a330 100644 --- a/src/dynamics/thermal/initialize_node.cpp +++ b/src/dynamics/thermal/initialize_node.cpp @@ -42,14 +42,14 @@ Node InitNode(const std::vector& node_str) { size_t node_str_size_defined = 11; // Correct size of node_str assert(node_str.size() == node_str_size_defined); // Check if size of node_str is correct - int node_id = 0; // node number - std::string node_label = "temp"; // node name - int node_type_int = 0; // node type - int heater_id = 0; // heater node index - double temperature_K = 0; // [K] - double capacity_J_K = 0; // [J/K] - double alpha = 0; // [] - double area_m2 = 0; // [m^2] + int node_id = 0; // node number + std::string node_label = "temp"; // node name + int node_type_int = 0; // node type + int heater_id = 0; // heater node index + double temperature_K = 0; // [K] + double capacity_J_K = 0; // [J/K] + double alpha = 0; // [] + double area_m2 = 0; // [m^2] // Index to read from node_str for each parameter int index_node_id = 0; From 85ab4b1ddea99b2cd4ceb03a93203125565b3905 Mon Sep 17 00:00:00 2001 From: tshibuk Date: Wed, 27 Sep 2023 10:38:24 +0900 Subject: [PATCH 6/8] resolve conflict in initialize_temperature.cpp --- src/dynamics/thermal/initialize_temperature.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/dynamics/thermal/initialize_temperature.cpp b/src/dynamics/thermal/initialize_temperature.cpp index 6e0a6ddc2..ca68bb419 100644 --- a/src/dynamics/thermal/initialize_temperature.cpp +++ b/src/dynamics/thermal/initialize_temperature.cpp @@ -70,9 +70,9 @@ Temperature* InitTemperature(const std::string file_name, const double rk_prop_s vector> node_str_list; // string vector of node property data vector> heater_str_list; // string vector of heater property data vector> heatload_str_list; // string vector of heatload property data - unsigned int node_num = 1; - unsigned int heatload_num = 1; - unsigned int heater_num = 1; + size_t node_num = 1; + size_t heatload_num = 1; + size_t heater_num = 1; bool is_calc_enabled = mainIni.ReadEnable("THERMAL", "calculation"); if (is_calc_enabled == false) { From 7d43f1180242fb49f84f58fbf0f3f006bc287bf3 Mon Sep 17 00:00:00 2001 From: tshibuk Date: Wed, 27 Sep 2023 11:28:57 +0900 Subject: [PATCH 7/8] modify CalcSolarRadiation_W based on comments --- src/dynamics/thermal/initialize_node.cpp | 9 ++++++++ src/dynamics/thermal/node.cpp | 5 ++--- src/dynamics/thermal/node.hpp | 4 ++-- src/dynamics/thermal/temperature.cpp | 26 ++++++++++++++---------- src/dynamics/thermal/temperature.hpp | 8 ++++---- 5 files changed, 32 insertions(+), 20 deletions(-) diff --git a/src/dynamics/thermal/initialize_node.cpp b/src/dynamics/thermal/initialize_node.cpp index c58b4a330..8624797ea 100644 --- a/src/dynamics/thermal/initialize_node.cpp +++ b/src/dynamics/thermal/initialize_node.cpp @@ -73,6 +73,15 @@ Node InitNode(const std::vector& node_str) { for (int i = 0; i < 3; i++) { normal_v_b[i] = stod(node_str[index_normal_v_b_head + i]); } + + // Normalize Norm Vector (Except for Boundary and Arithmetic Nodes) + if (node_type_int == 0) { + double norm = normal_v_b.CalcNorm(); + for (int i = 0; i < 3; i++) { + normal_v_b[i] = normal_v_b[i] / norm; + } + } + temperature_K = stod(node_str[index_temperature]); NodeType node_type = NodeType::kDiffusive; diff --git a/src/dynamics/thermal/node.cpp b/src/dynamics/thermal/node.cpp index 2a225610d..586bd6c51 100644 --- a/src/dynamics/thermal/node.cpp +++ b/src/dynamics/thermal/node.cpp @@ -29,9 +29,8 @@ Node::Node(const int node_id, const string node_name, const NodeType node_type, Node::~Node() {} -double Node::CalcSolarRadiation_W(libra::Vector<3> sun_position_b_m, double solar_flux_W_m2) { - double sun_distance_m = sun_position_b_m.CalcNorm(); - double cos_theta = InnerProduct(sun_position_b_m, normal_vector_b_) / sun_distance_m / normal_vector_b_.CalcNorm(); +double Node::CalcSolarRadiation_W(libra::Vector<3> sun_direction_b, double solar_flux_W_m2) { + double cos_theta = InnerProduct(sun_direction_b, normal_vector_b_); // calculate sun_power if (cos_theta > 0) diff --git a/src/dynamics/thermal/node.hpp b/src/dynamics/thermal/node.hpp index 8324bd8c4..6999b09ce 100644 --- a/src/dynamics/thermal/node.hpp +++ b/src/dynamics/thermal/node.hpp @@ -67,10 +67,10 @@ class Node { * @fn CalcSolarRadiation_W * @brief Calculate solar radiation [W] from sun direction, alpha, area, and normal vector * - * @param sun_position_b_m: Sun position in body frame [m] + * @param sun_direction_b: Sun direction in body frame * @return double: Solar Radiation [W] */ - double CalcSolarRadiation_W(libra::Vector<3> sun_position_b_m, double solar_flux_W_m2); + double CalcSolarRadiation_W(libra::Vector<3> sun_direction_b, double solar_flux_W_m2); // Getter /** diff --git a/src/dynamics/thermal/temperature.cpp b/src/dynamics/thermal/temperature.cpp index 51eac3bf6..1c31f7398 100644 --- a/src/dynamics/thermal/temperature.cpp +++ b/src/dynamics/thermal/temperature.cpp @@ -48,11 +48,16 @@ Temperature::~Temperature() {} void Temperature::Propagate(libra::Vector<3> sun_position_b_m, const double time_end_s) { if (!is_calc_enabled_) return; + double sun_distance_m = sun_position_b_m.CalcNorm(); + libra::Vector<3> sun_direction_b; + for (int i = 0; i < 3; i++) { + sun_direction_b[i] = sun_position_b_m[i] / sun_distance_m; + } while (time_end_s - propagation_time_s_ - propagation_step_s_ > 1.0e-6) { - CalcRungeOneStep(propagation_time_s_, propagation_step_s_, sun_position_b_m, node_num_); + CalcRungeOneStep(propagation_time_s_, propagation_step_s_, sun_direction_b, node_num_); propagation_time_s_ += propagation_step_s_; } - CalcRungeOneStep(propagation_time_s_, time_end_s - propagation_time_s_, sun_position_b_m, node_num_); + CalcRungeOneStep(propagation_time_s_, time_end_s - propagation_time_s_, sun_direction_b, node_num_); propagation_time_s_ = time_end_s; UpdateHeaterStatus(); @@ -67,9 +72,8 @@ void Temperature::Propagate(libra::Vector<3> sun_position_b_m, const double time cout << setprecision(4) << itr->GetSolarRadiation_W() << " "; } cout << "SunDir: "; - double sun_distance_m = sun_position_b_m.CalcNorm(); for (int i = 0; i < 3; i++) { - cout << setprecision(3) << sun_position_b_m[i] / sun_distance_m << " "; + cout << setprecision(3) << sun_direction_b[i] << " "; } cout << "Heatload: "; for (auto itr = heatloads_.begin(); itr != heatloads_.end(); ++itr) { @@ -79,7 +83,7 @@ void Temperature::Propagate(libra::Vector<3> sun_position_b_m, const double time } } -void Temperature::CalcRungeOneStep(double time_now_s, double time_step_s, libra::Vector<3> sun_position_b_m, int node_num) { +void Temperature::CalcRungeOneStep(double time_now_s, double time_step_s, libra::Vector<3> sun_direction_b, int node_num) { vector temperatures_now_K(node_num); for (int i = 0; i < node_num; i++) { temperatures_now_K[i] = nodes_[i].GetTemperature_K(); @@ -88,22 +92,22 @@ void Temperature::CalcRungeOneStep(double time_now_s, double time_step_s, libra: vector k1(node_num), k2(node_num), k3(node_num), k4(node_num); vector xk2(node_num), xk3(node_num), xk4(node_num); - k1 = CalcTemperatureDifferentials(temperatures_now_K, time_now_s, sun_position_b_m, node_num); + k1 = CalcTemperatureDifferentials(temperatures_now_K, time_now_s, sun_direction_b, node_num); for (int i = 0; i < node_num; i++) { xk2[i] = temperatures_now_K[i] + (time_step_s / 2.0) * k1[i]; } - k2 = CalcTemperatureDifferentials(xk2, (time_now_s + time_step_s / 2.0), sun_position_b_m, node_num); + k2 = CalcTemperatureDifferentials(xk2, (time_now_s + time_step_s / 2.0), sun_direction_b, node_num); for (int i = 0; i < node_num; i++) { xk3[i] = temperatures_now_K[i] + (time_step_s / 2.0) * k2[i]; } - k3 = CalcTemperatureDifferentials(xk3, (time_now_s + time_step_s / 2.0), sun_position_b_m, node_num); + k3 = CalcTemperatureDifferentials(xk3, (time_now_s + time_step_s / 2.0), sun_direction_b, node_num); for (int i = 0; i < node_num; i++) { xk4[i] = temperatures_now_K[i] + time_step_s * k3[i]; } - k4 = CalcTemperatureDifferentials(xk4, (time_now_s + time_step_s), sun_position_b_m, node_num); + k4 = CalcTemperatureDifferentials(xk4, (time_now_s + time_step_s), sun_direction_b, node_num); vector temperatures_next_K(node_num); for (int i = 0; i < node_num; i++) { @@ -115,7 +119,7 @@ void Temperature::CalcRungeOneStep(double time_now_s, double time_step_s, libra: } } -vector Temperature::CalcTemperatureDifferentials(vector temperatures_K, double t, libra::Vector<3> sun_position_b_m, int node_num) { +vector Temperature::CalcTemperatureDifferentials(vector temperatures_K, double t, libra::Vector<3> sun_direction_b, int node_num) { // TODO: consider the following unused arguments are really needed UNUSED(temperatures_K); @@ -125,7 +129,7 @@ vector Temperature::CalcTemperatureDifferentials(vector temperat if (nodes_[i].GetNodeType() == NodeType::kDiffusive) { double solar_flux_W_m2 = srp_environment_->GetPowerDensity_W_m2(); if (solar_calc_setting_ == SolarCalcSetting::kEnable) { - double solar_radiation_W = nodes_[i].CalcSolarRadiation_W(sun_position_b_m, solar_flux_W_m2); + double solar_radiation_W = nodes_[i].CalcSolarRadiation_W(sun_direction_b, solar_flux_W_m2); heatloads_[i].SetSolarHeatload_W(solar_radiation_W); } double heater_power_W = GetHeaterPower_W(i); diff --git a/src/dynamics/thermal/temperature.hpp b/src/dynamics/thermal/temperature.hpp index 2e71d1618..ca0af4265 100644 --- a/src/dynamics/thermal/temperature.hpp +++ b/src/dynamics/thermal/temperature.hpp @@ -51,21 +51,21 @@ class Temperature : public ILoggable { * * @param[in] time_now_s: Current elapsed time [s] * @param[in] time_step_s: Time step of RK4 [s] - * @param[in] sun_position_b_m: Sun position in body frame [m] + * @param[in] sun_direction_b: Sun position in body frame [m] * @param[in] node_num: Number of nodes */ - void CalcRungeOneStep(double time_now_s, double time_step_s, libra::Vector<3> sun_position_b_m, int node_num); + void CalcRungeOneStep(double time_now_s, double time_step_s, libra::Vector<3> sun_direction_b, int node_num); /** * @fn CalcTemperatureDifferentials * @brief Calculate differential of thermal equilibrium equation * * @param temperatures_K: [UNUSED] Temperatures of each node [K] * @param time_now_s: Current elapsed time [s] - * @param[in] sun_position_b_m: Sun position in body frame [m] + * @param[in] sun_direction_b: Sun direction in body frame * @param node_num: Number of nodes * @return std::vector: Differential of thermal equilibrium equation at time now */ - std::vector CalcTemperatureDifferentials(std::vector temperatures_K, double time_now_s, const libra::Vector<3> sun_position_b_m, + std::vector CalcTemperatureDifferentials(std::vector temperatures_K, double time_now_s, const libra::Vector<3> sun_direction_b, int node_num); public: From c6cad97eb53de7e3d9159a67cfe626af084af39f Mon Sep 17 00:00:00 2001 From: tshibuk Date: Wed, 27 Sep 2023 11:32:23 +0900 Subject: [PATCH 8/8] small fix --- src/dynamics/thermal/initialize_node.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/dynamics/thermal/initialize_node.cpp b/src/dynamics/thermal/initialize_node.cpp index 8624797ea..25b51775c 100644 --- a/src/dynamics/thermal/initialize_node.cpp +++ b/src/dynamics/thermal/initialize_node.cpp @@ -81,7 +81,7 @@ Node InitNode(const std::vector& node_str) { normal_v_b[i] = normal_v_b[i] / norm; } } - + temperature_K = stod(node_str[index_temperature]); NodeType node_type = NodeType::kDiffusive;