Merge pull request #12889 from KratosMultiphysics/fluid/weakly-compre…

…ssible-darcy-relative-velocity

[Fluid] Weakly compressible darcy relative velocity

applications/FluidDynamicsApplication/custom_python/fluid_dynamics_python_application.cpp

@@ -64,6 +64,7 @@ PYBIND11_MODULE(KratosFluidDynamicsApplication,m)
 
     // Darcy's flow variables
     KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, RESISTANCE);
+    KRATOS_REGISTER_IN_PYTHON_3D_VARIABLE_WITH_COMPONENTS(m,SOLID_FRACTION_VELOCITY);
 
     // Wall modelling
     KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, SLIP_TANGENTIAL_CORRECTION_SWITCH)

applications/FluidDynamicsApplication/custom_utilities/weakly_compressible_navier_stokes_data.h

@@ -58,6 +58,7 @@ NodalVectorData Velocity_OldStep1;
 NodalVectorData Velocity_OldStep2;
 NodalVectorData MeshVelocity;
 NodalVectorData BodyForce;
+NodalVectorData SolidFractionVelocity;
 
 NodalScalarData Pressure;
 NodalScalarData Pressure_OldStep1;
@@ -100,6 +101,7 @@ void Initialize(const Element& rElement, const ProcessInfo& rProcessInfo) overri
     this->FillFromHistoricalNodalData(Density,DENSITY,r_geometry);
     this->FillFromHistoricalNodalData(Pressure_OldStep1,PRESSURE,r_geometry,1);
     this->FillFromHistoricalNodalData(Pressure_OldStep2,PRESSURE,r_geometry,2);
+    this->FillFromHistoricalNodalData(SolidFractionVelocity, SOLID_FRACTION_VELOCITY, r_geometry);
     this->FillFromNonHistoricalNodalData(SoundVelocity, SOUND_VELOCITY, r_geometry);
     this->FillFromProperties(DynamicViscosity,DYNAMIC_VISCOSITY,r_properties); //TODO: This needs to be retrieved from the EffectiveViscosity of the constitutive law
     this->FillFromProcessInfo(DeltaTime,DELTA_TIME,rProcessInfo);
@@ -136,6 +138,7 @@ static int Check(const Element& rElement, const ProcessInfo& rProcessInfo)
         KRATOS_CHECK_VARIABLE_IN_NODAL_DATA(MESH_VELOCITY,r_geometry[i]);
         KRATOS_CHECK_VARIABLE_IN_NODAL_DATA(BODY_FORCE,r_geometry[i]);
         KRATOS_CHECK_VARIABLE_IN_NODAL_DATA(PRESSURE,r_geometry[i]);
+        KRATOS_CHECK_VARIABLE_IN_NODAL_DATA(DENSITY,r_geometry[i]);
     }
 
     return 0;

applications/FluidDynamicsApplication/fluid_dynamics_application_variables.cpp

@@ -40,6 +40,7 @@ KRATOS_CREATE_3D_VARIABLE_WITH_COMPONENTS(AUXILIAR_VECTOR_VELOCITY)
 
 // Darcy's flow variables
 KRATOS_CREATE_VARIABLE(double, RESISTANCE)
+KRATOS_CREATE_3D_VARIABLE_WITH_COMPONENTS(SOLID_FRACTION_VELOCITY)
 
 // Wall modelling
 KRATOS_CREATE_VARIABLE(bool, SLIP_TANGENTIAL_CORRECTION_SWITCH)

applications/FluidDynamicsApplication/fluid_dynamics_application_variables.h

@@ -50,6 +50,7 @@ KRATOS_DEFINE_APPLICATION_VARIABLE( FLUID_DYNAMICS_APPLICATION, double, FS_PRESS
 
 // Darcy's flow variables
 KRATOS_DEFINE_APPLICATION_VARIABLE(FLUID_DYNAMICS_APPLICATION, double, RESISTANCE)
+KRATOS_DEFINE_3D_APPLICATION_VARIABLE_WITH_COMPONENTS(FLUID_DYNAMICS_APPLICATION, SOLID_FRACTION_VELOCITY)
 
 // Wall modelling
 KRATOS_DEFINE_APPLICATION_VARIABLE(FLUID_DYNAMICS_APPLICATION, bool, SLIP_TANGENTIAL_CORRECTION_SWITCH)

applications/FluidDynamicsApplication/python_scripts/navier_stokes_embedded_solver.py

@@ -26,8 +26,8 @@ def __init__(self, formulation_settings):
         self.condition_name = None
         self.process_info_data = {}
         self.element_has_nodal_properties = False
-        self.historical_nodal_properties_variables_list = []
-        self.non_historical_nodal_properties_variables_list = []
+        self.historical_nodal_variables_list = []
+        self.non_historical_nodal_variables_list = []
 
         if formulation_settings.Has("element_type"):
             element_type = formulation_settings["element_type"].GetString()
@@ -96,8 +96,8 @@ def _SetUpEmbeddedWeaklyCompressibleNavierStokes(self, formulation_settings):
         self.level_set_type = formulation_settings["level_set_type"].GetString()
         self.element_integrates_in_time = True
         self.element_has_nodal_properties = True
-        self.historical_nodal_properties_variables_list = [KratosMultiphysics.DENSITY]
-        self.non_historical_nodal_properties_variables_list = [KratosMultiphysics.SOUND_VELOCITY]
+        self.historical_nodal_variables_list = [KratosMultiphysics.DENSITY, KratosCFD.SOLID_FRACTION_VELOCITY]
+        self.non_historical_nodal_variables_list = [KratosMultiphysics.SOUND_VELOCITY]
 
         self.process_info_data[KratosMultiphysics.DYNAMIC_TAU] = formulation_settings["dynamic_tau"].GetDouble()
         self.process_info_data[KratosMultiphysics.PENALTY_COEFFICIENT] = formulation_settings["penalty_coefficient"].GetDouble()
@@ -143,8 +143,8 @@ def _SetUpEmbeddedWeaklyCompressibleNavierStokesDiscontinuous(self, formulation_
         self.level_set_type = formulation_settings["level_set_type"].GetString()
         self.element_integrates_in_time = True
         self.element_has_nodal_properties = True
-        self.historical_nodal_properties_variables_list = [KratosMultiphysics.DENSITY]
-        self.non_historical_nodal_properties_variables_list = [KratosMultiphysics.SOUND_VELOCITY]
+        self.historical_nodal_variables_list = [KratosMultiphysics.DENSITY, KratosCFD.SOLID_FRACTION_VELOCITY]
+        self.non_historical_nodal_variables_list = [KratosMultiphysics.SOUND_VELOCITY]
 
         self.process_info_data[KratosMultiphysics.DYNAMIC_TAU] = formulation_settings["dynamic_tau"].GetDouble()
         self.process_info_data[KratosMultiphysics.PENALTY_COEFFICIENT] = formulation_settings["penalty_coefficient"].GetDouble()
@@ -323,8 +323,8 @@ def __init__(self, model, custom_settings):
         self.level_set_type = self.embedded_formulation.level_set_type
         self.element_integrates_in_time = self.embedded_formulation.element_integrates_in_time
         self.element_has_nodal_properties = self.embedded_formulation.element_has_nodal_properties
-        self.historical_nodal_properties_variables_list = self.embedded_formulation.historical_nodal_properties_variables_list
-        self.non_historical_nodal_properties_variables_list = self.embedded_formulation.non_historical_nodal_properties_variables_list
+        self.historical_nodal_variables_list = self.embedded_formulation.historical_nodal_variables_list
+        self.non_historical_nodal_variables_list = self.embedded_formulation.non_historical_nodal_variables_list
 
         ## Set the distance reading filename
         # TODO: remove the manual "distance_file_name" set as soon as the problem type one has been tested.
@@ -359,10 +359,9 @@ def AddVariables(self):
             self.main_model_part.AddNodalSolutionStepVariable(KratosMultiphysics.MESH_DISPLACEMENT)
             self.main_model_part.AddNodalSolutionStepVariable(KratosMultiphysics.MESH_REACTION)
 
-        # Adding variables required for the nodal material properties
-        if self.element_has_nodal_properties:
-            for variable in self.historical_nodal_properties_variables_list:
-                self.main_model_part.AddNodalSolutionStepVariable(variable)
+        # Adding variables required by the formulation (this includes the nodal material properties)
+        for variable in self.historical_nodal_variables_list:
+            self.main_model_part.AddNodalSolutionStepVariable(variable)
 
         KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__, "Fluid solver variables added correctly.")
 
@@ -498,8 +497,8 @@ def _SetPhysicalProperties(self):
         return materials_imported
 
     def _SetNodalProperties(self):
-        set_density = KratosMultiphysics.DENSITY in self.historical_nodal_properties_variables_list
-        set_sound_velocity = KratosMultiphysics.SOUND_VELOCITY in self.non_historical_nodal_properties_variables_list
+        set_density = KratosMultiphysics.DENSITY in self.historical_nodal_variables_list
+        set_sound_velocity = KratosMultiphysics.SOUND_VELOCITY in self.non_historical_nodal_variables_list
 
         # Get density and dynamic viscostity from the properties of the first element
         for el in self.main_model_part.Elements:

applications/FluidDynamicsApplication/python_scripts/navier_stokes_monolithic_solver.py

@@ -13,8 +13,8 @@ def __init__(self,settings):
         self.element_name = None
         self.element_integrates_in_time = False
         self.element_has_nodal_properties = False
-        self.historical_nodal_properties_variables_list = []
-        self.non_historical_nodal_properties_variables_list = []
+        self.historical_nodal_variables_list = []
+        self.non_historical_nodal_variables_list = []
         self.process_data = {}
 
         #TODO: Keep this until the MonolithicWallCondition is removed to ensure backwards compatibility in solvers with no defined condition_name
@@ -83,7 +83,7 @@ def _SetUpClassicVMS(self,settings):
 
         # set the nodal material properties flag
         self.element_has_nodal_properties = True
-        self.historical_nodal_properties_variables_list = [KratosMultiphysics.DENSITY, KratosMultiphysics.VISCOSITY]
+        self.historical_nodal_variables_list = [KratosMultiphysics.DENSITY, KratosMultiphysics.VISCOSITY]
 
         # validate the non-newtonian parameters if necessary
         if self.non_newtonian_option:
@@ -171,8 +171,8 @@ def _SetUpWeaklyCompressible(self,settings):
 
         # set the nodal material properties flag
         self.element_has_nodal_properties = True
-        self.historical_nodal_properties_variables_list = [KratosMultiphysics.DENSITY]
-        self.non_historical_nodal_properties_variables_list = [KratosMultiphysics.SOUND_VELOCITY]
+        self.historical_nodal_variables_list = [KratosMultiphysics.DENSITY, KratosCFD.SOLID_FRACTION_VELOCITY]
+        self.non_historical_nodal_variables_list = [KratosMultiphysics.SOUND_VELOCITY]
 
         self.process_data[KratosMultiphysics.DYNAMIC_TAU] = settings["dynamic_tau"].GetDouble()
         #TODO: Remove SOUND_VELOCITY from ProcessInfo. Should be obtained from the properties.
@@ -300,10 +300,9 @@ def AddVariables(self):
         self.main_model_part.AddNodalSolutionStepVariable(KratosMultiphysics.Y_WALL)
         self.main_model_part.AddNodalSolutionStepVariable(KratosCFD.Q_VALUE)
 
-        # Adding variables required for the nodal material properties
-        if self.element_has_nodal_properties:
-            for variable in self.historical_nodal_properties_variables_list:
-                self.main_model_part.AddNodalSolutionStepVariable(variable)
+        # Adding variables required by the formulation (this includes the nodal material properties)
+        for variable in self.historical_nodal_variables_list:
+            self.main_model_part.AddNodalSolutionStepVariable(variable)
 
         # Adding variables required for the periodic conditions
         if self.settings["consider_periodic_conditions"].GetBool() == True:
@@ -348,8 +347,8 @@ def _SetFormulation(self):
         self.condition_name = self.formulation.condition_name
         self.element_integrates_in_time = self.formulation.element_integrates_in_time
         self.element_has_nodal_properties = self.formulation.element_has_nodal_properties
-        self.historical_nodal_properties_variables_list = self.formulation.historical_nodal_properties_variables_list
-        self.non_historical_nodal_properties_variables_list = self.formulation.non_historical_nodal_properties_variables_list
+        self.historical_nodal_variables_list = self.formulation.historical_nodal_variables_list
+        self.non_historical_nodal_variables_list = self.formulation.non_historical_nodal_variables_list
 
     def _SetTimeSchemeBufferSize(self):
         scheme_type = self.settings["time_scheme"].GetString()
@@ -373,9 +372,9 @@ def _SetUpSteadySimulation(self):
             self.settings["formulation"]["dynamic_tau"].SetDouble(0.0)
 
     def _SetNodalProperties(self):
-        set_density = KratosMultiphysics.DENSITY in self.historical_nodal_properties_variables_list
-        set_viscosity = KratosMultiphysics.VISCOSITY in self.historical_nodal_properties_variables_list
-        set_sound_velocity = KratosMultiphysics.SOUND_VELOCITY in self.non_historical_nodal_properties_variables_list
+        set_density = KratosMultiphysics.DENSITY in self.historical_nodal_variables_list
+        set_viscosity = KratosMultiphysics.VISCOSITY in self.historical_nodal_variables_list
+        set_sound_velocity = KratosMultiphysics.SOUND_VELOCITY in self.non_historical_nodal_variables_list
 
         # Get density and dynamic viscostity from the properties of the first element
         for el in self.main_model_part.Elements:

applications/FluidDynamicsApplication/python_scripts/symbolic_generation/weakly_compressible_navier_stokes/generate_weakly_compressible_navier_stokes_element.py

@@ -118,6 +118,7 @@
 
     ## Other data definitions
     f = DefineMatrix('f',nnodes,dim)            # Forcing term
+    v_sol_frac = DefineMatrix('r_v_sol_frac',nnodes,dim) # Solid fraction velocity
 
     ## Constitutive matrix definition
     C = DefineSymmetricMatrix('C',strain_size,strain_size)
@@ -146,6 +147,7 @@
     ## Data interpolation to the Gauss points
     f_gauss = f.transpose()*N
     v_gauss = v.transpose()*N
+    v_sol_frac_gauss = v_sol_frac.transpose()*N
 
     ## Convective velocity definition
     if convective_term:
@@ -164,10 +166,10 @@
         for i in range(0, dim):
             stab_norm_a += vconv_gauss[i]**2
         stab_norm_a = sympy.sqrt(stab_norm_a)
-        tau1 = 1.0/(rho*dyn_tau/dt + stab_c2*rho*stab_norm_a/h + stab_c1*mu/h**2 + stab_c3*sigma/h**2) # Stabilization parameter 1
+        tau1 = 1.0/(rho*dyn_tau/dt + stab_c2*rho*stab_norm_a/h + stab_c1*mu/h**2 + stab_c3*sigma/h)    # Stabilization parameter 1
         tau2 = mu + (stab_c2*rho*stab_norm_a*h + stab_c3*sigma)/stab_c1                                # Stabilization parameter 2
     else:
-        tau1 = 1.0/(rho*dyn_tau/dt + stab_c1*mu/h**2 + stab_c3*sigma/h**2) # Stabilization parameter 1
+        tau1 = 1.0/(rho*dyn_tau/dt + stab_c1*mu/h**2 + stab_c3*sigma/h)    # Stabilization parameter 1
         tau2 = mu + stab_c3*sigma/stab_c1                                  # Stabilization parameter 2
 
     ## Compute the rest of magnitudes at the Gauss points
@@ -205,7 +207,7 @@
     ## Compute galerkin functional
     # Navier-Stokes functional
     if divide_by_rho:
-        rv_galerkin = rho*w_gauss.transpose()*f_gauss - rho*w_gauss.transpose()*accel_gauss - grad_w_voigt.transpose()*stress + div_w*p_gauss - sigma*w_gauss.transpose()*v_gauss - q_gauss*div_v
+        rv_galerkin = rho*w_gauss.transpose()*f_gauss - rho*w_gauss.transpose()*accel_gauss - grad_w_voigt.transpose()*stress + div_w*p_gauss - sigma*w_gauss.transpose()*(v_gauss-v_sol_frac_gauss) - q_gauss*div_v
         if artificial_compressibility:
             rv_galerkin -= (1/(rho*c*c))*q_gauss*pder_gauss
             if convective_term:
@@ -224,7 +226,7 @@
     ##  Stabilization functional terms
     # Momentum conservation residual
     # Note that the viscous stress term is dropped since linear elements are used
-    vel_residual = rho*f_gauss - rho*accel_gauss - grad_p - sigma*v_gauss
+    vel_residual = rho*f_gauss - rho*accel_gauss - grad_p - sigma*(v_gauss-v_sol_frac_gauss)
     if convective_term:
         vel_residual -= rho*convective_term_gauss.transpose()
 

applications/FluidDynamicsApplication/python_scripts/symbolic_generation/weakly_compressible_navier_stokes/weakly_compressible_navier_stokes_cpp_template.cpp

@@ -332,6 +332,7 @@ void WeaklyCompressibleNavierStokes<WeaklyCompressibleNavierStokesData<2,3>>::Co
     const BoundedMatrix<double,2,3>& vmesh = rData.MeshVelocity;
     const BoundedMatrix<double,2,3> vconv = v - vmesh;
     const BoundedMatrix<double,2,3>& f = rData.BodyForce;
+    const BoundedMatrix<double,2,3>& r_v_sol_frac = rData.SolidFractionVelocity;
     const array_1d<double,3>& p = rData.Pressure;
     const array_1d<double,3>& pn = rData.Pressure_OldStep1;
     const array_1d<double,3>& pnn = rData.Pressure_OldStep2;
@@ -376,6 +377,7 @@ void WeaklyCompressibleNavierStokes<WeaklyCompressibleNavierStokesData<3,4>>::Co
     const BoundedMatrix<double,3,4>& vmesh = rData.MeshVelocity;
     const BoundedMatrix<double,3,4> vconv = v - vmesh;
     const BoundedMatrix<double,3,4>& f = rData.BodyForce;
+    const BoundedMatrix<double,3,4>& r_v_sol_frac = rData.SolidFractionVelocity;
     const array_1d<double,4>& p = rData.Pressure;
     const array_1d<double,4>& pn = rData.Pressure_OldStep1;
     const array_1d<double,4>& pnn = rData.Pressure_OldStep2;

applications/FluidDynamicsApplication/tests/cpp_tests/test_drag_utils.cpp

@@ -53,6 +53,7 @@ namespace Kratos {
             rModelPart.AddNodalSolutionStepVariable(SOUND_VELOCITY);
             rModelPart.AddNodalSolutionStepVariable(DYNAMIC_VISCOSITY);
             rModelPart.AddNodalSolutionStepVariable(REACTION_WATER_PRESSURE);
+            rModelPart.AddNodalSolutionStepVariable(SOLID_FRACTION_VELOCITY);
 
             // Process info creation
             double delta_time = 0.1;

applications/FluidDynamicsApplication/tests/cpp_tests/test_embedded_fluid_element.cpp

@@ -42,6 +42,7 @@ KRATOS_TEST_CASE_IN_SUITE(EmbeddedElement2D3N, FluidDynamicsApplicationFastSuite
     model_part.AddNodalSolutionStepVariable(ACCELERATION);
     model_part.AddNodalSolutionStepVariable(EXTERNAL_PRESSURE);
     model_part.AddNodalSolutionStepVariable(REACTION_WATER_PRESSURE);
+    model_part.AddNodalSolutionStepVariable(SOLID_FRACTION_VELOCITY);
 
     // For VMS comparison
     model_part.AddNodalSolutionStepVariable(NODAL_AREA);

applications/FluidDynamicsApplication/tests/cpp_tests/test_embedded_fluid_element_discontinuous.cpp

@@ -42,6 +42,7 @@ KRATOS_TEST_CASE_IN_SUITE(EmbeddedElementDiscontinuous2D3N, FluidDynamicsApplica
     model_part.AddNodalSolutionStepVariable(ACCELERATION);
     model_part.AddNodalSolutionStepVariable(EXTERNAL_PRESSURE);
     model_part.AddNodalSolutionStepVariable(REACTION_WATER_PRESSURE);
+    model_part.AddNodalSolutionStepVariable(SOLID_FRACTION_VELOCITY);
 
     // For VMS comparison
     model_part.AddNodalSolutionStepVariable(NODAL_AREA);

applications/FluidDynamicsApplication/tests/cpp_tests/test_weakly_compressible_navier_stokes_element.cpp

@@ -43,6 +43,7 @@ KRATOS_TEST_CASE_IN_SUITE(WeaklyCompressibleNavierStokes2D3N, FluidDynamicsAppli
     r_model_part.AddNodalSolutionStepVariable(ACCELERATION);
     r_model_part.AddNodalSolutionStepVariable(REACTION);
     r_model_part.AddNodalSolutionStepVariable(REACTION_WATER_PRESSURE);
+    r_model_part.AddNodalSolutionStepVariable(SOLID_FRACTION_VELOCITY);
 
     // ProcessInfo container fill
     double delta_time = 0.1;
@@ -111,8 +112,8 @@ KRATOS_TEST_CASE_IN_SUITE(WeaklyCompressibleNavierStokes2D3N, FluidDynamicsAppli
     // KRATOS_WATCH(row(LHS,0))
 
     // Check values
-    const std::vector<double> rhs_ref = {34.2921172632, -19.1612962435, -0.0212935839184, -58.7701017724, 31.2334272413, -0.0186565395118, -8.37866202262, -77.7742128742, -0.0450498765698}; // WeaklyCompressibleNavierStokes2D3N
-    const std::vector<double> lhs_0_ref = {255.745969647, 0, 0.155863699867, -197.183916289, 250.033529091, 0.105861405575, 52.887840474, -250.033529091, 0.13527448821};                     // WeaklyCompressibleNavierStokes2D3N
+    const std::vector<double> rhs_ref = {35.4884272438, -16.3707126841, -0.0165099840565, -57.1752103975, 34.4226885377, -0.0120809017136, -6.38728003117, -74.1903403375, -0.0564091142299}; // WeaklyCompressibleNavierStokes2D3N
+    const std::vector<double> lhs_0_ref = {247.633101497, 0, 0.17185384644, -201.239603454, 250.033529091, 0.08787249068, 48.8263970168, -250.033529091, 0.137273256531};                     // WeaklyCompressibleNavierStokes2D3N
     KRATOS_EXPECT_VECTOR_NEAR(RHS, rhs_ref, 1.0e-10)
     KRATOS_EXPECT_VECTOR_NEAR(row(LHS,0), lhs_0_ref, 1.0e-8)
 }