Merge pull request #547 from Yashmandaokar/Incompressible-taper-chann…

…el-flow

Add New 3D Incompressible Channel Flow Test and Unstructured Mesh Reader

.github/workflows/ci.yml

@@ -540,4 +540,4 @@ jobs:
         if: ${{ steps.fourth-try.outcome == 'failure' }}
         run: |
           cd build 
-          ctest --rerun-failed --output-on-failure
+          ctest --rerun-failed --output-on-failure

src/shared/bodies/complex_bodies/mesh_helper.h

@@ -0,0 +1,86 @@
+/* ------------------------------------------------------------------------- *
+ *                                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-2023 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 	MeshHelper.h
+  * @brief 	Here, we define the functions used for reading the mesh data.
+  * @author	Yash Mandaokar Zhentong Wang and Xiangyu Hu
+  */
+#ifndef MESHTYPE_H
+#define MESHTYPE_H
+
+#include "unstructured_mesh.h"
+using namespace std;
+namespace SPH
+{
+
+    class MeshFileHelpers 
+    {
+        public:
+
+        MeshFileHelpers()
+        {
+
+        };
+        virtual ~MeshFileHelpers(){};
+
+        static void meshDimension(ifstream& mesh_file, size_t& dimension, string& text_line);
+        static void numberofNodes(ifstream& mesh_file, size_t& number_of_points, string& text_line);
+        static void nodeCoordinates(ifstream& mesh_file, StdLargeVec<Vecd>& node_coordinates_, string& text_line,size_t& dimension);
+        static void numberofElements(ifstream& mesh_file, size_t& number_of_elements, string& text_line);
+        static void dataStruct(vector<vector<vector<size_t>>>& mesh_topology_, StdLargeVec<StdVec<size_t>>& elements_nodes_connection_,
+                                size_t number_of_elements, size_t mesh_type, size_t dimension);
+        static size_t findBoundaryType(string& text_line, size_t boundary_type);
+        static Vecd nodeIndex(string& text_line);
+        static Vec2d cellIndex(string& text_line);
+        static void updateElementsNodesConnection(StdLargeVec<StdVec<size_t>>& elements_nodes_connection_, Vecd nodes, Vec2d cells,
+                                                    bool& check_neighbor_cell1, bool& check_neighbor_cell2);
+        static void updateCellLists(vector<vector<vector<size_t>>>& mesh_topology_, StdLargeVec<StdVec<size_t>>& elements_nodes_connection_, Vecd nodes,
+                                    Vec2d cells, bool& check_neighbor_cell1, bool& check_neighbor_cell2, size_t boundary_type);
+        static void updateBoundaryCellLists(vector<vector<vector<size_t>>>& mesh_topology_, StdLargeVec<StdVec<size_t>>& elements_nodes_connection_,
+                                            Vecd nodes, Vec2d cells, bool& check_neighbor_cell1, bool& check_neighbor_cell2, size_t boundary_type);
+        static void cellCenterCoordinates(StdLargeVec<StdVec<size_t>>& elements_nodes_connection_, std::size_t& element,
+                                        StdLargeVec<Vecd>& node_coordinates_, StdLargeVec<Vecd>& elements_center_coordinates_, Vecd& center_coordinate);
+        static void elementVolume(StdLargeVec<StdVec<size_t>>& elements_nodes_connection_, std::size_t& element,
+                                    StdLargeVec<Vecd>& node_coordinates_, StdLargeVec<Real>& elements_volumes_);
+        static void minimumdistance(vector<Real>& all_data_of_distance_between_nodes, StdLargeVec<Real>& elements_volumes_, 
+                                    vector<vector<vector<size_t>>>& mesh_topology_, StdLargeVec<Vecd>& node_coordinates_);
+        static void vtuFileHeader(std::ofstream& out_file);
+        static void vtuFileNodeCoordinates(std::ofstream& out_file, StdLargeVec<Vecd>& nodes_coordinates_, 
+                                    StdLargeVec<StdVec<size_t>>& elements_nodes_connection_, SPHBody& bounds_,Real& rangemax);
+        static void vtuFileInformationKey(std::ofstream& out_file, Real& rangemax);
+        static void vtuFileCellConnectivity(std::ofstream& out_file, StdLargeVec<StdVec<size_t>>& elements_nodes_connection_, StdLargeVec<Vecd>& node_coordinates_);
+        static void vtuFileOffsets(std::ofstream& out_file, StdLargeVec<StdVec<size_t>>& elements_nodes_connection_);
+        static void vtuFileTypeOfCell(std::ofstream& out_file, StdLargeVec<StdVec<size_t>>& elements_nodes_connection_);
+
+        /*Functions for .msh file from FLUENT*/
+        static void numberofNodesFluent(ifstream& mesh_file, size_t& number_of_points, string& text_line);
+        static void numberofElementsFluent(ifstream& mesh_file, size_t& number_of_elements, string& text_line);
+        static void nodeCoordinatesFluent(ifstream& mesh_file, StdLargeVec<Vecd>& node_coordinates_, string& text_line,
+                                            size_t& dimension);
+        static void updateBoundaryCellListsFluent(vector<vector<vector<size_t>>>& mesh_topology_, StdLargeVec<StdVec<size_t>>& elements_nodes_connection_,
+                                            Vecd nodes, Vec2d cells, bool& check_neighbor_cell1, bool& check_neighbor_cell2, size_t boundary_type);
+    };
+
+
+}
+#endif // MESHTYPE_H

src/shared/bodies/complex_bodies/unstructured_mesh.h

@@ -33,6 +33,7 @@
 #include "compressible_fluid.h"
 #include "fluid_body.h"
 #include "io_vtk.h"
+#include <Eigen/Dense>
 using namespace std;
 namespace SPH
 {
@@ -57,11 +58,9 @@ class ANSYSMesh
   protected:
     double min_distance_between_nodes_;
 
-    void readNodeCoordinate(const std::string &text_line, StdLargeVec<Vec2d> &node_coordinates);
-    void readNodeCoordinate(const std::string &text_line, StdLargeVec<Vec3d> &node_coordinates);
     void getDataFromMeshFile(const std::string &full_path);
     void getElementCenterCoordinates();
-    void computeMinimumDistanceBetweenNodes();
+    void getMinimumDistanceBetweenNodes();
 };
 
 /**
@@ -189,9 +188,20 @@ class BodyStatesRecordingInMeshToVtp : public BodyStatesRecording
   protected:
     virtual void writeWithFileName(const std::string &sequence) override;
     StdLargeVec<Vecd> &node_coordinates_;
-    StdLargeVec<StdVec<size_t>> &elements_nodes_connection_;
+    StdLargeVec<StdVec<size_t>>&elements_nodes_connection_;
+};
+class BodyStatesRecordingInMeshToVtu : public BodyStatesRecording
+{
+public:
+    BodyStatesRecordingInMeshToVtu(SPHBody& body, ANSYSMesh& ansys_mesh);
+    virtual ~BodyStatesRecordingInMeshToVtu() {};
+
+protected:
+    virtual void writeWithFileName(const std::string& sequence) override;
+    StdLargeVec<Vecd>& node_coordinates_;
+    StdLargeVec<StdVec<size_t>>& elements_nodes_connection_;
+    SPHBody& bounds_;
 };
-
 //----------------------------------------------------------------------
 //	BoundaryConditionSetupInFVM
 //----------------------------------------------------------------------
@@ -206,6 +216,9 @@ class BoundaryConditionSetupInFVM : public fluid_dynamics::FluidDataInner
     virtual void applyOutletBoundary(size_t ghost_index, size_t index_i){};
     virtual void applyTopBoundary(size_t ghost_index, size_t index_i){};
     virtual void applyFarFieldBoundary(size_t ghost_index){};
+    virtual void applyPressureOutletBC(size_t ghost_index, size_t index_i) {};
+    virtual void applySymmetryBoundary(size_t ghost_index, size_t index_i, Vecd e_ij) {};
+    virtual void applyVelocityInletFlow(size_t ghost_index, size_t index_i) {};
     // Common functionality for resetting boundary conditions
     void resetBoundaryConditions();
 

tests/3d_examples/test_3d_incompressible_channel_flow/CMakeLists.txt

@@ -0,0 +1,24 @@
+STRING(REGEX REPLACE ".*/(.*)" "\\1" CURRENT_FOLDER ${CMAKE_CURRENT_SOURCE_DIR})
+PROJECT("${CURRENT_FOLDER}")
+
+SET(LIBRARY_OUTPUT_PATH ${PROJECT_BINARY_DIR}/lib)
+SET(EXECUTABLE_OUTPUT_PATH "${PROJECT_BINARY_DIR}/bin/")
+SET(BUILD_INPUT_PATH "${EXECUTABLE_OUTPUT_PATH}/input")
+SET(BUILD_RELOAD_PATH "${EXECUTABLE_OUTPUT_PATH}/reload")
+
+file(MAKE_DIRECTORY ${BUILD_INPUT_PATH})
+file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/regression_test_tool/
+        DESTINATION ${BUILD_INPUT_PATH})
+file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/data/Channel_ICEM.msh
+        DESTINATION ${BUILD_INPUT_PATH})
+
+add_executable(${PROJECT_NAME})
+aux_source_directory(. DIR_SRCS)
+target_sources(${PROJECT_NAME} PRIVATE ${DIR_SRCS})
+target_link_libraries(${PROJECT_NAME} sphinxsys_3d)
+set_target_properties(${PROJECT_NAME} PROPERTIES VS_DEBUGGER_WORKING_DIRECTORY "${EXECUTABLE_OUTPUT_PATH}")
+
+add_test(NAME ${PROJECT_NAME}
+        COMMAND ${PROJECT_NAME} --state_recording=${TEST_STATE_RECORDING}
+        WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH})
+set_tests_properties(${PROJECT_NAME} PROPERTIES LABELS "FVM, Eulerian")

tests/3d_examples/test_3d_incompressible_channel_flow/common_weakly_compressible_FVM_classes.cpp

@@ -0,0 +1,53 @@
+
+#include "common_weakly_compressible_FVM_classes.h"
+
+using namespace std;
+namespace SPH
+{
+namespace fluid_dynamics
+{
+//=================================================================================================//
+WCAcousticTimeStepSizeInFVM::WCAcousticTimeStepSizeInFVM(SPHBody &sph_body, Real min_distance_between_nodes, Real acousticCFL)
+    : AcousticTimeStepSize(sph_body), rho_(particles_->rho_), p_(*particles_->getVariableByName<Real>("Pressure")),
+      vel_(particles_->vel_), fluid_(DynamicCast<WeaklyCompressibleFluid>(this, particles_->getBaseMaterial())),
+      min_distance_between_nodes_(min_distance_between_nodes), acousticCFL_(acousticCFL){};
+//=================================================================================================//
+Real WCAcousticTimeStepSizeInFVM::outputResult(Real reduced_value)
+{
+    // I chose a time-step size according to Eulerian method
+    return acousticCFL_ / Dimensions * min_distance_between_nodes_ / (reduced_value + TinyReal);
+}
+//=================================================================================================//
+BaseForceFromFluidInFVM::BaseForceFromFluidInFVM(BaseInnerRelation &inner_relation)
+    : LocalDynamics(inner_relation.getSPHBody()), fluid_dynamics::FluidDataInner(inner_relation), Vol_(particles_->Vol_){};
+//=================================================================================================//
+ViscousForceFromFluidInFVM::ViscousForceFromFluidInFVM(BaseInnerRelation &inner_relation, vector<vector<size_t>> each_boundary_type_contact_real_index)
+    : BaseForceFromFluidInFVM(inner_relation), fluid_(DynamicCast<WeaklyCompressibleFluid>(this, particles_->getBaseMaterial())),
+      vel_(particles_->vel_), mu_(fluid_.ReferenceViscosity()), each_boundary_type_contact_real_index_(each_boundary_type_contact_real_index)
+{
+    particles_->registerVariable(force_from_fluid_, "ViscousForceOnSolid");
+};
+//=================================================================================================//
+void ViscousForceFromFluidInFVM::interaction(size_t index_i, Real dt)
+{
+    for (size_t real_particle_num = 0; real_particle_num != each_boundary_type_contact_real_index_[3].size(); ++real_particle_num)
+    {
+        Vecd force = Vecd::Zero();
+        const Vecd &vel_i = vel_[index_i];
+        if (index_i == each_boundary_type_contact_real_index_[3][real_particle_num])
+        {
+            Real Vol_i = Vol_[index_i];
+            const Neighborhood &inner_neighborhood = inner_configuration_[index_i];
+
+
+            Vecd vel_j = -vel_i;
+            size_t index_j = inner_neighborhood.j_[2];
+            Vecd vel_derivative = (vel_j - vel_i) / (inner_neighborhood.r_ij_[2] + TinyReal);
+            force += 2.0 * mu_ * vel_derivative * Vol_i * inner_neighborhood.dW_ij_[2] * Vol_[index_j];
+            force_from_fluid_[index_i] = force;
+        }
+    }
+}
+//=================================================================================================//
+} // namespace fluid_dynamics
+} // namespace SPH

tests/3d_examples/test_3d_incompressible_channel_flow/common_weakly_compressible_FVM_classes.h

@@ -0,0 +1,140 @@
+/* ------------------------------------------------------------------------- *
+ *                                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-2023 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 	common_weakly_compressible_FVM_classes_3d.h
+ * @brief 	Here, we define the common weakly compressible classes for fluid dynamics in FVM.
+ * @author	Zhentong Wang and Xiangyu Hu
+ */
+
+#ifndef COMMON_WEAKLY_COMPRESSIBLE_FVM_CLASSES_3D_H
+#define COMMON_WEAKLY_COMPRESSIBLE_FVM_CLASSES_3D_H
+#include "unstructured_mesh.h"
+
+namespace SPH
+{
+namespace fluid_dynamics
+{
+/**
+ * @class WCAcousticTimeStepSizeInFVM
+ * @brief Computing the acoustic time step size
+ */
+class WCAcousticTimeStepSizeInFVM : public fluid_dynamics::AcousticTimeStepSize
+{
+  protected:
+    StdLargeVec<Real> &rho_, &p_;
+    StdLargeVec<Vecd> &vel_;
+    Fluid &fluid_;
+    Real min_distance_between_nodes_;
+
+  public:
+    explicit WCAcousticTimeStepSizeInFVM(SPHBody &sph_body, Real min_distance_between_nodes, Real acousticCFL = 0.6);
+    virtual ~WCAcousticTimeStepSizeInFVM(){};
+    virtual Real outputResult(Real reduced_value) override;
+    Real acousticCFL_;
+};
+
+/**
+ * @class BaseForceFromFluidInFVM
+ * @brief Base class for computing the forces from the fluid.
+ * Note that In FVM , we need FluidDataInner class to calculate force between solid and fluid.
+ */
+class BaseForceFromFluidInFVM : public LocalDynamics, public fluid_dynamics::FluidDataInner
+{
+  public:
+    explicit BaseForceFromFluidInFVM(BaseInnerRelation &inner_relation);
+    virtual ~BaseForceFromFluidInFVM(){};
+    StdLargeVec<Vecd> &getForceFromFluid() { return force_from_fluid_; };
+
+  protected:
+    StdLargeVec<Real> &Vol_;
+    StdLargeVec<Vecd> force_from_fluid_;
+};
+
+/**
+ * @class ViscousForceFromFluidInFVM
+ * @brief Computing the viscous force from the fluid
+ */
+class ViscousForceFromFluidInFVM : public BaseForceFromFluidInFVM
+{
+  public:
+    explicit ViscousForceFromFluidInFVM(BaseInnerRelation &inner_relation, vector<vector<size_t>> each_boundary_type_contact_real_index);
+    virtual ~ViscousForceFromFluidInFVM(){};
+    void interaction(size_t index_i, Real dt = 0.0);
+
+  protected:
+    Fluid &fluid_;
+    StdLargeVec<Vecd> &vel_;
+    Real mu_;
+    vector<vector<size_t>> each_boundary_type_contact_real_index_;
+};
+
+/**
+ * @class PressureForceFromFluidInFVM
+ * @brief Template class fro computing the pressure force from the fluid with different Riemann solvers in FVM.
+ * The pressure force is added on the viscous force of the latter is computed.
+ * time step size compared to the fluid dynamics
+ */
+template <class EulerianIntegration2ndHalfType>
+class PressureForceFromFluidInFVM : public BaseForceFromFluidInFVM
+{
+    using RiemannSolverType = typename EulerianIntegration2ndHalfType::RiemannSolver;
+
+  public:
+    explicit PressureForceFromFluidInFVM(BaseInnerRelation &inner_relation, vector<vector<size_t>> each_boundary_type_contact_real_index)
+        : BaseForceFromFluidInFVM(inner_relation), fluid_(DynamicCast<WeaklyCompressibleFluid>(this, particles_->getBaseMaterial())), vel_(particles_->vel_),
+          p_(*particles_->getVariableByName<Real>("Pressure")), rho_(particles_->rho_), riemann_solver_(fluid_, fluid_),
+          each_boundary_type_contact_real_index_(each_boundary_type_contact_real_index)
+    {
+        particles_->registerVariable(force_from_fluid_, "PressureForceOnSolid");
+    };
+    Fluid &fluid_;
+    StdLargeVec<Vecd> &vel_;
+    StdLargeVec<Real> &p_, &rho_;
+    RiemannSolverType riemann_solver_;
+    vector<vector<size_t>> each_boundary_type_contact_real_index_;
+    virtual ~PressureForceFromFluidInFVM(){};
+
+    void interaction(size_t index_i, Real dt = 0.0)
+    {
+        for (size_t real_particle_num = 0; real_particle_num != each_boundary_type_contact_real_index_[3].size(); ++real_particle_num)
+        {
+            Vecd force = Vecd::Zero();
+            if (index_i == each_boundary_type_contact_real_index_[3][real_particle_num])
+            {
+                Real Vol_i = Vol_[index_i];
+                FluidStateIn state_i(rho_[index_i], vel_[index_i], p_[index_i]);
+                const Neighborhood &inner_neighborhood = inner_configuration_[index_i];
+                size_t index_j = inner_neighborhood.j_[2];
+                Vecd e_ij = inner_neighborhood.e_ij_[2];
+                FluidStateIn state_j(rho_[index_j], vel_[index_j], p_[index_j]);
+                FluidStateOut interface_state = riemann_solver_.InterfaceState(state_i, state_j, e_ij);
+                force -= 2.0 * (-e_ij) * interface_state.p_ * Vol_i * inner_neighborhood.dW_ij_[2] * Vol_[index_j];
+                force_from_fluid_[index_i] = force;
+            }
+        }
+    };
+};
+
+} // namespace fluid_dynamics
+} // namespace SPH
+#endif //COMMON_WEAKLY_COMPRESSIBLE_FVM_CLASSES_3D_H