Kratos Rotating Frame Process

applications/FluidDynamicsApplication/CMakeLists.txt

@@ -12,6 +12,7 @@ set( KRATOS_FLUID_DYNAMICS_APPLICATION_CORE_SOURCES
   ${CMAKE_CURRENT_SOURCE_DIR}/fluid_dynamics_application_variables.cpp
 
   # utilities (compiled first because they are used within the element cpps)
+  ${CMAKE_CURRENT_SOURCE_DIR}/custom_utilities/rotating_frame_utility.cpp
   ${CMAKE_CURRENT_SOURCE_DIR}/custom_utilities/estimate_dt_utilities.cpp
   ${CMAKE_CURRENT_SOURCE_DIR}/custom_utilities/fluid_auxiliary_utilities.cpp
   ${CMAKE_CURRENT_SOURCE_DIR}/custom_utilities/fluid_calculation_utilities.cpp
@@ -176,4 +177,4 @@ endif((${USE_MPI} MATCHES ON) AND (${TRILINOS_FOUND}))
 
 # Define custom targets
 set(KRATOS_KERNEL "${KRATOS_KERNEL};KratosFluidDynamicsCore" PARENT_SCOPE)
-set(KRATOS_PYTHON_INTERFACE "${KRATOS_PYTHON_INTERFACE};KratosFluidDynamicsApplication" PARENT_SCOPE)
+set(KRATOS_PYTHON_INTERFACE "${KRATOS_PYTHON_INTERFACE};KratosFluidDynamicsApplication" PARENT_SCOPE)

applications/FluidDynamicsApplication/custom_python/add_custom_utilities_to_python.cpp

@@ -36,6 +36,7 @@
 #include "custom_utilities/compressible_element_rotation_utility.h"
 #include "custom_utilities/acceleration_limitation_utilities.h"
 #include "custom_utilities/fluid_test_utilities.h"
+#include "custom_utilities/rotating_frame_utility.h"
 
 #include "utilities/split_tetrahedra.h"
 
@@ -199,6 +200,12 @@ void  AddCustomUtilitiesToPython(pybind11::module& m)
         .def_static("RandomFillNonHistoricalVariable", py::overload_cast<ModelPart::ElementsContainerType&, const Variable<array_1d<double, 3>>&, const IndexType, const double, const double>(&FluidTestUtilities::RandomFillNonHistoricalVariable<ModelPart::ElementsContainerType, array_1d<double, 3>>))
         .def_static("RandomFillNonHistoricalVariable", py::overload_cast<ModelPart::ElementsContainerType&, const Variable<array_1d<double, 3>>&, const std::string&, const IndexType, const double, const double>(&FluidTestUtilities::RandomFillNonHistoricalVariable<ModelPart::ElementsContainerType, array_1d<double, 3>>))
         ;
+
+    // Rotating Frame Utility
+    py::class_<RotatingFrameUtility>(m, "RotatingFrameUtility")
+        .def_static("ApplyVelocityToRotatingObject", &RotatingFrameUtility::ApplyVelocityToRotatingObject)
+        .def_static("ApplyRotationAndMeshDisplacement", &RotatingFrameUtility::ApplyRotationAndMeshDisplacement)
+        ;
 }
 
 }  // namespace Python.

applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.cpp

@@ -0,0 +1,120 @@
+//    |  /           |
+//    ' /   __| _` | __|  _ \   __|
+//    . \  |   (   | |   (   |\__ `
+//   _|\_\_|  \__,_|\__|\___/ ____/
+//                   Multi-Physics
+//
+//  License:         BSD License
+//                   Kratos default license: kratos/license.txt
+//
+//  Main author:     Sebastian Ares de Parga Regalado
+//
+
+// System includes
+#include <string>
+#include <iostream>
+
+// External includes
+
+// Application includes
+#include "rotating_frame_utility.h"
+
+namespace Kratos 
+{
+
+void RotatingFrameUtility::ApplyVelocityToRotatingObject(
+    ModelPart& rModelPart,
+    const array_1d<double, 3>& rAxisOfRotation, 
+    const double& rOmega, 
+    const array_1d<double, 3>& rCenterOfRotation)
+{
+    KRATOS_TRY;
+
+    // Creating the angular velocity vector
+    array_1d<double, 3> angular_velocity_vector = rOmega * rAxisOfRotation;
+
+    thread_local array_1d<double, 3> position_vector;
+    thread_local array_1d<double, 3> velocity_vector;
+
+    // Apply the velocity calculations to each node in parallel
+    block_for_each(rModelPart.Nodes(), [&](Node& rNode) {
+        // Getting the current coordinates of the node
+        auto& r_point = rNode.Coordinates();
+
+        // Calculating the position vector (relative to the rotation center)
+        position_vector = r_point - rCenterOfRotation;
+
+        // Computing the velocity due to rotation (v = omega cross r)
+        velocity_vector = MathUtils<double>::CrossProduct(angular_velocity_vector, position_vector);
+
+        // Setting the node's velocity
+        rNode.FastGetSolutionStepValue(VELOCITY_X) = velocity_vector[0];
+        rNode.FastGetSolutionStepValue(VELOCITY_Y) = velocity_vector[1];
+        rNode.FastGetSolutionStepValue(VELOCITY_Z) = velocity_vector[2];
+
+        // Fix the velocity components
+        rNode.Fix(VELOCITY_X);
+        rNode.Fix(VELOCITY_Y);
+        rNode.Fix(VELOCITY_Z);
+    });
+
+    KRATOS_CATCH("");
+}
+
+void RotatingFrameUtility::ApplyRotationAndMeshDisplacement(
+    ModelPart& rRotatingFrameModelPart,
+    const array_1d<double, 3>& rAxisOfRotation, 
+    const double& rTheta, 
+    const array_1d<double, 3>& rCenterOfRotation)
+{
+    KRATOS_TRY;
+
+    // Quaternion constants
+    const double a = std::cos(rTheta / 2);
+    const double b = -rAxisOfRotation[0] * std::sin(rTheta / 2);
+    const double c = -rAxisOfRotation[1] * std::sin(rTheta / 2);
+    const double d = -rAxisOfRotation[2] * std::sin(rTheta / 2);
+
+    // Creating a quaternion rotation matrix
+    BoundedMatrix<double, 3, 3> rot_matrix;
+    rot_matrix(0, 0) = a*a+b*b-c*c-d*d;
+    rot_matrix(0, 1) = 2*(b*c-a*d);
+    rot_matrix(0, 2) = 2*(b*d+a*c);
+    rot_matrix(1, 0) = 2*(b*c+a*d);
+    rot_matrix(1, 1) = a*a+c*c-b*b-d*d;
+    rot_matrix(1, 2) = 2*(c*d-a*b);
+    rot_matrix(2, 0) = 2*(b*d-a*c);
+    rot_matrix(2, 1) = 2*(c*d+a*b);
+    rot_matrix(2, 2) = a*a+d*d-b*b-c*c;
+
+    thread_local array_1d<double, 3> rotated_point;
+
+    // Apply the rotation and mesh displacement calculations to each node in parallel
+    block_for_each(rRotatingFrameModelPart.Nodes(), [&](Node& rNode) {
+        // Getting the initial coordinates of the node
+        auto& r_point = rNode.GetInitialPosition().Coordinates();
+
+        // Shifting the rotation center to the origin
+        auto centered_point = r_point - rCenterOfRotation;
+
+        // Applying the rotation
+        rotated_point = prod(centered_point, rot_matrix);
+
+        // Shifting the point back and updating the node coordinates
+        rotated_point += rCenterOfRotation;
+        rNode.FastGetSolutionStepValue(MESH_DISPLACEMENT_X) = rotated_point[0] - rNode.X0();
+        rNode.FastGetSolutionStepValue(MESH_DISPLACEMENT_Y) = rotated_point[1] - rNode.Y0();
+        rNode.FastGetSolutionStepValue(MESH_DISPLACEMENT_Z) = rotated_point[2] - rNode.Z0();
+        rNode.Fix(MESH_DISPLACEMENT_X);
+        rNode.Fix(MESH_DISPLACEMENT_Y);
+        rNode.Fix(MESH_DISPLACEMENT_Z);
+        rNode.X() = rotated_point[0];
+        rNode.Y() = rotated_point[1];
+        rNode.Z() = rotated_point[2];
+    });
+
+    KRATOS_CATCH("");
+}
+
+}  // namespace Kratos.
+

applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.h

@@ -0,0 +1,113 @@
+//    |  /           |
+//    ' /   __| _` | __|  _ \   __|
+//    . \  |   (   | |   (   |\__ `
+//   _|\_\_|  \__,_|\__|\___/ ____/
+//                   Multi-Physics
+//
+//  License:         BSD License
+//                   Kratos default license: kratos/license.txt
+//
+//  Main author:     Sebastian Ares de Parga Regalado
+//
+
+#if !defined(KRATOS_ROTATING_FRAME_UTILITY_H)
+#define KRATOS_ROTATING_FRAME_UTILITY_H
+
+// System includes
+#include <string>
+#include <iostream>
+
+// External includes
+#include "includes/mesh_moving_variables.h"
+
+// Include kratos definitions
+#include "includes/define.h"
+
+// Application includes
+#include "utilities/variable_utils.h"
+
+// Project includes
+#include "utilities/builtin_timer.h"
+#include "utilities/math_utils.h"
+
+namespace Kratos {
+///@addtogroup FluidDynamicsApplication
+///@{
+
+///@name Kratos classes
+///@{
+
+/**
+ * @class RotateFrameUtility
+ * @ingroup Kratos Core
+ * @brief Utility for rotating a frame with a rotating object.
+ * @details A utility that rotates a frame (model part) with a rotating object (solid)
+ * based on the specified angular velocity and axis of rotation. It enables the
+ * transformation of coordinates and displacements between the rotating and
+ * non-rotating frames, ensuring consistent updates of the frame's entities
+ * during simulation. Designed for use within the Kratos framework for various
+ * applications such as sliding mesh problems and rotating machinery simulations.
+ * @note Thread-safe and supports parallel execution using Kratos block for each.
+ * @author Sebastian Ares de Parga Regalado
+*/
+
+class KRATOS_API(FLUID_DYNAMICS_APPLICATION) RotatingFrameUtility
+{
+    public:
+      ///@name Type Definitions
+      ///@{
+
+      ///@}
+      ///@name Static Operations
+      ///@{
+
+      /**
+       * @brief Applies an angular velocity to nodes in a rotating frame.
+       *
+       * This function applies an angular velocity to each node in a provided model part. 
+       * The angular velocity vector is calculated from a provided rotation axis and rotation speed. 
+       * The position vector for each node is determined relative to a provided rotation center. 
+       * The velocity due to rotation is then computed as the cross product of the angular velocity vector and the position vector.
+       * The calculated velocity is set as the solution step value for each node. 
+       * Finally, the function fixes the velocity components to ensure they remain constant during the simulation.
+       *
+       * @param rModelPart The model part in which to apply the velocity.
+       * @param rAxisOfRotation The axis of rotation.
+       * @param rOmega The rotation speed.
+       * @param rCenterOfRotation The center of rotation.
+       */
+      static void ApplyVelocityToRotatingObject(
+          ModelPart& rModelPart,
+          const array_1d<double, 3>& rAxisOfRotation, 
+          const double& rOmega, 
+          const array_1d<double, 3>& rCenterOfRotation);
+
+
+      /**
+       * @brief Apply rotation and mesh displacement to a rotating frame.
+       * This function applies rotation and mesh displacement to a rotating frame model 
+       * part based on the provided rotation axis, angle and rotation center. 
+       * It constructs a rotation matrix based on the rotation axis and angle. 
+       * The rotation is applied to each node, which is then displaced 
+       * by the difference between the rotated coordinates and the initial coordinates. 
+       * The solution step values for the mesh displacement components are set accordingly, 
+       * and the mesh displacement components are fixed to remain constant during the simulation. 
+       * Finally, the node coordinates are updated with the rotated coordinates.
+       * @param rRotatingFrameModelPart The rotating frame model part to which the rotation and mesh displacement will be applied.
+       * @param rAxisOfRotation The rotation axis vector.
+       * @param rTheta The rotation angle.
+       * @param rCenterOfRotation The center of rotation.
+       */
+      static void ApplyRotationAndMeshDisplacement(
+          ModelPart& rRotatingFrameModelPart,
+          const array_1d<double, 3>& rAxisOfRotation, 
+          const double& rTheta, 
+          const array_1d<double, 3>& rCenterOfRotation);
+
+}; // Class RotatingFrameUtility
+
+///@}
+
+}  // namespace Kratos.
+
+#endif // KRATOS_ROTATING_FRAME_UTILITY_H