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

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,110 @@
+//    |  /           |
+//    ' /   __| _` | __|  _ \   __|
+//    . \  |   (   | |   (   |\__ `
+//   _|\_\_|  \__,_|\__|\___/ ____/
+//                   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;
+
+    // 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)
+        array_1d<double, 3> position_vector = r_point - rCenterOfRotation;
+
+        // Computing the velocity due to rotation (v = omega cross r)
+        array_1d<double, 3> velocity_vector = MathUtils<double>::CrossProduct(angular_velocity_vector, position_vector);
+
+        // Setting the node's velocity
+        rNode.FastGetSolutionStepValue(VELOCITY) = velocity_vector;
+
+        // 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;
+
+    // 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
+        array_1d<double, 3> rotated_point = prod(centered_point, rot_matrix);
+
+        // Shifting the point back and updating the node coordinates
+        rotated_point += rCenterOfRotation;
+
+        rNode.FastGetSolutionStepValue(MESH_DISPLACEMENT) = rotated_point - rNode.GetInitialPosition();
+        rNode.Fix(MESH_DISPLACEMENT_X);
+        rNode.Fix(MESH_DISPLACEMENT_Y);
+        rNode.Fix(MESH_DISPLACEMENT_Z);
+        rNode.Coordinates() = rotated_point;
+    });
+
+    KRATOS_CATCH("");
+}
+
+}  // namespace Kratos.
+

applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.h

@@ -0,0 +1,111 @@
+//    |  /           |
+//    ' /   __| _` | __|  _ \   __|
+//    . \  |   (   | |   (   |\__ `
+//   _|\_\_|  \__,_|\__|\___/ ____/
+//                   Multi-Physics
+//
+//  License:         BSD License
+//                   Kratos default license: kratos/license.txt
+//
+//  Main author:     Sebastian Ares de Parga Regalado
+//
+
+#pragma once
+
+// 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.

applications/FluidDynamicsApplication/python_scripts/rotating_frame_process.py

@@ -0,0 +1,131 @@
+import KratosMultiphysics as KM
+import KratosMultiphysics.FluidDynamicsApplication as KratosCFD
+import numpy as np
+
+def Factory(settings, model):
+    if not isinstance(settings, KM.Parameters):
+        raise Exception("expected input shall be a Parameters object, encapsulating a json string")
+    return RotatingFrameProcess(model, settings["Parameters"])
+
+## All the processes python should be derived from "Process"
+class RotatingFrameProcess(KM.Process):
+    """A Kratos process that rotates a given RotatingFrameModelPart 
+    around a specified axis of rotation. This process uses the Kratos 
+    RotateNodesUtility to rotate the nodes of the RotatingFrameModelPart 
+    and assigns the corresponding MeshDisplacement to the rotated nodes. 
+    The process also assigns the corresponding rotational mesh velocity to 
+    the nodes in a RotatingObjectModelPart. The process takes as input the 
+    model part to be rotated, the axis of rotation vector, the final angular 
+    velocity in radians per second, and the accelerating time in seconds.
+
+    Public member variables:
+    model -- the container of the different model parts.
+    settings -- Kratos parameters containing process settings.
+    """
+
+    def __init__(self, model, settings):
+        """ The default constructor of the class
+
+        Keyword arguments:
+        self -- It signifies an instance of a class.
+        Model -- the container of the different model parts.
+        settings -- Kratos parameters containing process settings.
+        """
+        KM.Process.__init__(self) # calling the baseclass constructor
+
+        default_settings = KM.Parameters("""{
+            "rotating_frame_model_part_name": "",
+            "rotating_object_model_part_name": "",
+            "center_of_rotation": [0.0,0.0,0.0],
+            "axis_of_rotation": [1.0,0.0,0.0],
+            "target_angular_velocity_radians": 0.0,
+            "acceleration_time": 0.0
+        }""")
+
+        # Add missing settings that the user did not provide but that
+        # are necessary for this process
+        settings.ValidateAndAssignDefaults(default_settings)
+        # Alternative:
+        # settings.RecursivelyValidateAndAssignDefaults(default_settings)
+
+        #Assign settings
+        self.model = model
+
+        # Get the rotating frame model part name
+        if not settings["rotating_frame_model_part_name"].GetString():
+            raise Exception("\'rotating_frame_model_part_name\' not provided. Please specify the model part to rotate the frame.")
+        self.rotating_frame_model_part_name = settings["rotating_frame_model_part_name"].GetString() 
+
+        # Get the rotating object model part name
+        if not settings["rotating_object_model_part_name"].GetString():
+            raise Exception("\'rotating_object_model_part_name\' not provided. Please specify the slave model part.")
+        self.rotating_object_model_part_name = settings["rotating_object_model_part_name"].GetString()
+
+        ## Assign rotating frame and object model parts
+        self.rotating_frame_model_part = self.model.GetModelPart(self.rotating_frame_model_part_name)
+        self.rotating_object_model_part = self.model.GetModelPart(self.rotating_object_model_part_name)
+
+        # Get the center of rotation
+        if settings.Has("center_of_rotation"):
+            self.center_of_rotation = np.array(settings["center_of_rotation"].GetVector())
+        else:
+            raise Exception("The center_of_rotation parameter is missing from the settings.")
+
+        # Get the axis of rotation
+        if settings.Has("axis_of_rotation"):
+            self.axis_of_rotation = np.array(settings["axis_of_rotation"].GetVector())
+            if self.axis_of_rotation.size == 0:
+                raise Exception("The axis_of_rotation vector is empty.")
+            axis_norm = np.linalg.norm(self.axis_of_rotation)
+            if not np.isclose(np.linalg.norm(axis_norm), 1.0, rtol=1e-6):
+                KM.Logger.PrintWarning("RotatingFrameProcess", "The axis_of_rotation vector is not a unit vector... normalizing the vector.")
+                self.axis_of_rotation = self.axis_of_rotation / axis_norm
+        else:
+            raise Exception("The axis_of_rotation parameter is missing from the settings.")
+
+        # Get the angular velocity in radians
+        if settings.Has("target_angular_velocity_radians"):
+            self.target_angular_velocity_radians = settings["target_angular_velocity_radians"].GetDouble()
+        else:
+            raise Exception("The target_angular_velocity_radians parameter is missing from the settings.")
+
+        # Get the acceleration time
+        if settings.Has("acceleration_time"):
+            self.acceleration_time = settings["acceleration_time"].GetDouble()
+            if self.acceleration_time < 0.0:
+                raise Exception("The acceleration_time parameter must be non-negative.")
+        else:
+            raise Exception("The acceleration_time parameter is missing from the settings.")
+
+    def ExecuteInitializeSolutionStep(self):
+        """ This method is executed in order to initialize the current step
+
+        Keyword arguments:
+        self -- It signifies an instance of a class.
+        """
+        # Get Current time
+        self.time = self.rotating_frame_model_part.ProcessInfo[KM.TIME]
+
+        # Calculate the angle (theta) for the specified time
+        self.__GetThetaAndOmega() 
+
+        # Apply Rotation and Mesh Displacement to nodes.
+        KratosCFD.RotatingFrameUtility.ApplyRotationAndMeshDisplacement(self.rotating_frame_model_part, self.axis_of_rotation, self.theta, self.center_of_rotation)
+
+        # Apply Velocity to object.
+        KratosCFD.RotatingFrameUtility.ApplyVelocityToRotatingObject(self.rotating_object_model_part, self.axis_of_rotation, self.omega, self.center_of_rotation)
+
+    def __GetThetaAndOmega(self):
+        # Calculate the angular acceleration (alpha)
+        if np.isclose(self.acceleration_time, 0.0):
+            alpha = 0.0
+        else:
+            alpha = self.target_angular_velocity_radians / self.acceleration_time
+
+        # Calculate the angle (theta) based on current time
+        if self.time <= self.acceleration_time:
+            self.omega = alpha * self.time
+            self.theta = 0.5 * alpha * self.time**2
+        else:
+            self.omega = self.target_angular_velocity_radians
+            self.theta = 0.5 * alpha * self.acceleration_time**2 + self.omega * (self.time - self.acceleration_time)