Hi Luning, Thanks for the reply. I have looked into links and how they are used by simulating different options for a short amount of time. I am a little confused as testing different options did not seem to make much sense.
For a prismatic joint is defined as: prismatic1->Initialize(ground, chassis, ChCoordsys<>(chassis->GetPos(),* QUNIT*)); The motion is allowed in the Z-dir. However, in the SingleWheel demo, the joint is defined as: prismatic1->Initialize(ground, chassis, ChCoordsys<>(chassis->GetPos(), *Q_from_AngY*(CH_C_PI_2))); which results in the motion being allowed in both the x and z directions (two directions). However, when changing the input to this prismatic1->Initialize(ground, chassis, ChCoordsys<>(chassis->GetPos(), *Q_from_AngX*(CH_C_PI_2))); The motion becomes restricted in both x and y and only free in z. *Why is that?* To try to navigate through this, I started by removing the actuator completely from my implementation as I only added it to be able to extract the forces on the screw (*Is there a way to extract the forces on the screw without having to use an actuator*?). Then, I have only two prismatic joints left ( one between the ground and the chassis and one between the chassis and axle). I started my tests by having the input for my prismatic joint as *Q_from_AngY*(CH_C_PI_2). This resulted in my screw moving in the z and x directions as expected. Then, I changed the input to *Q_from_AngX*(CH_C_PI_2), and also as expected, the screw was allowed to move in z-direction but not in x and y. As a result, I started thinking that maybe the second joint is what is preventing my screw from moving in certain directions. As a result, I tried the following combinations: 1 - *prismatic1*->Initialize(ground, chassis, ChCoordsys<>(chassis->GetPos(), Q_from_AngY(CH_C_PI_2))); * prismatic2*->Initialize(chassis, axle, ChCoordsys<>(chassis->GetPos(), QUNIT)); *( screw moved in both x and z)* *2- prismatic1*->Initialize(ground, chassis, ChCoordsys<>(chassis->GetPos(), QUNIT)); *prismatic2*->Initialize(chassis, axle, ChCoordsys<>(chassis->GetPos(), Q_from_AngY(CH_C_PI_2))); *(screw moved in x and z)* *3- **prismatic1*->Initialize(ground, chassis, ChCoordsys<>(chassis->GetPos(), QUNIT)); *prismatic2*->Initialize(chassis, axle, ChCoordsys<>(chassis->GetPos(), QUNIT)); *(screw moved only in z)* 4 - *prismatic1*->Initialize(ground, chassis, ChCoordsys<>(chassis->GetPos(), Q_from_AngXCH_C_PI_2))); * prismatic2*->Initialize(chassis, axle, ChCoordsys<>(chassis->GetPos(), QUNIT)); *( screw moved only in z)* *5- **prismatic1*->Initialize(ground, chassis, ChCoordsys<>(chassis->GetPos(), QUNIT)); *prismatic2*->Initialize(chassis, axle, ChCoordsys<>(chassis->GetPos(), Q_from_AngX(CH_C_PI_2))); *( screw moved only in z)* *6- **prismatic1*->Initialize(ground, chassis, ChCoordsys<>(chassis->GetPos(), Q_from_AngX(CH_C_PI_2) )); *prismatic2*->Initialize(chassis, axle, ChCoordsys<>(chassis->GetPos(), Q_from_AngX(CH_C_PI_2))); *( no movement)* I then thought that maybe the type of link I was using was not correct. I looked into different links (Project Chrono: Links <https://api.projectchrono.org/links.html>) and could not find one that seems appropriate. I have also looked at the demos you specified to try to better understand (crank demo) and it does not look like the prismatic link is used there. However, the ChLinkLockPointLine link seems like something that can be used in my case but I cannot tell what the difference is between this link and the prismatic link. As for the motor, the motor seems to work well so far as the screw rotates about its axis as it is supposed to. One final note, I have multiplied the mass obtained with ComputeMassProperties() by the screw density and got 0.9888. This is still quite different from what my CAD provides (1.71kg). Any insight about this? The inertia values are also off. Sorry about the very long message and thank you so much for all your help, On Monday, November 20, 2023 at 12:22:10 PM UTC-7 [email protected] wrote: > Hello Mohammad, > > In the single wheel test demo for FSI, the longitudinal direction is > global X, the lateral direction is global Y, the wheel's motion is > restricted in Y direction, while the X direction has a prescribed motor > function. > > When you connect bodies with links, be careful with how you define the > degree of freedom of the constraint. Specifically, for prismatic joint, it > allows translation along the Z axis. For instance in this line of your code, > > prismatic1->Initialize(ground, chassis, ChCoordsys<>(chassis->GetPos(), > Q_from_AngX(CH_C_PI_2))); > > The joint frame is aligned with Q_from_AngX(Ch_C_PI_2), which rotates 90 > degrees around global X, since the ground is fixed, the chassis can move > freely in the global Z direction, but not X or Y direction. However, the > actuator is applying a motion in X direction, but restricting motion in Y > and Z. That's why the screw is not moving in any direction. > > I suggest you play with some demos in the mbs folder, such as > demo_MBS_motor.cpp or demo_MBS_crank.cpp. The visualization can help you > better understand how joint frames are defined in Chrono. You can also have > a mbs model of the screw setup only (without the FSI terrain, but the > chassis-axel-screw setup). All you need is the code in the > CreateSolidPhase(). Again, with the visualization, it's easier to figure > out what was wrong. > > For your second question, note that ComputeMassProperties() only has the > information of the obj mesh, which is the geometry of the object, whatever > computed has to be multiplied by the density to get mass and inertia of the > object. > > Thank you, > Luning > On Saturday, November 18, 2023 at 3:31:56 PM UTC-6 [email protected] > wrote: > >> Apologies, I uploaded the wrong cpp file. Attached is my latest one. >> >> On Saturday, November 18, 2023 at 2:30:51 PM UTC-7 Mohammad Wasfi wrote: >> >>> Thank you, Luning, >>> >>> I have a few more questions that I hope you can answer for me. >>> >>> 1- In the single-wheel demo, the wheel's movement is restricted in the >>> x-direction, and the wheel is allowed to move in the y-direction. However, >>> for this simulation, I am trying to change the movement restriction for the >>> screw from the x to the y direction and allow the screw to move in the >>> x-direction. I have tried to change multiple things but nothing seems to >>> work as I intended. Right now, it seems that my screw is restricted from >>> moving in both the x and y directions. I am just trying to remove any >>> restrictions put on the screw in the X-direction as the screw tried to >>> propel itself and move. I was wondering if you could take a look at my >>> code and tell me if you see such a restriction set. >>> >>> 2- Similar to the single wheel demo, I use the function >>> " ComputeMassProperties" to compute the screw mass, cog, and inertia >>> tensor. For some reason, the results returned by that function are >>> different from what my CAD program provides for the same material. I was >>> wondering why that is. >>> >>> Example: >>> *Computed mass from code:* 0.00078467 ( I am not even sure what the >>> units of this should be. All parameters I used had units of kg/m/etc. I >>> assume that this is supposed to be kg. However, It is a very small value >>> compared to what CAD provides). >>> *computed mass from CAD*: 1.71 kg >>> >>> >>> Thank you so much, >>> >>> On Monday, November 13, 2023 at 10:09:42 AM UTC-7 [email protected] >>> wrote: >>> >>>> Hello Mohammad, >>>> >>>> From the error message, the job was terminated for exceeding time >>>> limit, meaning the job is still running, it just takes a very long time. >>>> From the output message, the FSI solver did not get initialized either. >>>> Your FSI parameters look good to me >>>> >>>> Without running your simulation, I think it's likely that something is >>>> wrong with CreateMeshPoints(), which is for generating BCE markers that >>>> describe the screw geometry. Since decreasing resolution will take longer >>>> time to generate the markers, the program could still be computing the >>>> markers. The implementation right now doesn't scale well with higher >>>> resolution. It might also got stuck in some infinite loop ... I will need >>>> the screw geometry to look into this more. >>>> >>>> You can also generate the BCE markers yourself, and read the points in >>>> the FSI program, then attach them to the screw body using AddPointsBCE(). >>>> That way you have control over how to place the grid points so it can >>>> capture the geometry. You can also speed up the process by parallelize >>>> your >>>> implementation. You only need to do it once. >>>> >>>> Thank you, >>>> Luning >>>> On Sunday, November 12, 2023 at 10:53:15 PM UTC-6 [email protected] >>>> wrote: >>>> >>>>> Thank you so much, Professor. I have tested the code with some other >>>>> initial spacing values and it seems that the code starts working again >>>>> when >>>>> I set my initial spacing value to 1.5*particle diameter (originally, the >>>>> initial spacing value was set to be 0.0033m while the particle diam is >>>>> 0.003m). It would still be helpful to look into why the code froze and no >>>>> errors were outputted whenever you have time for it. >>>>> >>>>> Thank you so much for all the help, >>>>> >>>>> On Sunday, November 12, 2023 at 1:00:57 PM UTC-7 Dan Negrut wrote: >>>>> >>>>>> Mohammad – we’ll try to take a look, things are backed up, hopefully >>>>>> we can help in some way… >>>>>> >>>>>> Dan >>>>>> >>>>>> --------------------------------------------- >>>>>> >>>>>> Bernard A. and Frances M. Weideman Professor >>>>>> >>>>>> NVIDIA CUDA Fellow >>>>>> >>>>>> Department of Mechanical Engineering >>>>>> >>>>>> Department of Computer Science >>>>>> >>>>>> University of Wisconsin - Madison >>>>>> >>>>>> 4150ME, 1513 University Avenue >>>>>> >>>>>> Madison, WI 53706-1572 >>>>>> >>>>>> 608 772 0914 <(608)%20772-0914> >>>>>> >>>>>> http://sbel.wisc.edu/ >>>>>> >>>>>> http://projectchrono.org/ >>>>>> >>>>>> --------------------------------------------- >>>>>> >>>>>> >>>>>> >>>>>> *From:* [email protected] <[email protected]> *On >>>>>> Behalf Of *Mohammad Wasfi >>>>>> *Sent:* Sunday, November 12, 2023 1:42 AM >>>>>> *To:* ProjectChrono <[email protected]> >>>>>> *Subject:* [chrono] FSI module >>>>>> >>>>>> >>>>>> >>>>>> Hello, >>>>>> >>>>>> >>>>>> >>>>>> I am trying to create a simulation based on the single-wheel demo. In >>>>>> my simulation, I am using a screw instead of a wheel. I started with >>>>>> changing the way the material properties are imported into the >>>>>> simulation. >>>>>> Everything worked fine as shown in the attached video. However, now, I >>>>>> am >>>>>> trying to reduce particle diameter ( as well as initial spacing). The >>>>>> simulation seems to freeze at the beginning and no error is printed. I >>>>>> was >>>>>> wondering if more insight could be provided into why the simulation >>>>>> freezes >>>>>> when reducing the particle diameter with no error printed. >>>>>> >>>>>> >>>>>> >>>>>> Thank you in advance for your help, >>>>>> >>>>>> -- >>>>>> 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/239386b1-d155-40a8-bf1c-6999db12484fn%40googlegroups.com >>>>>> >>>>>> <https://groups.google.com/d/msgid/projectchrono/239386b1-d155-40a8-bf1c-6999db12484fn%40googlegroups.com?utm_medium=email&utm_source=footer> >>>>>> . >>>>>> >>>>> -- 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/ae724ee8-5509-4bd8-bc11-c3e6517ad071n%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 // Particle diameter double particle_diameter = 0.005; // Set initial spacing as 1.5 times the particle diameter double iniSpacing = 2 * particle_diameter;//2*sphere rad //Initial separation of the fluid particles (suggested value is 10*particle_diam) double kernelLength = 1.0*iniSpacing;/// one or two times the initial particle spacing, it is used to calculate R; the scaled length parameter double density = 1260.0; //kg/m3 screw // Dimension of the terrain container double smalldis = 1.0e-9; double bxDim = 0.38 + smalldis; double byDim = 0.8 + smalldis; double bzDim = 0.14 + smalldis; // Size of the screw double screw_radius = 0.1;//m Diameter/2 + thread depth double screw_slip = 0.0; double screw_AngVel = 6.0; double total_mass = 1.71;;//// screw mass + caps . float I_XX = 0.01885; float I_YY = 0.00801; float I_ZZ = 0.01885; float I_p =(I_XX,I_YY,I_ZZ); std::string screw_obj = "vehicle/hmmwv/screw.obj"; // Initial Position of screw ChVector<> screw_IniPos(0.0, 0.0, screw_radius + bzDim + iniSpacing); ChVector<> screw_IniVel(0.0, 0.0, 0.0); // Simulation time and stepsize double total_time = 5.0; double dT = 1.0e-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_Single_Screw_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 Screw to Paraview VTK files //------------------------------------------------------------------ void WritescrewVTK(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(0.1e8); cmaterial->SetFriction(0.45f); cmaterial->SetRestitution(0.5f); cmaterial->SetAdhesion(0); // screw contact material auto cmaterial2 = chrono_types::make_shared<ChMaterialSurfaceSMC>(); cmaterial2->SetYoungModulus(4.1e9); cmaterial2->SetFriction(0.5f); cmaterial2->SetRestitution(0.275f); cmaterial2->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); //(Body, material, frame, size, thickness, faces, visualization) 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 screw -- always SECOND body in the system: auto trimesh = chrono_types::make_shared<geometry::ChTriangleMeshConnected>(); double scale_ratio = 1.0; trimesh->LoadWavefrontMesh(GetChronoDataFile(screw_obj), false, true); trimesh->Transform(ChVector<>(0, 0, 0), ChMatrix33<>(scale_ratio)); // scale to a different size trimesh->RepairDuplicateVertexes(1e-9); // 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); std::cout << "Computed mass: " << mmass *density << std::endl; std::cout << "Computed I: " << minertia *density<< std::endl; std::cout << "Computed Principle I: " << principal_I *density<< std::endl; std::cout << "Computed Principle I Rot: " << principal_inertia_rot *density<< std::endl; // Set the abs orientation, position and velocity auto screw = chrono_types::make_shared<ChBodyAuxRef>(); // Set the COG coordinates to barycenter, without displacing the REF reference. // Make the COG frame a principal frame. screw->SetFrame_COG_to_REF(ChFrame<>(mcog, principal_inertia_rot)); // Set inertia screw->SetMass(total_mass); screw->SetInertiaXX(ChVector<>(I_XX,I_YY,I_ZZ)); screw->SetPos(screw_IniPos); screw->SetPos_dt(screw_IniVel); screw->SetWvel_loc(ChVector<>(0.0, 0.0, 0.0)); // set an initial anular velocity (rad/s) ChQuaternion<> screw_Rot = Q_from_Euler123(ChVector<double>(0, 0, 0)); // Set the absolute position of the body: screw->SetFrame_REF_to_abs(ChFrame<>(ChVector<>(screw_IniPos), ChQuaternion<>(screw_Rot))); sysMBS.AddBody(screw); screw->SetBodyFixed(false); screw->GetCollisionModel()->ClearModel(); screw->GetCollisionModel()->AddTriangleMesh(cmaterial2, trimesh, false, false, VNULL, ChMatrix33<>(1), 0.005); screw->GetCollisionModel()->BuildModel(); screw->SetCollide(false); // Add this body to the FSI system std::vector<ChVector<>> BCE_screw; sysFSI.CreateMeshPoints(*trimesh, iniSpacing, BCE_screw); sysFSI.AddPointsBCE(screw, BCE_screw, ChFrame<>(), true); sysFSI.AddFsiBody(screw); // Create the chassis -- always THIRD body in the system auto chassis = chrono_types::make_shared<ChBody>(); chassis->SetMass( 1.0 ); chassis->SetPos(screw->GetPos()); chassis->SetCollide(false); chassis->SetBodyFixed(false); // Add geometry of the chassis. chassis->GetCollisionModel()->ClearModel(); chrono::utils::AddBoxGeometry(chassis.get(), cmaterial2, 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(1.0); axle->SetPos(screw->GetPos()); axle->SetCollide(false); axle->SetBodyFixed(false); // Add geometry of the axle. axle->GetCollisionModel()->ClearModel(); chrono::utils::AddSphereGeometry(axle.get(), cmaterial2, 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)));// Do Not Understand prismatic1->SetName("prismatic_chassis_ground"); sysMBS.AddLink(prismatic1); // 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));// Do Not Understand prismatic2->SetName("prismatic_axle_chassis"); sysMBS.AddLink(prismatic2); // Connect the screw to the axle through a engine joint. motor->SetName("engine_screw"); motor->Initialize(screw, axle, ChFrame<>(screw->GetPos(), chrono::Q_from_AngAxis(-CH_C_PI / 2.0, ChVector<>(1, 0, 0))));//here motor->SetAngleFunction(chrono_types::make_shared<ChFunction_Ramp>(0, screw_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 a physic system to handle multibody dynamics ChSystemSMC sysMBS; sysMBS.SetContactForceModel(ChSystemSMC::ContactForceModel::Hooke); sysMBS.SetTimestepperType(ChTimestepper::Type::EULER_IMPLICIT); // Create an FSI system to handle fluid dynamics ChSystemFsi sysFSI(&sysMBS); ChVector<> gravity = ChVector<>(0, 0, -9.81); sysMBS.Set_G_acc(gravity); sysFSI.Set_G_acc(gravity); sysFSI.SetVerbose(verbose_fsi); // Set the initial particle spacing sysFSI.SetInitialSpacing(iniSpacing); // Set the SPH kernel length sysFSI.SetKernelLength(kernelLength); ChSystemFsi::ElasticMaterialProperties mat_props; mat_props.Young_modulus = 0.10e8; mat_props.Poisson_ratio = 0.45; mat_props.viscosity_alpha = 0.5;//Artifical viscosity coefficient mat_props.viscosity_beta = 0.0;//Artifical viscosity coefficient mat_props.mu_I0 = 0.04;//Reference Inertia number mat_props.mu_fric_s = 0.75;//friction mu_s mat_props.mu_fric_2 = 0.75;//mu_2 constant in mu=mu(I) mat_props.average_diam = particle_diameter; mat_props.friction_angle = CH_C_PI / 10; // default//Frictional angle of granular material mat_props.dilation_angle = CH_C_PI / 10; // default//Dilate angle of granular material mat_props.cohesion_coeff = 0; //?? // default mat_props.kernel_threshold = 0.8;//Threshold of the integration of the kernel function sysFSI.SetElasticSPH(mat_props); sysFSI.SetInitialSpacing(iniSpacing); // Set the SPH kernel length sysFSI.SetKernelLength(kernelLength); // one or two times the initial particle spacing // Set the terrain density sysFSI.SetDensity(2.5e3);//https://link.springer.com/chapter/10.1007/978-3-030-22818-7_3 // 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);///< Lower limit point ChVector<> cMax(bxDim / 2 * 10, byDim / 2 + 0.5 * iniSpacing, bzDim * 10);///< Upper limit point 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 screw = sysMBS.Get_bodylist()[1]; // ChVector<> force = actuator->Get_react_force(); // ChVector<> torque = motor->Get_react_torque(); ChVector<> w_pos = screw->GetPos(); ChVector<> w_vel = screw->GetPos_dt(); ChVector<> angvel = screw->GetWvel_loc(); // Save screw mesh ChTriangleMeshConnected screw_mesh; screw_mesh.LoadWavefrontMesh(GetChronoDataFile(screw_obj), false, true); screw_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 single screw 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 screw //force = actuator->Get_react_force(); //torque = motor->Get_react_torque(); w_pos = screw->GetPos(); w_vel = screw->GetPos_dt(); angvel = screw->GetWvel_loc(); if (verbose) { std::cout << "time: " << time << std::endl; std::cout << " screw position: " << w_pos << std::endl; std::cout << " screw linear velocity: " << w_vel << std::endl; std::cout << " screw angular velocity: " << angvel << std::endl; // std::cout << " drawbar pull: " << force << std::endl; // std::cout << " screw 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/screw." + std::to_string(counter++) + ".vtk"; WritescrewVTK(filename, screw_mesh, screw->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; }
