Hi, I had some questions that I was hoping you could help me with. I have been studying the FSI demos to try to better understand the FSI module. I have mostly been using the DEM engine, so I have noticed some significant differences between the DEM module and the FSI module. I am trying to duplicate a simulation that I have done in the DEM ( https://drive.google.com/file/d/1L9Br1Vf9hsFkjXqbhAZGNpgVep3p5YK_/view?usp=drive_link). In this simulation, I first initialize a bed with specific materials' properties (Young's modulus, possession ratio, static friction, rolling friction, coefficient of restitution, and cohesion energy density (the cohesion module is implemented by me)). Then I settle the bed. After the bed is settled, I drop the screw in a drop phase. After the screw is settled on the bed, I start rotating the screw at 1 rad/s while restricting its movement in the x direction and applying a down force on it.
To duplicate this simulation, I started with the DEMO_FSI_SingleWheelTest where I started switching the wheel geometry to my geometry and changing the bed size to what is in the DEM simulation( https://drive.google.com/file/d/18KhOFfXNEIlgcoLTvxTi27XL5soYB4yz/view?usp=drive_link - FSI_ScrewTest.cpp). The simulation looks okay but I got stuck on many things after. I was wondering if you could help me with this stuff: 1- the bed looked too rigid. I have tried to figure out what to change to get a similar bed to what is in my DEM simulation but I could not. My bed in the DEM sim is generated using a PD sampler and then is allowed to settle. I have noticed that there is no settling phase in the FSI simulation. Is this normal for FSI simulations? What is controlling the rigidity of the bed in this simulation? 2- Bed properties? I have noticed that you get some of the parameters from a JSON file and some of them are defined in the simulation. However, it seems that FSI simulations use some different properties for the bed than what is used in DEM. Is it possible to define the bed using ONLY the material properties that I used in my DEM simulation (mentioned above)? Also, I could not find where the particle radius is defined in the simulation and was wondering if the kernelLength parameter is the same as the particle radius. Finally, I would like to set all my simulation parameters in my .cpp file instead of using a JSON file and was wondering about the appropriate way to set such parameters (such as Young's modulus, Can I have a code example?). 3- I was wondering about the best way to restrict motion in some of the directions. Also, what is the best way to apply forces to the mesh? Thank you so much for your help in advance, -- You received this message because you are subscribed to the Google Groups "ProjectChrono" group. To unsubscribe from this group and stop receiving emails from it, send an email to [email protected]. To view this discussion on the web visit https://groups.google.com/d/msgid/projectchrono/2fceb906-5ad8-41ee-addb-bf81c7e75958n%40googlegroups.com.
// ============================================================================= // PROJECT CHRONO - http://projectchrono.org // // Copyright (c) 2014 projectchrono.org // All rights reserved. // // Use of this source code is governed by a BSD-style license that can be found // in the LICENSE file at the top level of the distribution and at // http://projectchrono.org/license-chrono.txt. // // ============================================================================= // Author: Wei Hu, Radu Serban // ============================================================================= #include <cassert> #include <cstdlib> #include <iostream> #include <vector> #include <string> #include <sstream> #include <cmath> #include "chrono/ChConfig.h" #include "chrono/physics/ChSystemSMC.h" #include "chrono/physics/ChBody.h" #include "chrono/physics/ChInertiaUtils.h" #include "chrono/physics/ChLinkMotorRotationAngle.h" #include "chrono/utils/ChUtilsGeometry.h" #include "chrono/utils/ChUtilsInputOutput.h" #include "chrono/assets/ChTriangleMeshShape.h" #include "chrono/core/ChTimer.h" #include "chrono_fsi/ChSystemFsi.h" #include "chrono_fsi/ChVisualizationFsi.h" #include "chrono_thirdparty/filesystem/path.h" // Chrono namespaces using namespace chrono; using namespace chrono::fsi; using namespace chrono::geometry; // Physical properties of terrain particles double iniSpacing = 0.01; double kernelLength = 0.01;//what is kernel length? double density = 1700.0; // Dimension of the terrain container double smalldis = 1.0e-9; //why do you add small distance? double bxDim = 0.38 + smalldis; double byDim = 0.8 + smalldis; double bzDim = 0.25 + smalldis; // Size of the wheel double wheel_radius = 0.1;//?? double wheel_slip = 0.0; double wheel_AngVel = 1.0; double total_mass = 1.143; std::string wheel_obj = "vehicle/hmmwv/screw17.obj"; // Initial Position of wheel ChVector<> wheel_IniPos(0.0 , 0.0, wheel_radius + bzDim + iniSpacing); ChVector<> wheel_IniVel(0.0, 0.0, 0.0); // Simulation time and stepsize double total_time = 10.0; double dT = 2.5e-4; // linear actuator and angular actuator auto actuator = chrono_types::make_shared<ChLinkLinActuator>(); auto motor = chrono_types::make_shared<ChLinkMotorRotationAngle>(); // Save data as csv files to see the results off-line using Paraview bool output = true; int out_fps = 20; // Output directories and settings const std::string out_dir = GetChronoOutputPath() + "FSI_wheel_Test/"; // Enable/disable run-time visualization (if Chrono::OpenGL is available) bool render = true; float render_fps = 100; // Verbose terminal output bool verbose_fsi = true; bool verbose = true; //------------------------------------------------------------------ // Function to save wheel to Paraview VTK files //------------------------------------------------------------------ void WritewheelVTK(const std::string& filename, ChTriangleMeshConnected& mesh, const ChFrame<>& frame) { std::ofstream outf; outf.open(filename); outf << "# vtk DataFile Version 2.0" << std::endl; outf << "VTK from simulation" << std::endl; outf << "ASCII" << std::endl; outf << "DATASET UNSTRUCTURED_GRID" << std::endl; outf << "POINTS " << mesh.getCoordsVertices().size() << " " << "float" << std::endl; for (auto& v : mesh.getCoordsVertices()) { auto w = frame.TransformPointLocalToParent(v); outf << w.x() << " " << w.y() << " " << w.z() << std::endl; } auto nf = mesh.getIndicesVertexes().size(); outf << "CELLS " << nf << " " << 4 * nf << std::endl; for (auto& f : mesh.getIndicesVertexes()) { outf << "3 " << f.x() << " " << f.y() << " " << f.z() << std::endl; } outf << "CELL_TYPES " << nf << std::endl; for (int i = 0; i < nf; i++) { outf << "5 " << std::endl; } outf.close(); } //------------------------------------------------------------------ // Create the objects of the MBD system. Rigid bodies, and if FSI, // their BCE representation are created and added to the systems //------------------------------------------------------------------ void CreateSolidPhase(ChSystemSMC& sysMBS, ChSystemFsi& sysFSI) { // Common contact material auto cmaterial = chrono_types::make_shared<ChMaterialSurfaceSMC>(); cmaterial->SetYoungModulus(100e9); cmaterial->SetFriction(0.75f); cmaterial->SetRestitution(0.475f); cmaterial->SetAdhesion(0); // Create a container -- always FIRST body in the system auto ground = chrono_types::make_shared<ChBody>(); ground->SetIdentifier(-1); ground->SetBodyFixed(true); sysMBS.AddBody(ground); ground->GetCollisionModel()->ClearModel(); chrono::utils::AddBoxContainer(ground, cmaterial, ChFrame<>(), ChVector<>(bxDim, byDim, bzDim), 0.1, ChVector<int>(0, 0, -1), false); ground->GetCollisionModel()->BuildModel(); ground->SetCollide(true); // Add BCE particles attached on the walls into FSI system sysFSI.AddContainerBCE(ground, ChFrame<>(), ChVector<>(bxDim, byDim, 2 * bzDim), ChVector<int>(2, 0, -1)); // Create the wheel -- always SECOND body in the system auto trimesh = chrono_types::make_shared<ChTriangleMeshConnected>(); double scale_ratio = 1.0; trimesh->LoadWavefrontMesh(GetChronoDataFile(wheel_obj), false, true); trimesh->Transform(ChVector<>(0, 0, 0), ChMatrix33<>(scale_ratio)); // scale to a different size trimesh->RepairDuplicateVertexes(1e-9); // if meshes are not watertight // Compute mass inertia from mesh double mmass; double mdensity = density; ChVector<> mcog; ChMatrix33<> minertia; trimesh->ComputeMassProperties(true, mmass, mcog, minertia); ChMatrix33<> principal_inertia_rot; ChVector<> principal_I; ChInertiaUtils::PrincipalInertia(minertia, principal_I, principal_inertia_rot); mcog = ChVector<>(0.0, 0.0, 0.0); // Set the abs orientation, position and velocity auto wheel = chrono_types::make_shared<ChBodyAuxRef>(); ChQuaternion<> wheel_Rot = Q_from_Euler123(ChVector<double>(0, 0, 0)); // Set the COG coordinates to barycenter, without displacing the REF reference. // Make the COG frame a principal frame. wheel->SetFrame_COG_to_REF(ChFrame<>(mcog, principal_inertia_rot)); // Set inertia wheel->SetMass(total_mass * 1.0 / 2.0); wheel->SetInertiaXX(mdensity * principal_I); wheel->SetPos_dt(wheel_IniVel); wheel->SetWvel_loc(ChVector<>(0.0, 0.0, 0.0)); // set an initial anular velocity (rad/s) // Set the absolute position of the body: wheel->SetFrame_REF_to_abs(ChFrame<>(ChVector<>(wheel_IniPos), ChQuaternion<>(wheel_Rot))); sysMBS.AddBody(wheel); wheel->SetBodyFixed(false); wheel->GetCollisionModel()->ClearModel(); wheel->GetCollisionModel()->AddTriangleMesh(cmaterial, trimesh, false, false, VNULL, ChMatrix33<>(1), 0.005); wheel->GetCollisionModel()->BuildModel(); wheel->SetCollide(false); // Add this body to the FSI system std::vector<ChVector<>> BCE_wheel; sysFSI.CreateMeshPoints(*trimesh, iniSpacing, BCE_wheel); sysFSI.AddPointsBCE(wheel, BCE_wheel, ChFrame<>(), true); sysFSI.AddFsiBody(wheel); // Create the chassis -- always THIRD body in the system auto chassis = chrono_types::make_shared<ChBody>(); chassis->SetMass(total_mass * 1.0 / 2.0); chassis->SetPos(wheel->GetPos()); chassis->SetCollide(false); chassis->SetBodyFixed(false); // Add geometry of the chassis. chassis->GetCollisionModel()->ClearModel(); chrono::utils::AddBoxGeometry(chassis.get(), cmaterial, ChVector<>(0.1, 0.1, 0.1), ChVector<>(0, 0, 0)); chassis->GetCollisionModel()->BuildModel(); sysMBS.AddBody(chassis); // Create the axle -- always FOURTH body in the system auto axle = chrono_types::make_shared<ChBody>(); axle->SetMass(total_mass * 1.0 / 2.0); axle->SetPos(wheel->GetPos()); axle->SetCollide(false); axle->SetBodyFixed(false); // Add geometry of the axle. axle->GetCollisionModel()->ClearModel(); chrono::utils::AddSphereGeometry(axle.get(), cmaterial, 0.5, ChVector<>(0, 0, 0)); axle->GetCollisionModel()->BuildModel(); sysMBS.AddBody(axle); // Connect the chassis to the containing bin (ground) through a translational joint and create a linear actuator. auto prismatic1 = chrono_types::make_shared<ChLinkLockPrismatic>(); prismatic1->Initialize(ground, chassis, ChCoordsys<>(chassis->GetPos(), Q_from_AngY(CH_C_PI_2))); prismatic1->SetName("prismatic_chassis_ground"); sysMBS.AddLink(prismatic1); double velocity = wheel_AngVel * wheel_radius * (1.0 - wheel_slip); auto actuator_fun = chrono_types::make_shared<ChFunction_Ramp>(0.0, velocity); actuator->Initialize(ground, chassis, false, ChCoordsys<>(chassis->GetPos(), QUNIT), ChCoordsys<>(chassis->GetPos() + ChVector<>(1, 0, 0), QUNIT)); actuator->SetName("actuator"); actuator->SetDistanceOffset(1); actuator->SetActuatorFunction(actuator_fun); sysMBS.AddLink(actuator); // Connect the axle to the chassis through a vertical translational joint. auto prismatic2 = chrono_types::make_shared<ChLinkLockPrismatic>(); prismatic2->Initialize(chassis, axle, ChCoordsys<>(chassis->GetPos(), QUNIT)); prismatic2->SetName("prismatic_axle_chassis"); sysMBS.AddLink(prismatic2); // Connect the wheel to the axle through a engine joint. motor->SetName("engine_wheel_axle"); motor->Initialize(wheel, axle, ChFrame<>(wheel->GetPos(), chrono::Q_from_AngAxis(-CH_C_PI / 2.0, ChVector<>(1, 0, 0)))); motor->SetAngleFunction(chrono_types::make_shared<ChFunction_Ramp>(0, wheel_AngVel)); sysMBS.AddLink(motor); } // ============================================================================= int main(int argc, char* argv[]) { // Create oputput directories if (!filesystem::create_directory(filesystem::path(out_dir))) { std::cerr << "Error creating directory " << out_dir << std::endl; return 1; } if (!filesystem::create_directory(filesystem::path(out_dir + "/particles"))) { std::cerr << "Error creating directory " << out_dir + "/particles" << std::endl; return 1; } if (!filesystem::create_directory(filesystem::path(out_dir + "/fsi"))) { std::cerr << "Error creating directory " << out_dir + "/fsi" << std::endl; return 1; } if (!filesystem::create_directory(filesystem::path(out_dir + "/vtk"))) { std::cerr << "Error creating directory " << out_dir + "/vtk" << std::endl; return 1; } // Create the MBS and FSI systems ChSystemSMC sysMBS; ChSystemFsi sysFSI(&sysMBS); ChVector<> gravity = ChVector<>(0, 0, -9.81); sysMBS.Set_G_acc(gravity); sysFSI.Set_G_acc(gravity); sysFSI.SetVerbose(verbose_fsi); // Use the default input file or you may enter your input parameters as a command line argument std::string inputJson = GetChronoDataFile("fsi/input_json/demo_FSI_SingleWheelTest.json"); if (argc == 3) { inputJson = std::string(argv[1]); wheel_slip = std::stod(argv[2]); } else if (argc != 1) { std::cout << "usage: ./FSI_ScrewTest <json_file> <wheel_slip>" << std::endl; std::cout << "or to use default input parameters ./FSI_ScrewTest " << std::endl; return 1; } sysFSI.ReadParametersFromFile(inputJson); // Set the initial particle spacing sysFSI.SetInitialSpacing(iniSpacing); // Set the SPH kernel length sysFSI.SetKernelLength(kernelLength); // Set the terrain density sysFSI.SetDensity(density); // Set the simulation stepsize sysFSI.SetStepSize(dT); // Set the terrain container size sysFSI.SetContainerDim(ChVector<>(bxDim, byDim, bzDim)); // Set SPH discretization type, consistent or inconsistent sysFSI.SetDiscreType(false, false); // Set wall boundary condition sysFSI.SetWallBC(BceVersion::ADAMI); // Set rigid body boundary condition sysFSI.SetRigidBodyBC(BceVersion::ADAMI); // Set cohsion of the granular material sysFSI.SetCohesionForce(1.0e2); // Setup the SPH method sysFSI.SetSPHMethod(FluidDynamics::WCSPH); // Set up the periodic boundary condition (if not, set relative larger values) ChVector<> cMin(-bxDim / 2 * 10, -byDim / 2 - 0.5 * iniSpacing, -bzDim * 10); ChVector<> cMax(bxDim / 2 * 10, byDim / 2 + 0.5 * iniSpacing, bzDim * 10); sysFSI.SetBoundaries(cMin, cMax); // Initialize the SPH particles ChVector<> boxCenter(0.0, 0.0, bzDim / 2); ChVector<> boxHalfDim(bxDim / 2, byDim / 2, bzDim / 2); sysFSI.AddBoxSPH(boxCenter, boxHalfDim); // Create Solid region and attach BCE SPH particles CreateSolidPhase(sysMBS, sysFSI); // Set simulation data output length sysFSI.SetOutputLength(0); // Construction of the FSI system must be finalized before running sysFSI.Initialize(); auto wheel = sysMBS.Get_bodylist()[1]; ChVector<> force = actuator->Get_react_force(); ChVector<> torque = motor->Get_react_torque(); ChVector<> w_pos = wheel->GetPos(); ChVector<> w_vel = wheel->GetPos_dt(); ChVector<> angvel = wheel->GetWvel_loc(); // Save wheel mesh ChTriangleMeshConnected wheel_mesh; wheel_mesh.LoadWavefrontMesh(GetChronoDataFile(wheel_obj), false, true); wheel_mesh.RepairDuplicateVertexes(1e-9); // Write the information into a txt file std::ofstream myFile; std::ofstream myDBP_Torque; if (output) { myFile.open(out_dir + "/results.txt", std::ios::trunc); myDBP_Torque.open(out_dir + "/DBP_Torque.txt", std::ios::trunc); } // Create a run-tme visualizer ChVisualizationFsi fsi_vis(&sysFSI); if (render) { fsi_vis.SetTitle("Chrono::FSI wheel demo"); fsi_vis.SetCameraPosition(ChVector<>(0, -5 * byDim, 5 * bzDim), ChVector<>(0, 0, 0)); fsi_vis.SetCameraMoveScale(0.05f); fsi_vis.EnableBoundaryMarkers(true); fsi_vis.Initialize(); } // Start the simulation unsigned int output_steps = (unsigned int)round(1 / (out_fps * dT)); unsigned int render_steps = (unsigned int)round(1 / (render_fps * dT)); double time = 0.0; int current_step = 0; ChTimer<> timer; timer.start(); while (time < total_time) { // Get the infomation of the wheel force = actuator->Get_react_force(); torque = motor->Get_react_torque(); w_pos = wheel->GetPos(); w_vel = wheel->GetPos_dt(); angvel = wheel->GetWvel_loc(); if (verbose) { std::cout << "time: " << time << std::endl; std::cout << " wheel position: " << w_pos << std::endl; std::cout << " wheel linear velocity: " << w_vel << std::endl; std::cout << " wheel angular velocity: " << angvel << std::endl; std::cout << " drawbar pull: " << force << std::endl; std::cout << " wheel torque: " << torque << std::endl; } if (output) { myFile << time << "\t" << w_pos.x() << "\t" << w_pos.y() << "\t" << w_pos.z() << "\t" << w_vel.x() << "\t" << w_vel.y() << "\t" << w_vel.z() << "\t" << angvel.x() << "\t" << angvel.y() << "\t" << angvel.z() << "\t" << force.x() << "\t" << force.y() << "\t" << force.z() << "\t" << torque.x() << "\t" << torque.y() << "\t" << torque.z() << "\n"; myDBP_Torque << time << "\t" << force.x() << "\t" << torque.z() << "\n"; } if (output && current_step % output_steps == 0) { std::cout << "-------- Output" << std::endl; sysFSI.PrintParticleToFile(out_dir + "/particles"); sysFSI.PrintFsiInfoToFile(out_dir + "/fsi", time); static int counter = 0; std::string filename = out_dir + "/vtk/wheel." + std::to_string(counter++) + ".vtk"; WritewheelVTK(filename, wheel_mesh, wheel->GetFrame_REF_to_abs()); } // Render SPH particles if (render && current_step % render_steps == 0) { if (!fsi_vis.Render()) break; } // Call the FSI solver sysFSI.DoStepDynamics_FSI(); time += dT; current_step++; } timer.stop(); std::cout << "\nSimulation time: " << timer() << " seconds\n" << std::endl; if (output) { myFile.close(); myDBP_Torque.close(); } return 0; }
