[FluidDynamicsHydraulicsApplication] Adding new hydraulic utilities; inflow detection and artificial viscosity.

applications/FluidDynamicsHydraulicsApplication/custom_python/add_custom_utilities_to_python.cpp

@@ -31,7 +31,9 @@ void AddCustomUtilitiesToPython(pybind11::module& m)
         .def_static("InitialWaterDepth", [](ModelPart &rModelPart){return HydraulicFluidAuxiliaryUtilities::InitialWaterDepth(rModelPart);})
         .def_static("SetInletVelocity", [](ModelPart &rModelPart, double InletVelocity, const Variable<double> &rDistanceVariable){return HydraulicFluidAuxiliaryUtilities::SetInletVelocity(rModelPart, InletVelocity, rDistanceVariable);})
         .def_static("FreeInlet", [](ModelPart &rModelPart){return HydraulicFluidAuxiliaryUtilities::FreeInlet(rModelPart);})
-        .def_static("SetInletFreeSurface", [](ModelPart &rModelPart, const Flags &rSkinFlag, const Variable<double> &rDistanceVariable){return HydraulicFluidAuxiliaryUtilities::SetInletFreeSurface(rModelPart, rSkinFlag, rDistanceVariable);});
+        .def_static("SetInletFreeSurface", [](ModelPart &rModelPart, const Flags &rSkinFlag, const Variable<double> &rDistanceVariable){return HydraulicFluidAuxiliaryUtilities::SetInletFreeSurface(rModelPart, rSkinFlag, rDistanceVariable);})
+        .def_static("CalculateArtificialViscosity", [](ModelPart &rModelPart, double WaterDynamicViscosityMax){ return HydraulicFluidAuxiliaryUtilities::CalculateArtificialViscosity(rModelPart, WaterDynamicViscosityMax); })
+        .def_static("ApplyOutletInflowLimiter", [](ModelPart &rModelPart,Variable<array_1d<double, 3>>& rVariable,Variable<array_1d<double, 3>>& rVariableNormal){ return HydraulicFluidAuxiliaryUtilities::ApplyOutletInflowLimiter(rModelPart, rVariable, rVariableNormal); });
 }
 
 } // namespace Kratos::Python

applications/FluidDynamicsHydraulicsApplication/custom_utilities/hydraulic_fluid_auxiliary_utilities.cpp

@@ -307,4 +307,60 @@ void HydraulicFluidAuxiliaryUtilities::SetInletFreeSurface(ModelPart &rModelPart
     });
 }
 
+void HydraulicFluidAuxiliaryUtilities::CalculateArtificialViscosity(ModelPart &rModelPart,double WaterDynamicViscosityMax){
+    const auto &r_process_info = rModelPart.GetProcessInfo();
+
+    block_for_each(rModelPart.Elements(), [&](Element &rElement)
+    {
+    double artificial_viscosity;
+
+    rElement.Calculate(ARTIFICIAL_DYNAMIC_VISCOSITY,artificial_viscosity ,r_process_info);
+        if (artificial_viscosity > WaterDynamicViscosityMax)
+        {
+            artificial_viscosity = WaterDynamicViscosityMax;
+        }
+        double neg_nodes = 0.0;
+        double pos_nodes=0.0;
+        for (auto &r_node : rElement.GetGeometry())
+        {
+            double distance = r_node.FastGetSolutionStepValue(DISTANCE);
+
+            if (distance > 0)
+            {
+                pos_nodes += 1;
+            }
+            else
+            {
+                neg_nodes += 1;
+            }
+        }
+        if (neg_nodes > 0 && pos_nodes > 0)
+        {
+            artificial_viscosity = 0.0;
+        }
+        rElement.SetValue(ARTIFICIAL_DYNAMIC_VISCOSITY, artificial_viscosity);
+    });
+}
+
+void HydraulicFluidAuxiliaryUtilities::ApplyOutletInflowLimiter(ModelPart &rModelPart,const Variable<array_1d<double, 3>>& rVariable,const Variable<array_1d<double, 3>>& rVariableNormal)
+{
+
+    block_for_each(rModelPart.Nodes(), [&](NodeType& rNode)
+    {
+        array_1d<double, 3>& r_velocity = rNode.FastGetSolutionStepValue(rVariable);
+        // We use a non-historical variable in case rVariableNormal is not the NORMAL variable and an auxiliary variable is used
+        const array_1d<double, 3>& r_normal = rNode.GetValue(rVariableNormal);
+        const double norm_n = norm_2(r_normal);
+        if (norm_n > 0.0) 
+        {
+            array_1d<double, 3> n_unit = r_normal / norm_n;
+            const double aux = inner_prod(r_velocity, n_unit);
+            if (aux < 0.0)
+            {
+                noalias(r_velocity) -= aux * n_unit;
+            }
+        }
+    });
+}
+
 } // namespace Kratos

applications/FluidDynamicsHydraulicsApplication/custom_utilities/hydraulic_fluid_auxiliary_utilities.h

@@ -116,7 +116,23 @@ class KRATOS_API(FLUID_DYNAMICS_HYDRAULICS_APPLICATION) HydraulicFluidAuxiliaryU
      * @param  rDistancesVariable Variable name of the inlet distance.
      */
     static void SetInletFreeSurface(ModelPart &rModelPart, const Flags &rSkinFlag, const Variable<double> &rDistanceVariable);
-
+
+    /**
+     * @brief  Artificial viscosity is calculated. The purpose of adding this artificial viscosity is to avoid non-physical spikes in velocities.
+     * @param  rModelPart Fluid Model Part
+     * @param  WaterDynamicViscosityMax It is a threshold value to prevent adding excessive artificial numerical viscosity and thereby losing the real physics.
+     */
+    static void CalculateArtificialViscosity(ModelPart &rModelPart, double WaterDynamicViscosityMax);
+
+    /**
+     * @brief  When there is inflow on a boundary considered as an outlet, this function retains only the tangential component, preventing inflows that cause instabilities.
+     * @param  rModelPart Fluid Model Part
+     * @param  rVariable it is possible to use the variable VELOCITY_FRACTIONAL or VELOCITY
+     * @param  rVariableNormal it is possible to use an auxiliar normal such as INLET_NORMAL 
+     */
+    static void ApplyOutletInflowLimiter(ModelPart &rModelPart,
+    const Variable<array_1d<double,3>>& rVariable,
+    const Variable<array_1d<double, 3>>& rVariableNormal);    
 
     ///@}
 

applications/FluidDynamicsHydraulicsApplication/tests/test_hydraulic_fluid_auxiliary_utilities.py

@@ -1,7 +1,6 @@
 import KratosMultiphysics as Kratos
 import KratosMultiphysics.FluidDynamicsHydraulicsApplication as KratosFluidHydraulics
 import KratosMultiphysics.KratosUnittest as UnitTest
-import KratosMultiphysics.kratos_utilities as KratosUtils
 
 class HydraulicFluidAuxiliaryUtilitiesTest(UnitTest.TestCase):
 
@@ -76,5 +75,18 @@ def testCalculateWettedArea(self):
         theoretical_wetted_area= 2.0*(level_set_z)
         self.assertAlmostEqual(WettedArea, theoretical_wetted_area, 12)
 
+    def testAvoidOutletInflow(self):
+        inlet_skin_model_part = self.model.GetModelPart("TestModelPart.Inlet")
+        for node in inlet_skin_model_part.Nodes:
+            node.SetSolutionStepValue(Kratos.VELOCITY, [0.0,0.0,1.0])
+            node.SetValue(Kratos.NORMAL,[0.0,0.0,-1.0])
+
+        KratosFluidHydraulics.HydraulicFluidAuxiliaryUtilities.ApplyOutletInflowLimiter(inlet_skin_model_part,Kratos.VELOCITY,Kratos.NORMAL)
+
+        for node in inlet_skin_model_part.Nodes:
+            if node.Id==5:
+               v=node.GetSolutionStepValue(Kratos.VELOCITY_Z)
+        self.assertAlmostEqual(v, 0.0, 12)
+
 if __name__ == '__main__':
     UnitTest.main()