Hi Ruochun,
Thank you so much for your help. This actually helped a lot. I have some
questions that I hope you could answer. I am trying to control the mesh
motion using *ApplyMeshMotion *as you mentioned. However, I feel like I do
not fully understand everything related to that so I wanted to ask some
questions and hopefully, this will help others as well. I have looked at
the documentation and they are really brief so I hope you could answer some
of these questions here. I will use examples from my code(attached here). I
might be using some stuff incorrectly so I think it would be a good way to
see if what I am doing is correct or not.
1- gpu_sys.AddMesh("../Screw_CM.obj", *ChVector<float>(0,0,-90)*,
*ChMatrix33<float>(1)*, screw_mass);
- *ChVector<float>(0,0,-90): *does this represent the
initial position where the mesh is initiated, and does it affect other
parts of the simulation?
- *ChMatrix33<float>(1): *what does this represent? I
looked through the documentation and it says that this is the rotscale.
Could you give a little more explanation as to
what that is?
2- ChVector<> ball_initial_pos(0, 0, 0);
std::shared_ptr<ChBody> screw_body(sys_screw.NewBody());
screw_body->SetMass(screw_mass);
screw_body->SetInertiaXX(ChVector<>(screw_inertia_x,
screw_inertia_y, screw_inertia_z));
screw_body->*SetPos(ball_initial_pos)*;
screw_body->SetBodyFixed(false);
sys_screw.AddBody(screw_body);
- When creating a system using the imported mesh, what does
screw_body->*SetPos(ball_initial_pos) *do? Clearly, it is a way to set the
position of the mesh. However, since the
*ChVector<float>(0,0,-90)
*was used before when adding the mesh, is there a need to use *SetPos?*
*3- * ChVector<float> *mesh_pos(0, 0, 0); *
float rev_per_sec = 10.f;
float ang_vel_Z = rev_per_sec * 2 * (float)CH_C_PI;
ChVector<> *mesh_lin_vel*(0);
ChQuaternion<> *mesh_rot *= Q_from_AngZ(t * ang_vel_Z);
ChVector<> *mesh_ang_vel*(0, 0, ang_vel_Z);
gpu_sys.ApplyMeshMotion(0, mesh_pos, mesh_rot, mesh_lin_vel,
mesh_ang_vel);
- *mesh_pos(0, 0, 0):* to me, it seems that we are setting
the position of the mesh a third time. could you explain how this affects
the simulation?
- *mesh_lin_vel*(0): does this represent the initial
velocity of the mesh? or is it applying a constant velocity throughout the
entire simulation?
- *mesh_rot : *Could you explain what this exactly is
and how changing it affects the simulation?
* - **mesh_ang_vel*(0, 0, ang_vel_Z): is this the
initial rotational speed or it is a constant rotational speed that the mesh
will have throughout the entire
simulation?
Thank you so much in advance for taking the time to help us with
everything,
On Saturday, September 3, 2022 at 12:09:24 PM UTC-7 Ruochun Zhang wrote:
> Hi Mohammad,
>
> I think that *demo_GPU_mixer.cpp* demo does exactly what you needed.
> While you are looking into this demo, I just want to mention one thing to
> avoid confusion: in terms of controlling the mesh motion in C::GPU, the
> *Transform
> *call before *AddMesh *should only be used to "align" the mesh with its
> CoM before the simulation starts, so that when the mesh is loaded it is
> expressed in its centroid and principal frame; it should not be used to set
> the initial location of the mesh. Instead, like in the demo, use
> *ApplyMeshMotion
> *to control the mesh motion throughout the simulation.
>
> Hope it helps!
> Ruochun
>
> On Saturday, September 3, 2022 at 10:58:19 AM UTC-5 [email protected]
> wrote:
>
>> Hello there,
>>
>> I am trying to perform a simulation using the gpu module where I want to
>> import a cylindrical screw (picture attached) into a bed on granular media
>> and then have the screw spin on top of the bed and then obtain the
>> forces/torques on the screw.
>>
>> I am not sure about what would be the best way to import the screw in a
>> way that it's set on the bed at the start of the simulation. Also, I am not
>> sure how I can impose a type of motion that would force the screw to spin
>> on top of that bed since I am importing it as an .obj file. Finally, due to
>> the shape of the screw, I assume that when the screw moves it will tend to
>> turn to the left/right. Is there a way to restrict the motion of the screw
>> to be just straight while spinning? I was wondering if there are any demos
>> regarding this or if someone has done something similar before. Any help
>> will be appreciated.
>>
>> Thank you so much 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/a6a3153c-6b84-470f-9be2-caaa3ad008d6n%40googlegroups.com.
#include <cmath>
#include <iomanip>
#include <iostream>
#include <string>
#include "chrono/core/ChGlobal.h"
#include "chrono/physics/ChBody.h"
#include "chrono/physics/ChSystemSMC.h"
#include "chrono/utils/ChUtilsSamplers.h"
#include "chrono_gpu/physics/ChSystemGpu.h"
#include "chrono_gpu/utils/ChGpuJsonParser.h"
#include "chrono_gpu/ChGpuData.h"
// added for motor:
#include "chrono/motion_functions/ChFunction.h"
#include "chrono/physics/ChForce.h"
#include "chrono/timestepper/ChTimestepper.h"
#include "chrono/utils/ChUtilsCreators.h"
#include "chrono/physics/ChLinkMotorLinearForce.h"
#include "chrono/physics/ChLinkMotorLinearPosition.h"
#include "chrono/physics/ChLinkMotorLinearSpeed.h"
#include "chrono/physics/ChLinkMotorRotationAngle.h"
#include "chrono/physics/ChLinkMotorRotationSpeed.h"
#include "chrono/physics/ChLinkMotorRotationTorque.h"
#include "chrono_gpu/ChGpuData.h"
#include "chrono_thirdparty/filesystem/path.h"
using namespace chrono;
using namespace chrono::gpu;
//////////////////////////////////////////////////////////////////////////////////////////////////////////
// Simulation parameters (same as those in the json file)
double step_size = 5e-7; // time step s
double sphere_radius = 0.2; // cm sphere radius (6 mm glass beads same as
experiments )
double sphere_density = 0.92; // g/cm^3 density of sphere particle
double time_settle =0.75;// time needed for bed to settle
double time_end = 0.1;// time needed for the drop phase
// Define Big domain size
double box_X = 32;
double box_Y = 32;
double box_Z = 200;
//int run_mode_id = 0;//settling phase
int run_mode_id = 1;//drop phase
double youngs_modulus1 = 1e10;
// Gravity
double grav_X = 0.0f;
double grav_Y = 0.0f;
double grav_Z = -981.0f;
float screw_mass = 13522; //g value from Ansys
///////////////////////////////////////////////////////////////////////////////////////////////////////
double getVoidRatio( ChSystemGpuMesh& gpu_sys, double sphere_radius){
double vol_sphere = 4.0f/3.0f * chrono::CH_C_PI * std::pow(sphere_radius,
3) * gpu_sys.GetNumParticles () ;
double vol_box = (chrono::CH_C_PI*box_X*box_Y/4)* (gpu_sys.GetMaxParticleZ
() - gpu_sys.GetMinParticleZ ());
double eta = vol_sphere/vol_box;
double voidRatio = (1.0f - eta)/eta;
return voidRatio;
}
///////////////////////////////////////////////////////////////////////////////////////////////////////
void write_csv(std::string filename, std::vector<std::pair<std::string,
std::vector<float>>> dataset){
// Make a CSV file with one or more columns of integer values
// Each column of data is represented by the pair <column name, column data>
// as std::pair<std::string, std::vector<int>>
// The dataset is represented as a vector of these columns
// Note that all columns should be the same size
// Create an output filestream object
std::ofstream myFile(filename);
// Send column names to the stream
for(int j = 0; j < dataset.size(); ++j)
{
myFile << dataset.at(j).first;
if(j != dataset.size() - 1) myFile << ","; // No comma at end of line
}
myFile << "\n";
// Send data to the stream
for(int i = 0; i < dataset.at(0).second.size(); ++i)
{
for(int j = 0; j < dataset.size(); ++j)
{
myFile << dataset.at(j).second.at(i);
if(j != dataset.size() - 1) myFile << ","; // No comma at end of
line
}
myFile << "\n";
}
// Close the file
myFile.close();
}
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// Set up GPU system -bed Particles
void SetupGPUSystem(ChSystemGpuMesh& gpu_sys) {
double cor_p = 0.36; // sphere restitution
double cor_w = 0.36; // mesh & wall restitution (value from Thoesen et
al., 2019)
double cor_m = 0.8; // mesh restitution
double youngs_modulus_p = 1e10; //(cms^2)
double youngs_modulus_m = 1.93e12;
double youngs_modulus_w = 1e10;
double poisson_ratio_m = 0.265;
double poisson_ratio_p = 0.32;
double poisson_ratio_w = 0.32;
double mu_s2s = 0.5; // 0.5 static friction coefficient sphere to sphere
double mu_s2w = 0.5; //0.5 static friction coefficient sphere to wall &
mesh (value from Thoesen et al., 2019)
double mu_s2m = 0.04;
float roll_s2s = 1; //1
float roll_s2w = 1;
float roll_s2m = 0.1;
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
gpu_sys.UseMaterialBasedModel(true); // use material based model (Young's
modulus, restitution, Poisson's ratio)
gpu_sys.SetYoungModulus_SPH(youngs_modulus_p);
gpu_sys.SetYoungModulus_WALL(youngs_modulus_w);
gpu_sys.SetYoungModulus_MESH(youngs_modulus_m);
gpu_sys.SetRestitution_SPH(cor_p);
gpu_sys.SetRestitution_WALL(cor_w);
gpu_sys.SetRestitution_MESH(cor_m);
gpu_sys.SetPoissonRatio_SPH(poisson_ratio_p);
gpu_sys.SetPoissonRatio_WALL(poisson_ratio_w);
gpu_sys.SetPoissonRatio_MESH(poisson_ratio_m);
gpu_sys.SetFrictionMode(chrono::gpu::CHGPU_FRICTION_MODE::MULTI_STEP);
gpu_sys.SetStaticFrictionCoeff_SPH2SPH(mu_s2s);
gpu_sys.SetStaticFrictionCoeff_SPH2WALL(mu_s2w);
gpu_sys.SetStaticFrictionCoeff_SPH2MESH(mu_s2m);
gpu_sys.SetRollingMode(CHGPU_ROLLING_MODE::SCHWARTZ);
gpu_sys.SetRollingCoeff_SPH2SPH(roll_s2s);
gpu_sys.SetRollingCoeff_SPH2WALL(roll_s2w);
gpu_sys.SetRollingCoeff_SPH2MESH(roll_s2m);
gpu_sys.SetParticleOutputMode(CHGPU_OUTPUT_MODE::CSV);
gpu_sys.SetTimeIntegrator(CHGPU_TIME_INTEGRATOR::CENTERED_DIFFERENCE);
gpu_sys.SetFixedStepSize(step_size);
gpu_sys.SetBDFixed(true);
/////////////////////////////////////////////////////////////////////////////////////////////
// Rain particles
double fill_top;
double fill_bottom = -5.0f;
double box_xy = 32; // 32 cm by 32 cm box
double box_r = box_xy / 2;
double spacing = 2.001 * sphere_radius;
std::vector<ChVector<float>> body_points;
utils::PDSampler<float> sampler(spacing);
fill_top = box_Z / 2 - spacing;
ChVector<> hdims(box_r - sphere_radius, box_r - sphere_radius, 0);
int counter = 0;
for (double z = fill_bottom + spacing; z < fill_top; z += spacing) {
ChVector<> center(0, 0, z);
auto points = sampler.SampleBox(center, hdims);
body_points.insert(body_points.end(), points.begin(), points.end());
counter = counter + points.size();
if (counter > 5000) {
break;
}
}
std::cout << "Created " << body_points.size() << "spheres" << std::endl;
gpu_sys.SetParticles(body_points);
gpu_sys.SetGravitationalAcceleration(ChVector<float>(grav_X, grav_Y,
grav_Z));
}
///////////////////////////////////////////////////////////////////////////////////////////////////////////////
void SetupGPUParams(ChSystemGpuMesh& gpu_sys) {
double cor_p = 0.36; // sphere restitution
double cor_w = 0.36; // mesh & wall restitution (value from Thoesen et
al., 2019)
double cor_m = 0.8; // mesh restitution
double youngs_modulus_p = 1e10; // 40Mpa = 4e7Pa = 4e8 g/(cms^2)
double youngs_modulus_m = 1.93e12;
double youngs_modulus_w = 1e10;
double poisson_ratio_m = 0.265;
double poisson_ratio_p = 0.32;
double poisson_ratio_w = 0.32;
double mu_s2s = 0.5; // 0.5 static friction coefficient sphere to sphere
double mu_s2w = 0.5; //0.5 static friction coefficient sphere to wall &
mesh (value from Thoesen et al., 2019)
double mu_s2m = 0.04;
float roll_s2s = 1; //1
float roll_s2w = 1;
float roll_s2m = 0.1;
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
gpu_sys.UseMaterialBasedModel(true); // use material based model (Young's
modulus, restitution, Poisson's ratio)
gpu_sys.SetYoungModulus_SPH(youngs_modulus_p);
gpu_sys.SetYoungModulus_WALL(youngs_modulus_w);
gpu_sys.SetYoungModulus_MESH(youngs_modulus_m);
gpu_sys.SetRestitution_SPH(cor_p);
gpu_sys.SetRestitution_WALL(cor_w);
gpu_sys.SetRestitution_MESH(cor_m);
gpu_sys.SetPoissonRatio_SPH(poisson_ratio_p);
gpu_sys.SetPoissonRatio_WALL(poisson_ratio_w);
gpu_sys.SetPoissonRatio_MESH(poisson_ratio_m);
gpu_sys.SetGravitationalAcceleration(ChVector<float>(grav_X, grav_Y,
grav_Z));
gpu_sys.SetFrictionMode(chrono::gpu::CHGPU_FRICTION_MODE::MULTI_STEP);
gpu_sys.SetStaticFrictionCoeff_SPH2SPH(mu_s2s);
gpu_sys.SetStaticFrictionCoeff_SPH2WALL(mu_s2w);
gpu_sys.SetStaticFrictionCoeff_SPH2MESH(mu_s2m);
gpu_sys.SetRollingMode(CHGPU_ROLLING_MODE::SCHWARTZ);
gpu_sys.SetRollingCoeff_SPH2SPH(roll_s2s);
gpu_sys.SetRollingCoeff_SPH2WALL(roll_s2w);
gpu_sys.SetRollingCoeff_SPH2MESH(roll_s2m);
//gpu_sys.SetCohesionEnergyDensity(CED);
gpu_sys.SetParticleOutputMode(CHGPU_OUTPUT_MODE::CSV);
gpu_sys.SetTimeIntegrator(CHGPU_TIME_INTEGRATOR::CENTERED_DIFFERENCE);
gpu_sys.SetFixedStepSize(step_size);
gpu_sys.SetBDFixed(true);
}
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
int main(int argc, char* argv[]) {
// Output directory
std::string out_dir = GetChronoOutputPath() + "GPU/";
filesystem::create_directory(filesystem::path(out_dir));
out_dir = out_dir + "ballcosim";
filesystem::create_directory(filesystem::path(out_dir));
std::string checkpoint_file = out_dir + "checkpoint.dat";
// run_mode_id == 1, this is simulation phase
if (run_mode_id == 1)
{
// Load Checkpoint file
ChSystemGpuMesh gpu_sys(checkpoint_file);
// Material based model is not included in the checkpoint file
// Setup GPU system simulation parameters
SetupGPUParams(gpu_sys);
//float inital_location = gpu_sys.GetMinParticleZ+10.0f;
// Setup meshes
gpu_sys.AddMesh("../Screw_CM.obj", ChVector<float>(0,0,-90),
ChMatrix33<float>(1), screw_mass);
gpu_sys.EnableMeshCollision(true);
gpu_sys.Initialize();
std::cout << gpu_sys.GetNumMeshes() << " meshes" << std::endl;
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
ChSystemSMC sys_screw;
sys_screw.SetContactForceModel(ChSystemSMC::ContactForceModel::Hooke);
sys_screw.SetTimestepperType(ChTimestepper::Type::EULER_IMPLICIT);
sys_screw.Set_G_acc(ChVector<>(grav_X, grav_Y, grav_Z));
float screw_inertia_x = 2.0462e6;//g·cm² 0.1004
float screw_inertia_y = 3.4689e6;//g·cm² 0.1004
float screw_inertia_z = 2.0462e6;//g·cm² 0.1
ChVector<> ball_initial_pos(0, 0, 0);
//screw body
std::shared_ptr<ChBody> screw_body(sys_screw.NewBody());
screw_body->SetMass(screw_mass);
screw_body->SetInertiaXX(ChVector<>(screw_inertia_x, screw_inertia_y,
screw_inertia_z));
screw_body->SetPos(ball_initial_pos);
screw_body->SetBodyFixed(false);
sys_screw.AddBody(screw_body);
// Setup detailed output options
float out_fps = 200;
std::cout << "Output at " << out_fps << " FPS" << std::endl;
unsigned int out_steps = (unsigned int)(1.0 / (out_fps * step_size));
unsigned int total_frames = (unsigned int)(time_end * out_fps);
int currframe = 0;
unsigned int curr_step = 0;
std::vector<float> screw_fx(total_frames);
std::vector<float> screw_fy(total_frames);
std::vector<float> screw_fz(total_frames);
clock_t start = std::clock();
for (double t = 0; t < (double)time_end;
t += step_size, curr_step++) {
//// Apply motion to the meshes
// initial positions and velocity (not used)
ChVector<float> mesh_pos(0, 0, 0);
float rev_per_sec = 150.f;
float ang_vel_Z = rev_per_sec * 2 * (float)CH_C_PI;
ChVector<> mesh_lin_vel(0);
ChQuaternion<> mesh_rot = Q_from_AngZ(t * ang_vel_Z);
ChVector<> mesh_ang_vel(0, 0, ang_vel_Z);
gpu_sys.ApplyMeshMotion(0, mesh_pos, mesh_rot, mesh_lin_vel,
mesh_ang_vel);
ChVector<float> screw_pos(0, 0, 0);
screw_pos.z() = screw_body->GetPos().z();
screw_pos.y() = screw_body->GetPos().y();
screw_pos.x() = screw_body->GetPos().x();
gpu_sys.ApplyMeshMotion(0, screw_body->GetPos(), screw_body->GetRot(),
screw_body->GetPos_dt(), screw_body->GetWvel_par());
//// calculate forces on the tip
ChVector<> screw_force;
ChVector<> screw_torque;
gpu_sys.CollectMeshContactForces(0, screw_force, screw_torque);
screw_body->Empty_forces_accumulators();
screw_body->Accumulate_force(screw_force, screw_body->GetPos(),
false);
screw_body->Accumulate_torque(screw_torque, false);
if (curr_step % out_steps == 0) {
std::cout << "Output frame " << currframe + 1 << " of " <<
total_frames << std::endl;
char filename[100];
char mesh_filename[100];
const ChVector<>& pos= screw_body->GetPos();
std::cout << "Position of Screw in z direction" << pos <<
std::endl;
screw_fx[currframe] = screw_force.x();
screw_fy[currframe] = screw_force.y();
screw_fz[currframe] = screw_force.z();
sprintf(filename, "%s/step%06d.csv", out_dir.c_str(), currframe);
sprintf(mesh_filename, "%s/step%06d_mesh", out_dir.c_str(),
currframe++);
gpu_sys.WriteParticleFile(std::string(filename));
gpu_sys.WriteMeshes(std::string(mesh_filename));
}
gpu_sys.AdvanceSimulation(step_size);
sys_screw.DoStepDynamics(step_size);
}
std::vector<std::pair<std::string, std::vector<float>>> vals = {{"Fx",
screw_fx}, {"Fy", screw_fy}, {"Fz", screw_fz}};
write_csv("screw_forces.csv",vals);
clock_t end = std::clock();
double total_time = ((double)(end - start)) / CLOCKS_PER_SEC;
std::cout << "Time: " << total_time << " seconds" << std::endl;
return 0;
}
// run_mode_id == 0, this is the sample preparation phase
ChSystemGpuMesh gpu_sys(sphere_radius, sphere_density,
ChVector<float>(box_X, box_Y, box_Z));
SetupGPUSystem(gpu_sys);
gpu_sys.EnableMeshCollision(false);
gpu_sys.Initialize();
float curr_timee = 0;
unsigned int currframe = 0;
float out_fps = 100;
unsigned int out_steps = (unsigned int)(1 / (out_fps * step_size));
unsigned int total_frames = (unsigned int)(time_settle * out_fps);
unsigned int curr_step = 0;
for (double t = 0; t < (double)time_settle; t += step_size, curr_step++) {
if (curr_step % out_steps == 0) {
std::cout << "Output frame " << currframe + 1 << " of " <<
total_frames << std::endl;
char filename[100];
sprintf(filename, "%s/step%06d.csv", out_dir.c_str(), currframe++);
gpu_sys.WriteParticleFile(std::string(filename));
}
gpu_sys.AdvanceSimulation(step_size);
curr_timee += step_size;
}
double vr1 = getVoidRatio( gpu_sys, sphere_radius);
double pr1 = vr1/(1.0f+vr1);
std::cout << "Void_Ratio " << vr1 << std::endl;
std::cout << "Porosity " << pr1<< std::endl;
std::cout << "MAX_z " << gpu_sys.GetMaxParticleZ () << std::endl;
std::cout << "Min_z " << gpu_sys.GetMinParticleZ () << std::endl;
std::cout << "youngs_modulus " << youngs_modulus1 << std::endl;
//////////////////////////////////////////////////////////////////////////////////////////////////
// This is settling phase, so output a checkpoint file
gpu_sys.WriteCheckpointFile(checkpoint_file);
return 0;
}
