fluid shell interaction model.

CMakeLists.txt

@@ -54,10 +54,11 @@ target_link_libraries(sphinxsys_core INTERFACE ${Simbody_LIBS})
 target_compile_definitions(sphinxsys_core INTERFACE _SILENCE_CXX17_ITERATOR_BASE_CLASS_DEPRECATION_WARNING) # Because of Simbody headers
 ### Eigen
 find_package(Eigen3 CONFIG REQUIRED)
-target_link_libraries(sphinxsys_core INTERFACE Eigen3::Eigen)
+target_link_libraries (sphinxsys_core INTERFACE Eigen3::Eigen)
 ### TBB
 find_package(TBB CONFIG REQUIRED)
 target_link_libraries(sphinxsys_core INTERFACE TBB::tbb TBB::tbbmalloc $<$<PLATFORM_ID:Windows>:TBB::tbbmalloc_proxy>) 
+add_compile_definitions(TBB_SUPPRESS_DEPRECATED_MESSAGES=1)
 ### Threads
 find_package(Threads REQUIRED)
 target_link_libraries(sphinxsys_core INTERFACE Threads::Threads)

PythonScriptStore/RegressionTest/regression_test_base_tool.py

@@ -16,30 +16,30 @@ def __init__(self, casename, bodyname, parametername):
         self.sphinxsys_case_name = casename
         self.sphinxsys_body_name = bodyname
         self.sphinxsys_parameter_name = parametername
-        self.enter_sphinxsys_exec_folder = f"cd {self.sphinxsys_exec_path};"
-        self.enter_sphinxsys_case_folder = f"cd {self.sphinxsys_case_path};"
+        self.enter_sphinxsys_exec_folder = f"{self.sphinxsys_exec_path}"
+        self.enter_sphinxsys_case_folder = f"{self.sphinxsys_case_path}"
         self.input_file_path = os.path.join(self.sphinxsys_exec_path, "input")
         self.condition_file_path = os.path.join(self.input_file_path, f"{bodyname}_{parametername}_runtimes.dat")
 
     def compile_case(self) -> None:
         print('Start compiling test case....')
         command = "make -j8"
-        os.system(self.enter_sphinxsys_case_folder)
+        os.chdir(self.enter_sphinxsys_case_folder)
         os.system(command)
         print('Compiling test case is finished...')
 
     def run_particle_relaxation(self) -> None:
         print('Start particle relaxation for the simulation...')
         command = f".{os.sep}{self.sphinxsys_case_name} --r=true"
-        os.system(self.enter_sphinxsys_exec_folder)
+        os.chdir(self.enter_sphinxsys_exec_folder)
         os.system(command)
         print('Simulating case is finished...')
 
     def run_case(self) -> None:
         print('Start case simulation...')
         print(self.enter_sphinxsys_exec_folder)
         command = f".{os.sep}{self.sphinxsys_case_name} --r=false --i=true --rt=true"
-        os.system(self.enter_sphinxsys_exec_folder)
+        os.chdir(self.enter_sphinxsys_exec_folder)
         os.system(command)
         print('Simulating case is finished...')
 
@@ -60,44 +60,44 @@ def __init__(self, casename, bodyname, parametername):
         self.sphinxsys_case_name = casename
         self.sphinxsys_body_name = bodyname
         self.sphinxsys_parameter_name = parametername
-        self.enter_sphinxsys_exec_folder = f"cd {self.sphinxsys_exec_path};"
-        self.enter_sphinxsys_case_folder = f"cd {self.sphinxsys_case_path};"
+        self.enter_sphinxsys_exec_folder = f"{self.sphinxsys_exec_path}"
+        self.enter_sphinxsys_case_folder = f"{self.sphinxsys_case_path}"
         self.input_file_path = os.path.join(self.sphinxsys_exec_path, "input")
         self.condition_file_path = os.path.join(self.input_file_path, f"{bodyname}_{parametername}_runtimes.dat")
 
     def compile_case(self) -> None:
         print('Start compiling test case....')
         command = "make -j8"
-        os.system(self.enter_sphinxsys_case_folder)
+        os.chdir(self.enter_sphinxsys_case_folder)
         os.system(command)
         print('Compiling test case is finished...')
 
     def test_case(self) -> None:
         print('Start test case...')
         command = "make test"
-        os.system(self.enter_sphinxsys_case_folder)
+        os.chdir(self.enter_sphinxsys_case_folder)
         os.system(command)
         print('Testing case is finished...')
 
     def copy_reload(self) -> None:
         print('Start copy the reload file...')
         command = "cp -r reload bin"
-        os.system(self.enter_sphinxsys_case_folder)
+        os.chdir(self.enter_sphinxsys_case_folder)
         os.system(command)
         print('Copying the reload file is finished...')
 
     def run_particle_relaxation(self) -> None:
         print('Start particle relaxation for the simulation...')
         command = f".{os.sep}{self.sphinxsys_case_name} --r=true"
-        os.system(self.enter_sphinxsys_exec_folder)
+        os.chdir(self.enter_sphinxsys_exec_folder)
         os.system(command)
         print('Simulating case is finished...')
 
     def run_case(self) -> None:
         print('Start case simulation...')
         print(self.enter_sphinxsys_exec_folder)
         command = f".{os.sep}{self.sphinxsys_case_name} --r=false --i=true --rt=true"
-        os.system(self.enter_sphinxsys_exec_folder)
+        os.chdir(self.enter_sphinxsys_exec_folder)
         os.system(command)
         print('Simulating case is finished...')
 

SPHINXsys/src/shared/particle_dynamics/dissipation_dynamics/particle_dynamics_dissipation.h

@@ -182,11 +182,36 @@ namespace SPH
 		void interaction(size_t index_i, Real dt = 0.0);
 
 	private:
-		Real eta_; /**< damping coefficient */
 		StdLargeVec<Real> &Vol_, &mass_;
 		StdLargeVec<VariableType> &variable_;
 		StdVec<StdLargeVec<Real> *> wall_Vol_;
 		StdVec<StdLargeVec<VariableType> *> wall_variable_;
+
+	protected:
+		Real eta_; /**< damping coefficient */
+	};
+
+	/**
+	 * @class DampingPairwiseFromShell
+	 * @brief Damping to wall by which the wall velocity is not updated
+	 * and the mass of wall particle is not considered.
+	 */
+	template <typename VariableType>
+	class DampingPairwiseFromShell : public LocalDynamics,
+									public DataDelegateContact<BaseParticles, ShellParticles>
+	{
+	public:
+		DampingPairwiseFromShell(BaseContactRelation &contact_relation, const std::string &variable_name, Real eta);
+		virtual ~DampingPairwiseFromShell(){};
+		void interaction(size_t index_i, Real dt = 0.0);
+
+	private:
+		StdLargeVec<Real> &Vol_, &mass_;
+		StdLargeVec<VariableType> &variable_;
+		StdVec<StdLargeVec<VariableType> *> shell_variable_;
+
+	protected:
+		Real eta_; /**< damping coefficient */
 	};
 
 	/**

SPHINXsys/src/shared/particle_dynamics/dissipation_dynamics/particle_dynamics_dissipation.hpp

@@ -427,6 +427,55 @@ namespace SPH
 		}
 	}
 	//=================================================================================================//
+	template <typename VariableType>
+	DampingPairwiseFromShell<VariableType>::
+		DampingPairwiseFromShell(BaseContactRelation &contact_relation, const std::string &variable_name, Real eta)
+		: LocalDynamics(contact_relation.sph_body_),
+		  DataDelegateContact<BaseParticles, ShellParticles>(contact_relation),
+		  eta_(eta), Vol_(particles_->Vol_), mass_(particles_->mass_),
+		  variable_(*particles_->getVariableByName<VariableType>(variable_name))
+	{
+		for (size_t k = 0; k != contact_particles_.size(); ++k)
+		{
+			shell_variable_.push_back(contact_particles_[k]->template getVariableByName<VariableType>(variable_name));
+		}
+	}
+	//=================================================================================================//
+	template <typename VariableType>
+	void DampingPairwiseFromShell<VariableType>::interaction(size_t index_i, Real dt)
+	{
+		Real Vol_i = Vol_[index_i];
+		Real mass_i = mass_[index_i];
+		VariableType &variable_i = variable_[index_i];
+
+		std::array<Real, MaximumNeighborhoodSize> parameter_b;
+
+		/** Contact interaction. */
+		for (size_t k = 0; k < contact_configuration_.size(); ++k)
+		{
+			StdLargeVec<VariableType> &variable_k = *(shell_variable_[k]);
+			Neighborhood &contact_neighborhood = (*contact_configuration_[k])[index_i];
+			// forward sweep
+			for (size_t n = 0; n != contact_neighborhood.current_size_; ++n)
+			{
+				size_t index_j = contact_neighborhood.j_[n];
+
+				parameter_b[n] = contact_particles_[k]->DegeneratedSpacing(index_j) * eta_ * contact_neighborhood.dW_ijV_j_[n] * Vol_i * dt / contact_neighborhood.r_ij_[n];
+
+				// only update particle i
+				variable_[index_i] += parameter_b[n] * (variable_i - variable_k[index_j]) / (mass_i - 2.0 * parameter_b[n]);
+			}
+			// backward sweep
+			for (size_t n = contact_neighborhood.current_size_; n != 0; --n)
+			{
+				size_t index_j = contact_neighborhood.j_[n - 1];
+
+				// only update particle i
+				variable_[index_i] += parameter_b[n - 1] * (variable_i - variable_k[index_j]) / (mass_i - 2.0 * parameter_b[n - 1]);
+			}
+		}
+	}
+	//=================================================================================================//
 	template <class DampingAlgorithmType>
 	template <typename... ConstructorArgs>
 	DampingWithRandomChoice<DampingAlgorithmType>::

SPHINXsys/src/shared/particle_dynamics/fluid_dynamics/all_fluid_dynamics.h

@@ -40,3 +40,5 @@
 #include "fluid_boundary.h"
 #include "all_eulerian_compressible_fluid_dynamics.h"
 #include "all_eulerian_weakly_compressible_fluid_dynamics.h"
+#include "fluid_shell_interaction.h"
+#include "fluid_shell_interaction.hpp"

SPHINXsys/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_complex.cpp

@@ -36,7 +36,7 @@ namespace SPH
 					Vecd nablaW_ijV_j = contact_neighborhood.dW_ijV_j_[n] * contact_neighborhood.e_ij_[n];
 
 					// acceleration for transport velocity
-					acceleration_trans -= 2.0 * nablaW_ijV_j;
+					acceleration_trans -= 2.0 * nablaW_ijV_j * this->contact_particles_[k]->DegeneratedSpacing(index_j);
 				}
 			}
 
@@ -59,7 +59,7 @@ namespace SPH
 					Vecd nablaW_ijV_j = contact_neighborhood.dW_ijV_j_[n] * contact_neighborhood.e_ij_[n];
 
 					// acceleration for transport velocity
-					acceleration_trans -= 2.0 * nablaW_ijV_j;
+					acceleration_trans -= 2.0 * nablaW_ijV_j * this->contact_particles_[k]->DegeneratedSpacing(index_j);
 				}
 			}
 

SPHINXsys/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_complex.hpp

@@ -88,7 +88,8 @@ namespace SPH
 				Neighborhood &contact_neighborhood = (*this->contact_configuration_[k])[index_i];
 				for (size_t n = 0; n != contact_neighborhood.current_size_; ++n)
 				{
-					sigma += contact_neighborhood.W_ij_[n] * contact_inv_rho0_k * contact_mass_k[contact_neighborhood.j_[n]];
+					Real mass = this->contact_particles_[k]->ParticleMass(contact_neighborhood.j_[n]);
+					sigma += contact_neighborhood.W_ij_[n] * contact_inv_rho0_k * mass;
 				}
 			}
 			return sigma;

SPHINXsys/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_inner.cpp

@@ -14,6 +14,7 @@ namespace SPH
 		//=================================================================================================//
 		BaseDensitySummationInner::BaseDensitySummationInner(BaseInnerRelation &inner_relation)
 			: LocalDynamics(inner_relation.sph_body_), FluidDataInner(inner_relation),
+			  spacing_ref_(sph_body_.sph_adaptation_->ReferenceSpacing()),
 			  rho_(particles_->rho_), rho_sum_(particles_->rho_sum_), mass_(particles_->mass_),
 			  rho0_(sph_body_.base_material_->ReferenceDensity()) {}
 		//=================================================================================================//

SPHINXsys/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_inner.h

@@ -75,6 +75,7 @@ namespace SPH
 
 		protected:
 			Real rho0_;
+			Real spacing_ref_;
 			StdLargeVec<Real> &rho_, &rho_sum_, &mass_;
 		};
 

SPHINXsys/src/shared/particle_dynamics/fluid_dynamics/fluid_shell_interaction.cpp

@@ -0,0 +1,12 @@
+#include "fluid_shell_interaction.h"
+
+//=================================================================================================//
+namespace SPH
+{
+	//=================================================================================================//
+	namespace fluid_dynamics
+	{
+		//=================================================================================================//
+    }
+//=================================================================================================//
+}

SPHINXsys/src/shared/particle_dynamics/fluid_dynamics/fluid_shell_interaction.h

@@ -0,0 +1,136 @@
+/* -------------------------------------------------------------------------*
+ *								SPHinXsys									*
+ * -------------------------------------------------------------------------*
+ * SPHinXsys (pronunciation: s'finksis) is an acronym from Smoothed Particle*
+ * Hydrodynamics for industrial compleX systems. It provides C++ APIs for	*
+ * physical accurate simulation and aims to model coupled industrial dynamic*
+ * systems including fluid, solid, multi-body dynamics and beyond with SPH	*
+ * (smoothed particle hydrodynamics), a meshless computational method using	*
+ * particle discretization.													*
+ *																			*
+ * SPHinXsys is partially funded by German Research Foundation				*
+ * (Deutsche Forschungsgemeinschaft) DFG HU1527/6-1, HU1527/10-1,			*
+ *  HU1527/12-1 and HU1527/12-4												*
+ *                                                                          *
+ * Portions copyright (c) 2017-2022 Technical University of Munich and		*
+ * the authors' affiliations.												*
+ *                                                                          *
+ * Licensed under the Apache License, Version 2.0 (the "License"); you may  *
+ * not use this file except in compliance with the License. You may obtain a*
+ * copy of the License at http://www.apache.org/licenses/LICENSE-2.0.       *
+ *                                                                          *
+ * ------------------------------------------------------------------------*/
+/**
+* @file     fluid_shell_interaction.h
+* @brief    Here, we define the algorithm classes for the interaction between fluid and
+*			thin structure, plate and shell.
+* @author	Chi ZHang and Xiangyu Hu
+*/
+
+
+#ifndef FLUID_SHELL_INTERACTION_H
+#define FLUID_SHELL_INTERACTION_H
+
+#include "fluid_dynamics_complex.h"
+#include "fluid_dynamics_complex.hpp"
+
+namespace SPH
+{
+	namespace fluid_dynamics
+	{		
+		typedef DataDelegateContact<FluidParticles, ShellParticles> FluidShellContactData;
+		typedef DataDelegateContact<FluidParticles, ShellParticles, DataDelegateEmptyBase> FluidShellData;
+
+		/**
+		* @class InteractionWithShell
+		* @brief Abstract base class for general fluid-shell interaction model. 
+		* @note  Here, 
+		*/
+		template<class BaseIntegrationType>
+		class InteractionWithShell : public BaseIntegrationType, public FluidShellData
+		{
+		public:
+			template <class BaseBodyRelationType, typename... Args>
+			InteractionWithShell(BaseContactRelation &contact_relation, BaseBodyRelationType &base_body_relation, Args &&...args);
+
+			template <typename... Args>
+			InteractionWithShell(ComplexRelation &fluid_shell_relation, Args &&...args)
+				: InteractionWithShell(fluid_shell_relation.getContactRelation(), fluid_shell_relation.getInnerRelation(), std::forward<Args>(args)...) 
+			{}
+
+			virtual ~InteractionWithShell(){};
+
+		protected:
+			Real spacing_ref_;
+			StdVec<StdLargeVec<Real> *> shell_thickness_;
+			StdVec<StdLargeVec<Vecd> *> shell_n_, shell_vel_ave_, shell_acc_ave_;
+		};
+
+		/**
+		 * @class BaseFluidShellIntegration1stHalf
+		 * @brief  template class for fluid-shell pressure relaxation scheme
+		 */
+		template <class Integration1stHalfType>
+		class BaseFluidShellIntegration1stHalf : public InteractionWithShell<Integration1stHalfType>
+		{
+		public:
+			template <typename... Args>
+			explicit BaseFluidShellIntegration1stHalf(Args &&...args) 
+				: InteractionWithShell<Integration1stHalfType>(std::forward<Args>(args)...){};
+			virtual ~BaseFluidShellIntegration1stHalf(){};
+			void interaction(size_t index_i, Real dt = 0.0);
+
+		protected:
+			virtual Vecd computeNonConservativeAcceleration(size_t index_i) override;
+		};
+
+		using FluidShellIntegration1stHalf = BaseFluidShellIntegration1stHalf<Integration1stHalf>;
+		using FluidShellIntegration1stHalfRiemann = BaseFluidShellIntegration1stHalf<Integration1stHalfRiemann>;
+
+		using FluidShellandWallIntegration1stHalf = BaseIntegration1stHalfWithWall<FluidShellIntegration1stHalf>;
+		using FluidShellandWallIntegration1stHalfRiemann = BaseIntegration1stHalfWithWall<FluidShellIntegration1stHalfRiemann>;
+		using ExtendFluidShellandWallIntegration1stHalfRiemann = BaseExtendIntegration1stHalfWithWall<FluidShellIntegration1stHalfRiemann>;
+
+		/**
+		 * @class BaseFluidShellIntegration2ndHalf
+		 * @brief  template class pressure relaxation scheme with wall boundary
+		 */
+		template <class Integration2ndHalfType>
+		class BaseFluidShellIntegration2ndHalf : public InteractionWithShell<Integration2ndHalfType>
+		{
+		public:
+			template <typename... Args>
+			explicit BaseFluidShellIntegration2ndHalf(Args &&...args) 
+				: InteractionWithShell<Integration2ndHalfType>(std::forward<Args>(args)...)
+			{};
+			virtual ~BaseFluidShellIntegration2ndHalf(){};
+			void interaction(size_t index_i, Real dt = 0.0);
+		};
+		using FluidShellIntegration2ndHalf = BaseFluidShellIntegration2ndHalf<Integration2ndHalf>;
+		using FluidShellIntegration2ndHalfRiemann = BaseFluidShellIntegration2ndHalf<Integration2ndHalfRiemann>;
+
+		using FluidShellandWallIntegration2ndHalf = BaseIntegration2ndHalfWithWall<FluidShellIntegration2ndHalf>;
+		using FluidShellandWallIntegration2ndHalfRiemann = BaseIntegration2ndHalfWithWall<FluidShellIntegration2ndHalfRiemann>;
+
+		/**
+		 * @class ViscousWithWall
+		 * @brief  template class viscous acceleration with wall boundary
+		 */
+		template <class ViscousAccelerationInnerType>
+		class BaseViscousAccelerationWithShell : public InteractionWithShell<ViscousAccelerationInnerType>
+		{
+		public:
+			template <typename... Args>
+			BaseViscousAccelerationWithShell(Args &&...args)
+				: InteractionWithShell<ViscousAccelerationInnerType>(std::forward<Args>(args)...)
+			{};
+			virtual ~BaseViscousAccelerationWithShell(){};
+			void interaction(size_t index_i, Real dt = 0.0);
+		};
+
+		using ViscousAccelerationWithShell = BaseViscousAccelerationWithShell<ViscousAccelerationInner>;
+		using ViscousAccelerationWithShellandWall = BaseViscousAccelerationWithWall<ViscousAccelerationWithShell>;
+
+    }
+}
+#endif //FLUID_SHELL_INTERACTION_H