diff --git a/applications/FluidDynamicsApplication/custom_processes/compute_aerodynamic_coefficients_process.cpp b/applications/FluidDynamicsApplication/custom_processes/compute_aerodynamic_coefficients_process.cpp new file mode 100644 index 000000000000..08ddc31b73a7 --- /dev/null +++ b/applications/FluidDynamicsApplication/custom_processes/compute_aerodynamic_coefficients_process.cpp @@ -0,0 +1,195 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main authors: Juan I. Camarotti +// +// + +// System includes + +// External includes + +// Project includes + +// Application includes +#include "compute_aerodynamic_coefficients_process.h" +#include "fluid_dynamics_application_variables.h" +#include "utilities/variable_utils.h" + + +namespace Kratos +{ + +ComputeAerodynamicCoefficientsProcess::ComputeAerodynamicCoefficientsProcess( + Model &rModel, + Parameters Params) + : Process(), + mrModelPart(rModel.GetModelPart(Params["model_part_name"].GetString())) +{ + // Check default settings + KRATOS_TRY + + Params.ValidateAndAssignDefaults(GetDefaultParameters()); + + ReadProcessParameters(Params); + + KRATOS_CATCH("") +} + + +const Parameters ComputeAerodynamicCoefficientsProcess::GetDefaultParameters() const +{ + return Parameters( R"( + { + "model_part_name" : "PLEASE_PROVIDE_A_MODELPART_NAME", + "reference_surface" : 0.0, + "reference_chord" : 1.0, + "reference_span" : 1.0, + "moment_reference_point" : [0.0, 0.0, 0.0], + "freestream_dynamic_pressure" : 1.0, + "angle_of_attack" : 0.0, + "sideslip_angle" : 0.0 + })" ); +} + +void ComputeAerodynamicCoefficientsProcess::ReadProcessParameters(const Parameters& rParams) +{ + constexpr double tol = 1e-12; + + mReference_Surface = rParams["reference_surface"].GetDouble(); + KRATOS_ERROR_IF(std::abs(mReference_Surface) <= tol) + << "Invalid value for 'reference_surface' = " << mReference_Surface << "."; + + mReference_Chord = rParams["reference_chord"].GetDouble(); + KRATOS_ERROR_IF(std::abs(mReference_Chord) <= tol) + << "Invalid value for 'reference_chord' = " << mReference_Chord << "."; + + mReference_Span = rParams["reference_span"].GetDouble(); + KRATOS_ERROR_IF(std::abs(mReference_Span) <= tol) + << "Invalid value for 'reference_span' = " << mReference_Span << "."; + + mQInf = rParams["freestream_dynamic_pressure"].GetDouble(); + KRATOS_ERROR_IF(std::abs(mQInf) <= tol) + << "Invalid value for 'freestream_dynamic_pressure' = " << mQInf << "."; + + mAngleOfAttack = rParams["angle_of_attack"].GetDouble(); + mSideslipAngle = rParams["sideslip_angle"].GetDouble(); + + const auto& r_mrp = rParams["moment_reference_point"]; + KRATOS_ERROR_IF_NOT(r_mrp.IsArray() && r_mrp.size() == 3) + << "'moment_reference_point' must be an array of size 3."; + mMomentReferencePoint[0] = r_mrp[0].GetDouble(); + mMomentReferencePoint[1] = r_mrp[1].GetDouble(); + mMomentReferencePoint[2] = r_mrp[2].GetDouble(); +} + + +void ComputeAerodynamicCoefficientsProcess::ExecuteBeforeOutputStep() +{ + Execute(); +} + + +void ComputeAerodynamicCoefficientsProcess::Execute() +{ + KRATOS_TRY; + + // Build wind axes from alpha (AoA) and beta (sideslip) + const double alpha = mAngleOfAttack * Globals::Pi / 180.0; + const double beta = mSideslipAngle * Globals::Pi / 180.0; + + // Wind-x (drag axis): freestream direction (unit) + array_1d e_drag = ZeroVector(3); + e_drag[0] = std::cos(alpha) * std::cos(beta); + e_drag[1] = std::sin(beta); + e_drag[2] = std::sin(alpha) * std::cos(beta); + + const double n_drag = norm_2(e_drag); + KRATOS_ERROR_IF(n_drag < 1e-15) << "Invalid wind axis: e_drag near zero."; + e_drag /= n_drag; + + // Build a reference "up" vector to construct orthonormal basis + array_1d up = ZeroVector(3); + up[0] = 0.0; up[1] = 0.0; up[2] = 1.0; + if (std::abs(inner_prod(e_drag, up)) > 0.99) { + up[0] = 0.0; up[1] = 1.0; up[2] = 0.0; + } + + // Wind-y (side axis) + array_1d e_side = MathUtils::CrossProduct(up, e_drag); + const double n_side = norm_2(e_side); + KRATOS_ERROR_IF(n_side < 1e-15) << "Invalid wind axis: e_side near zero."; + e_side /= n_side; + + // Wind-z (lift axis) + array_1d e_lift = MathUtils::CrossProduct(e_drag, e_side); + const double n_lift = norm_2(e_lift); + KRATOS_ERROR_IF(n_lift < 1e-15) << "Invalid wind axis: e_lift near zero."; + e_lift /= n_lift; + + // Accumulate aerodynamic force and moment from nodal reactions + array_1d F_aero = ZeroVector(3); + array_1d M_aero = ZeroVector(3); + + for (auto& r_node : mrModelPart.Nodes()) + { + KRATOS_ERROR_IF_NOT(r_node.SolutionStepsDataHas(REACTION)) + << "Node " << r_node.Id() << " does not have REACTION in SolutionStepData."; + + const array_1d& r_R = r_node.FastGetSolutionStepValue(REACTION); + const array_1d f_node = -r_R; + + F_aero += f_node; + + array_1d r = r_node.Coordinates() - mMomentReferencePoint; + M_aero += MathUtils::CrossProduct(r, f_node); + } + + // Projections in wind axes + const double D = inner_prod(F_aero, e_drag); // drag + const double Y = inner_prod(F_aero, e_side); // side force + const double L = inner_prod(F_aero, e_lift); // lift + + const double Mx = inner_prod(M_aero, e_drag); // moment about wind-x (roll in wind frame) + const double My = inner_prod(M_aero, e_side); // moment about wind-y (pitch in wind frame) + const double Mz = inner_prod(M_aero, e_lift); // moment about wind-z (yaw in wind frame) + + // Denominator for non-dimensional coefficients + const double denom_F = mQInf * mReference_Surface; + + KRATOS_ERROR_IF(std::abs(denom_F) < 1e-15) << "Invalid denominator q_inf*S_ref."; + + const double C_D = D / denom_F; + const double C_Y = Y / denom_F; + const double C_L = L / denom_F; + + const double C_l = Mx / (denom_F * mReference_Span); + const double C_m = My / (denom_F * mReference_Chord); + const double C_n = Mz / (denom_F * mReference_Span); + + // Store on ModelPart + mrModelPart.SetValue(LIFT_COEFFICIENT, C_L); + mrModelPart.SetValue(DRAG_COEFFICIENT, C_D); + mrModelPart.SetValue(LATERAL_FORCE_COEFFICIENT, C_Y); + + mrModelPart.SetValue(ROLLING_MOMENT_COEFFICIENT, C_l); + mrModelPart.SetValue(PITCHING_MOMENT_COEFFICIENT, C_m); + mrModelPart.SetValue(YAWING_MOMENT_COEFFICIENT, C_n); + + std::cout << "Aero coeffs from nodal reactions | " + << "alpha=" << mAngleOfAttack << " deg, beta=" << mSideslipAngle << " deg : " + << "CL=" << C_L << ", CD=" << C_D << ", CY=" << C_Y + << ", Cl=" << C_l << ", Cm=" << C_m << ", Cn=" << C_n + << std::endl; + + KRATOS_CATCH(""); +} + + +} \ No newline at end of file diff --git a/applications/FluidDynamicsApplication/custom_processes/compute_aerodynamic_coefficients_process.h b/applications/FluidDynamicsApplication/custom_processes/compute_aerodynamic_coefficients_process.h new file mode 100644 index 000000000000..67a6674cb1ff --- /dev/null +++ b/applications/FluidDynamicsApplication/custom_processes/compute_aerodynamic_coefficients_process.h @@ -0,0 +1,193 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main authors: Juan I. Camarotti +// +// + +#ifndef KRATOS_COMPUTE_AERODYNAMIC_COEFFICIENTS_PROCESS_H +#define KRATOS_COMPUTE_AERODYNAMIC_COEFFICIENTS_PROCESS_H + +// System includes + +// External includes + +// Project includes +#include "includes/model_part.h" +#include "includes/define.h" +#include "processes/process.h" +#include "includes/kratos_parameters.h" + +// Application includes + + +namespace Kratos +{ +///@addtogroup FluidDynamicsApplication +///@{ + +///@name Kratos Globals +///@{ + +///@} +///@name Type Definitions +///@{ + + +///@} +///@name Enum's +///@{ + +///@} +///@name Functions +///@{ + +///@} +///@name Kratos Classes +///@{ + +/// This process computes the lift coefficient as a function of reference fluid properties +class KRATOS_API(FLUID_DYNAMICS_APPLICATION) ComputeAerodynamicCoefficientsProcess : public Process +{ +public: + ///@name Type Definitions + ///@{ + using NodeType = ModelPart::NodeType; + + /// Pointer definition of ComputeAerodynamicCoefficientsProcess + KRATOS_CLASS_POINTER_DEFINITION(ComputeAerodynamicCoefficientsProcess); + + ///@} + ///@name Life Cycle + ///@{ + + /// Constructor with Kratos model + ComputeAerodynamicCoefficientsProcess( + Model& rModel, + Parameters Params); + + /// Destructor. + ~ComputeAerodynamicCoefficientsProcess() override {} + + ///@} + ///@name Operators + ///@{ + + ///@} + ///@name Operations + ///@{ + + const Parameters GetDefaultParameters() const override; + + void Execute() override; + + void ExecuteBeforeOutputStep() override; + + ///@} + ///@name Access + ///@{ + + ///@} + ///@name Inquiry + ///@{ + + ///@} + ///@name Input and output + ///@{ + + /// Turn back information as a string. + std::string Info() const override + { + std::stringstream buffer; + buffer << "ComputeAerodynamicCoefficientsProcess" ; + return buffer.str(); + } + + /// Print information about this object. + void PrintInfo(std::ostream& rOStream) const override {rOStream << "ComputeAerodynamicCoefficientsProcess";} + + /// Print object's data. + void PrintData(std::ostream& rOStream) const override {} + + + ///@} + ///@name Friends + ///@{ + + ///@} + +private: + ///@name Static Member Variables + ///@{ + + ///@} + ///@name Member Variables + ///@{ + + ModelPart& mrModelPart; + double mReference_Surface = 0.0; + double mReference_Chord = 1.0; + double mReference_Span = 1.0; + double mQInf = 1.0; + array_1d mMomentReferencePoint = ZeroVector(3); + double mAngleOfAttack = 0.0; + double mSideslipAngle = 0.0; + + ///@} + ///@name Private Operators + ///@{ + + ///@} + ///@name Private Operations + ///@{ + + void ReadProcessParameters(const Parameters& Params); + + ///@} + ///@name Private Access + ///@{ + + + ///@} + ///@name Private Inquiry + ///@{ + + + ///@} + ///@name Un accessible methods + ///@{ + + /// Default constructor. + ComputeAerodynamicCoefficientsProcess() = delete; + + /// Assignment operator. + ComputeAerodynamicCoefficientsProcess& operator=(ComputeAerodynamicCoefficientsProcess const& rOther) = delete; + + /// Copy constructor. + ComputeAerodynamicCoefficientsProcess(ComputeAerodynamicCoefficientsProcess const& rOther) = delete; + + ///@} + +}; // Class ComputeAerodynamicCoefficientsProcess + +///@} +///@name Type Definitions +///@{ + +///@} +///@name Input and output +///@{ + +///@} + +///@} addtogroup block + +}; // namespace Kratos. + +#endif // KRATOS_COMPUTE_AERODYNAMIC_COEFFICIENTS_PROCESS_H diff --git a/applications/FluidDynamicsApplication/custom_python/add_custom_processes_to_python.cpp b/applications/FluidDynamicsApplication/custom_python/add_custom_processes_to_python.cpp index e85de405bee5..c1b49b9848ca 100644 --- a/applications/FluidDynamicsApplication/custom_python/add_custom_processes_to_python.cpp +++ b/applications/FluidDynamicsApplication/custom_python/add_custom_processes_to_python.cpp @@ -28,6 +28,7 @@ #include "custom_processes/boussinesq_force_process.h" #include "custom_processes/calculate_levelset_consistent_nodal_gradient_process.h" #include "custom_processes/compute_pressure_coefficient_process.h" +#include "custom_processes/compute_aerodynamic_coefficients_process.h" #include "custom_processes/distance_modification_process.h" #include "custom_processes/distance_smoothing_process.h" #include "custom_processes/embedded_nodes_initialization_process.h" @@ -181,6 +182,10 @@ void AddCustomProcessesToPython(pybind11::module& m) .def(py::init()) ; + py::class_(m, "ComputeAerodynamicCoefficientsProcess") + .def(py::init()) + ; + py::class_(m, "ComputeYPlusProcess") .def(py::init()) ; diff --git a/applications/FluidDynamicsApplication/custom_python/fluid_dynamics_python_application.cpp b/applications/FluidDynamicsApplication/custom_python/fluid_dynamics_python_application.cpp index 81d1d7c493e0..bdca9db22258 100644 --- a/applications/FluidDynamicsApplication/custom_python/fluid_dynamics_python_application.cpp +++ b/applications/FluidDynamicsApplication/custom_python/fluid_dynamics_python_application.cpp @@ -148,6 +148,14 @@ PYBIND11_MODULE(KratosFluidDynamicsApplication,m) KRATOS_REGISTER_IN_PYTHON_VARIABLE(m,SPECTRAL_RADIUS_LIMIT) KRATOS_REGISTER_IN_PYTHON_3D_VARIABLE_WITH_COMPONENTS(m, INLET_NORMAL); KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, Y_PLUS) + + // Aerodynamic coefficients + KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, LIFT_COEFFICIENT); + KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, DRAG_COEFFICIENT); + KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, LATERAL_FORCE_COEFFICIENT); + KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, ROLLING_MOMENT_COEFFICIENT); + KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, PITCHING_MOMENT_COEFFICIENT); + KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, YAWING_MOMENT_COEFFICIENT); } diff --git a/applications/FluidDynamicsApplication/fluid_dynamics_application.cpp b/applications/FluidDynamicsApplication/fluid_dynamics_application.cpp index b44ac8900fcb..4a6ed7000340 100644 --- a/applications/FluidDynamicsApplication/fluid_dynamics_application.cpp +++ b/applications/FluidDynamicsApplication/fluid_dynamics_application.cpp @@ -311,6 +311,15 @@ void KratosFluidDynamicsApplication::Register() { KRATOS_REGISTER_VARIABLE( Y_PLUS ) + // Aerodynamic coefficients + KRATOS_REGISTER_VARIABLE( LIFT_COEFFICIENT ) + KRATOS_REGISTER_VARIABLE( DRAG_COEFFICIENT ) + KRATOS_REGISTER_VARIABLE( LATERAL_FORCE_COEFFICIENT ) + KRATOS_REGISTER_VARIABLE( ROLLING_MOMENT_COEFFICIENT ) + KRATOS_REGISTER_VARIABLE( PITCHING_MOMENT_COEFFICIENT ) + KRATOS_REGISTER_VARIABLE( YAWING_MOMENT_COEFFICIENT ) + + // Register Elements KRATOS_REGISTER_ELEMENT("VMS2D3N",mVMS2D); //this is the name the element should have according to the naming convention KRATOS_REGISTER_ELEMENT("VMS3D4N",mVMS3D); //this is the name the element should have according to the naming convention diff --git a/applications/FluidDynamicsApplication/fluid_dynamics_application_variables.cpp b/applications/FluidDynamicsApplication/fluid_dynamics_application_variables.cpp index d91f92fa0625..e9c2d62a6273 100644 --- a/applications/FluidDynamicsApplication/fluid_dynamics_application_variables.cpp +++ b/applications/FluidDynamicsApplication/fluid_dynamics_application_variables.cpp @@ -134,4 +134,11 @@ KRATOS_CREATE_VARIABLE( double, DISTANCE_CORRECTION ) KRATOS_CREATE_VARIABLE( double, Y_PLUS ) +// Aerodynamic coefficients +KRATOS_CREATE_VARIABLE( double, LIFT_COEFFICIENT ) +KRATOS_CREATE_VARIABLE( double, DRAG_COEFFICIENT ) +KRATOS_CREATE_VARIABLE( double, LATERAL_FORCE_COEFFICIENT ) +KRATOS_CREATE_VARIABLE( double, ROLLING_MOMENT_COEFFICIENT ) +KRATOS_CREATE_VARIABLE( double, PITCHING_MOMENT_COEFFICIENT ) +KRATOS_CREATE_VARIABLE( double, YAWING_MOMENT_COEFFICIENT ) } diff --git a/applications/FluidDynamicsApplication/fluid_dynamics_application_variables.h b/applications/FluidDynamicsApplication/fluid_dynamics_application_variables.h index 378e41fabf49..0b8043347431 100644 --- a/applications/FluidDynamicsApplication/fluid_dynamics_application_variables.h +++ b/applications/FluidDynamicsApplication/fluid_dynamics_application_variables.h @@ -142,4 +142,12 @@ KRATOS_DEFINE_APPLICATION_VARIABLE( FLUID_DYNAMICS_APPLICATION, bool, MOMENTUM_C KRATOS_DEFINE_APPLICATION_VARIABLE( FLUID_DYNAMICS_APPLICATION, double, DISTANCE_CORRECTION ) KRATOS_DEFINE_APPLICATION_VARIABLE( FLUID_DYNAMICS_APPLICATION, double, Y_PLUS ) + +// Aerodynamic coefficients +KRATOS_DEFINE_APPLICATION_VARIABLE( FLUID_DYNAMICS_APPLICATION, double, LIFT_COEFFICIENT ) +KRATOS_DEFINE_APPLICATION_VARIABLE( FLUID_DYNAMICS_APPLICATION, double, DRAG_COEFFICIENT ) +KRATOS_DEFINE_APPLICATION_VARIABLE( FLUID_DYNAMICS_APPLICATION, double, LATERAL_FORCE_COEFFICIENT ) +KRATOS_DEFINE_APPLICATION_VARIABLE( FLUID_DYNAMICS_APPLICATION, double, ROLLING_MOMENT_COEFFICIENT ) +KRATOS_DEFINE_APPLICATION_VARIABLE( FLUID_DYNAMICS_APPLICATION, double, PITCHING_MOMENT_COEFFICIENT ) +KRATOS_DEFINE_APPLICATION_VARIABLE( FLUID_DYNAMICS_APPLICATION, double, YAWING_MOMENT_COEFFICIENT ) } diff --git a/applications/FluidDynamicsApplication/python_scripts/compute_aerodynamic_coefficients_process.py b/applications/FluidDynamicsApplication/python_scripts/compute_aerodynamic_coefficients_process.py new file mode 100644 index 000000000000..98acc52ecbea --- /dev/null +++ b/applications/FluidDynamicsApplication/python_scripts/compute_aerodynamic_coefficients_process.py @@ -0,0 +1,8 @@ +import KratosMultiphysics +from KratosMultiphysics.FluidDynamicsApplication import ComputeAerodynamicCoefficientsProcess + + +def Factory(settings, Model): + if not isinstance(settings, KratosMultiphysics.Parameters): + raise Exception("expected input shall be a Parameters object, encapsulating a json string") + return ComputeAerodynamicCoefficientsProcess(Model, settings["Parameters"])