1 // =============================================================================
2 // PROJECT CHRONO - http://projectchrono.org
3 //
4 // Copyright (c) 2014 projectchrono.org
5 // All rights reserved.
6 //
7 // Use of this source code is governed by a BSD-style license that can be found
8 // in the LICENSE file at the top level of the distribution and at
9 // http://projectchrono.org/license-chrono.txt.
10 //
11 // =============================================================================
12 // Author: Milad Rakhsha, Wei Hu
13 // =============================================================================
14
15 // General Includes
16 #include <assert.h>
17 #include <stdlib.h>
18 #include <ctime>
19
20 // Chrono includes
21 #include "chrono/physics/ChSystemSMC.h"
22 #include "chrono/utils/ChUtilsCreators.h"
23 #include "chrono/utils/ChUtilsGenerators.h"
24 #include "chrono/utils/ChUtilsGeometry.h"
25
26 // Chrono fsi includes
27 #include "chrono_fsi/ChSystemFsi.h"
28
29 // Chrono namespaces
30 using namespace chrono;
31 using namespace chrono::collision;
32 using namespace chrono::fsi;
33
34 // Output directories and settings
35 const std::string out_dir = GetChronoOutputPath() + "FSI_POISEUILLE_FLOW/";
36 std::string demo_dir;
37
38 // Save data as csv files to see the results off-line using Paraview
39 bool save_output = true;
40
41 // Dimension of the space domain
42 double bxDim = 0.2;
43 double byDim = 0.1;
44 double bzDim = 0.2;
45 // Dimension of the fluid domain
46 double fxDim = 0.2;
47 double fyDim = 0.1;
48 double fzDim = 0.2;
49
ShowUsage()50 void ShowUsage() {
51 std::cout << "usage: ./demo_FSI_Poiseuille_flow <json_file>" << std::endl;
52 }
53
54 //------------------------------------------------------------------
55 // Function to save the paraview files
56 //------------------------------------------------------------------
SaveParaViewFilesMBD(ChSystemFsi & myFsiSystem,ChSystemSMC & mphysicalSystem,std::shared_ptr<fsi::SimParams> paramsH,int next_frame,double mTime)57 void SaveParaViewFilesMBD(ChSystemFsi& myFsiSystem,
58 ChSystemSMC& mphysicalSystem,
59 std::shared_ptr<fsi::SimParams> paramsH,
60 int next_frame,
61 double mTime) {
62 // Simulation time between two output frames
63 double frame_time = 1.0 / paramsH->out_fps;
64
65 // Output data to files
66 if (save_output && std::abs(mTime - (next_frame)*frame_time) < 1e-6) {
67 myFsiSystem.PrintParticleToFile(demo_dir);
68
69 std::cout << "\n--------------------------------\n" << std::endl;
70 std::cout << "------------ Output Frame: " << next_frame << std::endl;
71 std::cout << "------------ Sim Time: " << mTime << " (s)\n" << std::endl;
72 std::cout << "--------------------------------\n" << std::endl;
73 }
74 }
75
76 //------------------------------------------------------------------
77 // Create the objects of the MBD system. Rigid bodies, and if FSI,
78 // their BCE representation are created and added to the systems
79 //------------------------------------------------------------------
CreateSolidPhase(ChSystemSMC & mphysicalSystem,ChSystemFsi & myFsiSystem,std::shared_ptr<fsi::SimParams> paramsH)80 void CreateSolidPhase(ChSystemSMC& mphysicalSystem,
81 ChSystemFsi& myFsiSystem,
82 std::shared_ptr<fsi::SimParams> paramsH) {
83 // Set common material Properties
84 auto mysurfmaterial = chrono_types::make_shared<ChMaterialSurfaceSMC>();
85 mysurfmaterial->SetYoungModulus(6e4);
86 mysurfmaterial->SetFriction(0.3f);
87 mysurfmaterial->SetRestitution(0.2f);
88 mysurfmaterial->SetAdhesion(0);
89
90 // General setting of ground body
91 auto ground = chrono_types::make_shared<ChBody>();
92 ground->SetIdentifier(-1);
93 ground->SetBodyFixed(true);
94 ground->SetCollide(true);
95 ground->GetCollisionModel()->ClearModel();
96
97 // Create the geometry of the boundaries
98 auto initSpace0 = paramsH->MULT_INITSPACE * paramsH->HSML;
99
100 // Bottom and top wall - size and position
101 ChVector<> sizeBottom(bxDim / 2, byDim / 2 + 0 * initSpace0, 2 * initSpace0);
102 ChVector<> sizeTop = sizeBottom;
103 ChVector<> posBottom(0, 0, -3 * initSpace0);
104 ChVector<> posTop(0, 0, bzDim + 1 * initSpace0);
105
106 // Add the walls into chrono system
107 chrono::utils::AddBoxGeometry(ground.get(), mysurfmaterial, sizeBottom, posBottom, QUNIT, true);
108 chrono::utils::AddBoxGeometry(ground.get(), mysurfmaterial, sizeTop, posTop, QUNIT, true);
109
110 ground->GetCollisionModel()->BuildModel();
111 mphysicalSystem.AddBody(ground);
112
113 // Add BCE particles attached on the walls into FSI system
114 myFsiSystem.AddBceBox(paramsH, ground, posBottom, QUNIT, sizeBottom);
115 myFsiSystem.AddBceBox(paramsH, ground, posTop, QUNIT, sizeBottom);
116
117 }
118
119 // =============================================================================
main(int argc,char * argv[])120 int main(int argc, char* argv[]) {
121 // Create a physics system and an FSI system
122 ChSystemSMC mphysicalSystem;
123 ChSystemFsi myFsiSystem(mphysicalSystem);
124
125 // Get the pointer to the system parameter and use a JSON file to fill it out with the user parameters
126 std::shared_ptr<fsi::SimParams> paramsH = myFsiSystem.GetSimParams();
127
128 // Use the default input file or you may enter your input parameters as a command line argument
129 std::string inputJson = GetChronoDataFile("fsi/input_json/demo_FSI_Poiseuille_flow_Explicit.json");
130 if (argc == 1) {
131 std::cout << "Use the default JSON file" << std::endl;
132 } else if (argc == 2) {
133 std::cout << "Use the specified JSON file" << std::endl;
134 std::string my_inputJson = std::string(argv[1]);
135 inputJson = my_inputJson;
136 } else {
137 ShowUsage();
138 return 1;
139 }
140 myFsiSystem.SetSimParameter(inputJson, paramsH, ChVector<>(bxDim, byDim, bzDim));
141
142 // Dimension of the space domain
143 bxDim = paramsH->boxDimX;
144 byDim = paramsH->boxDimY;
145 bzDim = paramsH->boxDimZ;
146 // Dimension of the fluid domain
147 fxDim = paramsH->fluidDimX;
148 fyDim = paramsH->fluidDimY;
149 fzDim = paramsH->fluidDimZ;
150
151 // Set the time integration type and the linear solver type (only for ISPH)
152 myFsiSystem.SetFluidDynamics(paramsH->fluid_dynamic_type);
153 myFsiSystem.SetFluidSystemLinearSolver(paramsH->LinearSolver);
154
155 // Set the periodic boundary condition (in X and Y direction)
156 auto initSpace0 = paramsH->MULT_INITSPACE * paramsH->HSML;
157 ChVector<> cMin = ChVector<>(-bxDim / 2 - initSpace0 / 2, -byDim / 2 - initSpace0 / 2, -5.0 * initSpace0);
158 ChVector<> cMax = ChVector<>( bxDim / 2 + initSpace0 / 2, byDim / 2 + initSpace0 / 2, bzDim + 5.0 * initSpace0);
159 myFsiSystem.SetBoundaries(cMin, cMax, paramsH);
160
161 // Setup sub doamins for a faster neighbor particle searching
162 myFsiSystem.SetSubDomain(paramsH);
163
164 // Setup the output directory for FSI data
165 myFsiSystem.SetFsiOutputDir(paramsH, demo_dir, out_dir, inputJson);
166
167 // Create Fluid region and discretize with SPH particles
168 ChVector<> boxCenter(0.0, 0.0, fzDim / 2);
169 ChVector<> boxHalfDim(fxDim / 2, fyDim / 2, fzDim / 2);
170
171 // Use a chrono sampler to create a bucket of points
172 chrono::utils::GridSampler<> sampler(initSpace0);
173 chrono::utils::Generator::PointVector points = sampler.SampleBox(boxCenter, boxHalfDim);
174
175 // Add fluid particles from the sampler points to the FSI system
176 size_t numPart = points.size();
177 for (int i = 0; i < numPart; i++) {
178 myFsiSystem.AddSphMarker(points[i], paramsH->rho0, paramsH->BASEPRES, paramsH->mu0, paramsH->HSML, -1);
179 }
180 myFsiSystem.AddRefArray(0, (int)numPart, -1, -1);
181
182 // Create Solid region and attach BCE SPH particles
183 CreateSolidPhase(mphysicalSystem, myFsiSystem, paramsH);
184
185 // Construction of the FSI system must be finalized before running
186 myFsiSystem.Finalize();
187
188 // Start the simulation
189 double time = 0;
190 int stepEnd = int(paramsH->tFinal / paramsH->dT);
191 double TIMING_sta = clock();
192 for (int tStep = 0; tStep < stepEnd + 1; tStep++) {
193 printf("\nstep : %d, time= : %f (s) \n", tStep, time);
194 double frame_time = 1.0 / paramsH->out_fps;
195 int next_frame = (int)floor((time + 1e-6) / frame_time) + 1;
196
197 // Call the FSI solver
198 myFsiSystem.DoStepDynamics_FSI();
199 time += paramsH->dT;
200
201 // Save data of the simulation
202 SaveParaViewFilesMBD(myFsiSystem, mphysicalSystem, paramsH, next_frame, time);
203 }
204
205 // Total computational cost
206 double TIMING_end = (clock() - TIMING_sta) / (double)CLOCKS_PER_SEC;
207 printf("\nSimulation Finished in %f (s)\n", TIMING_end);
208
209 return 0;
210 }
211