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/75172f6c-5e60-4fd4-98f9-1f10ed5dce9an%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 = 1.0; double total_mass = 1.71;;//// screw mass + caps . float I_XX = 0.01885; float I_YY = 0.00801; float I_ZZ = 0.01885; 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 = 10.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<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); //trimesh -> SetMass(screw_mass); // 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);// Why are values differnt than what cad has ChMatrix33<> principal_inertia_rot; ChVector<> principal_I; ChInertiaUtils::PrincipalInertia(minertia, principal_I, principal_inertia_rot); std::cout << "Computed COG: " << mcog << std::endl; //mcog = ChVector<>(0.0, 0.0, 0.0); std::cout << "Computed mass: " << mmass << std::endl; std::cout << "Computed I: " << minertia << std::endl; std::cout << "Computed Principle I: " << principal_I << std::endl; // Set the abs orientation, position and velocity auto screw = chrono_types::make_shared<ChBodyAuxRef>(); ChQuaternion<> screw_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. screw->SetFrame_COG_to_REF(ChFrame<>(mcog, principal_inertia_rot)); // Set inertia screw->SetMass(total_mass); screw->SetInertiaXX(mdensity * principal_I); screw->SetPos_dt(screw_IniVel); screw->SetWvel_loc(ChVector<>(0.0, 0.0, 0.0)); // set an initial anular velocity (rad/s) // 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_AngX(CH_C_PI_2))); prismatic1->SetName("prismatic_chassis_ground"); sysMBS.AddLink(prismatic1); double velocity = screw_AngVel * screw_radius * (1.0 - screw_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<>(0, 1, 0), QUNIT));//here 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 screw to the axle through a engine joint. motor->SetName("engine_screw_axle"); 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; // 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; }
