Hi,
Thank you so much for your help. I really appreciate it. I will definitely
rebuild Chrono again using the latest version and will let you know if any
changes occur. For the cone, I have implemented a motor that ensures a
constant velocity of the cone. My understanding is that the motor applies
as much force as needed to maintain a constant velocity of the cone which
is the reason why the cone does not stop in the bed. I have generated
another video using ParaView
(https://drive.google.com/file/d/1q-hgTcJY62kLblqr37sFoHX-4goOwbzR/view?usp=sharing)
and it shows that the cone pushes the particles away as it moves. Also, non
of the particles get inside the cone. I have included the .cpp file needed
to run this simulation to avoid any confusion.
Forgive me I had to upload the video as a link due to size limitations for
the post.
Thank you so much for your help and please let me know if there is anything
else you want me to provide to you,
Mohammad,
--
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/6d42e6d0-d9f3-4612-b5ea-ad853ac2b082n%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;
double time_end = 1.75;
// Define Big domain size
double box_X = 10;
double box_Y = 10;
double box_Z = 37;
//int run_mode_id = 0;
int run_mode_id = 1;
// Gravity
double grav_X = 0.0f;
double grav_Y = 0.0f;
double grav_Z = -981.0f;
//double ISE = 0;
//double CED = 0;
double ISE = 100000;
double CED = 200000000;
float cone_density = 8; //g/cm^3
float cone_mass = 47.436 ;
float tip_mass = 1.8137 ;
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
void SetupGPUSystem(ChSystemGpuMesh& gpu_sys) {
// Set GPU system parameters
// Currently, parameters are given by Luning
// After calibration, they may be modified
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.88; // mesh restitution
double youngs_modulus_p = 1e8; // 40Mpa = 4e7Pa = 4e8 g/(cms^2)
double youngs_modulus_m = 2e6;
double youngs_modulus_w = 1e8;
double poisson_ratio_m = 0.27;
double poisson_ratio_p = 0.32;
double poisson_ratio_w = 0.32;
// double mu_s2s = 0.1; // 0.5 static friction coefficient sphere to sphere
// double mu_s2w = 0.1; //0.5 static friction coefficient sphere to wall &
mesh (value from Thoesen et al., 2019)
// double mu_s2m = 0.1;
// float roll_s2s = 0.1; //1
// float roll_s2w = 0.1;
// float roll_s2m = 0.1;
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.5;
float roll_s2s = 1; //1
float roll_s2w = 1;
float roll_s2m = 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.SetInterfacialSurfaceEnergy(ISE);
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);
//gpu_sys.CreateBCPlane(ChVector<>(0, 0, -6.5), ChVector<>(0, 0, 1), false);
// Rain particles
// create cylinder containter
ChVector<float> cyl_center(0.0f, 0.0f, 0.0f);
float cyl_rad = std::min(box_X, box_Y) / 2.0f;
gpu_sys.CreateBCCylinderZ(cyl_center, cyl_rad, false, true);
// generate a cloud of particles
std::vector<chrono::ChVector<float>> body_points;
utils::PDSampler<float> sampler(2.001 * sphere_radius);
ChVector<float> sampler_center(0.0f, 0.0f, 0.0f);
body_points = sampler.SampleCylinderZ(sampler_center, cyl_rad - 4 *
sphere_radius, box_Z / 2 - 4 * sphere_radius);
int numSpheres = body_points.size();
std::cout << "Numbers of particles created1: " << numSpheres << std::endl;
gpu_sys.SetParticles(body_points);
std::cout << "Created " << body_points.size() << "spheres" << std::endl;
gpu_sys.SetParticles(body_points);
}
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.88; // mesh restitution
double youngs_modulus_p = 1e8; // 40Mpa = 4e7Pa = 4e8 g/(cms^2)
double youngs_modulus_m = 2e6;
double youngs_modulus_w = 1e8;
double poisson_ratio_m = 0.27;
double poisson_ratio_p = 0.32;
double poisson_ratio_w = 0.32;
// double mu_s2s = 0.1; // 0.5 static friction coefficient sphere to sphere
// double mu_s2w = 0.1; //0.5 static friction coefficient sphere to wall &
mesh (value from Thoesen et al., 2019)
// double mu_s2m = 0.1;
// float roll_s2s = 0.1; //1
// float roll_s2w = 0.1;
// float roll_s2m = 0.1;
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.5;
float roll_s2s = 1; //1
float roll_s2w = 1;
float roll_s2m = 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.SetInterfacialSurfaceEnergy(ISE);
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);
// Setup meshes
// Load the cone
gpu_sys.AddMesh("../cone_tip.obj", ChVector<float>(0, 0, -12 ),
ChMatrix33<float>(1), tip_mass);
gpu_sys.EnableMeshCollision(true);
gpu_sys.Initialize();
std::cout << gpu_sys.GetNumMeshes() << " meshes" << std::endl;
// Create rigid ball_body simulation
ChSystemSMC sys_cone;
sys_cone.SetContactForceModel(ChSystemSMC::ContactForceModel::Hooke);
sys_cone.SetTimestepperType(ChTimestepper::Type::EULER_IMPLICIT);
sys_cone.Set_G_acc(ChVector<>(0, 0, -980));
//cone Tip --> change position
gpu_sys.AddMesh("../cone_shaft.obj", ChVector<float>(0, 0, -12),
ChMatrix33<float>(1), cone_mass);
gpu_sys.EnableMeshCollision(true);
gpu_sys.Initialize();
std::cout << gpu_sys.GetNumMeshes() << " meshes" << std::endl;
// Create rigid ball_body simulation
ChSystemSMC sys_cone_2;
sys_cone_2.SetContactForceModel(ChSystemSMC::ContactForceModel::Hooke);
sys_cone_2.SetTimestepperType(ChTimestepper::Type::EULER_IMPLICIT);
sys_cone_2.Set_G_acc(ChVector<>(0, 0, -980));
// set inertia values
float cone_inertia_x = 673.76 ; // (3. / 10.) * cone_mass * 5.0 * 5.0;
float cone_inertia_y = 673.76; // cone_mass * ((3. / 20.) * 5.0 * 5.0 +
(3. / 80.) * 5.0 * 5.0);
float cone_inertia_z = 3.4199; // cone_inertia_y;
float tip_inertia_x = 0.11666; // (3. / 10.) * cone_mass * 5.0 * 5.0;
float tip_inertia_y = 0.11666; // cone_mass * ((3. / 20.) * 5.0 * 5.0 +
(3. / 80.) * 5.0 * 5.0);
float tip_inertia_z = 0.13513; // cone_inertia_y;
// Creation of slider(slave) and guide (master) bodies
// guide body only used when using a motor
ChVector<> ball_initial_pos(0, 0, 0); // for slider
ChVector<> ball_initial_pos2(0, 0, 0); // for master
// Creation of rotation matrix. it is used to set the orientaion of the
master with respect to the slider
ChQuaternion<> mesh_rot = Q_from_AngAxis(-90 * CH_C_DEG_TO_RAD, VECT_Y);
ChVector<> mesh_trans(ball_initial_pos+ball_initial_pos2);
ChFrame<> Xb(mesh_trans,mesh_rot);
//shaft:
std::shared_ptr<ChBody> cone_body(sys_cone_2.NewBody());
cone_body->SetMass(cone_mass);
cone_body->SetInertiaXX(ChVector<>(cone_inertia_x, cone_inertia_y,
cone_inertia_z));
cone_body->SetPos(ball_initial_pos);
cone_body->SetBodyFixed(false);
sys_cone_2.AddBody(cone_body);
//Tip
std::shared_ptr<ChBody> tip_body(sys_cone.NewBody());
tip_body->SetMass(tip_mass);
tip_body->SetInertiaXX(ChVector<>(tip_inertia_x, tip_inertia_y,
tip_inertia_z));
tip_body->SetPos(ball_initial_pos);
tip_body->SetBodyFixed(false);
sys_cone.AddBody(tip_body);
//guide1
std::shared_ptr<ChBody> guide_body(sys_cone_2.NewBody());
guide_body->SetMass(cone_mass);
guide_body->SetInertiaXX(ChVector<>(cone_inertia_x, cone_inertia_y,
cone_inertia_z));
guide_body->SetPos(ball_initial_pos2);
guide_body->SetBodyFixed(true); //must be fixed
sys_cone_2.AddBody(guide_body);
//guide2
std::shared_ptr<ChBody> guide_body_2(sys_cone.NewBody());
guide_body_2->SetMass(tip_mass);
guide_body_2->SetInertiaXX(ChVector<>(tip_inertia_x, tip_inertia_y,
tip_inertia_z));
guide_body_2->SetPos(ball_initial_pos2);
guide_body_2->SetBodyFixed(true); //must be fixed
sys_cone.AddBody(guide_body_2);
// 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);
float lin_freq = 0.5f;
float lin_period = 1.f / lin_freq;
float lin_amp = 1.0; //??
float lin_vel = 10;
// Create and initialize the linear motor
auto motor_lin_0 = chrono_types::make_shared<ChLinkMotorLinearSpeed>();
motor_lin_0->Initialize(
tip_body, // body A (slave)
guide_body_2, // body B (master)
Xb); // motor frame, in abs. coords
// Add the motor to the system
sys_cone.AddLink(motor_lin_0);
// Create a ChFunction to be used for the ChLinkMotorLinearSpeed
// lin_rep_seq is a step function for the linear velocity control
auto lin_part0 = chrono_types::make_shared<ChFunction_Const>();
lin_part0->Set_yconst(-lin_vel);
auto lin_seq_0 = chrono_types::make_shared<ChFunction_Sequence>();
lin_seq_0->InsertFunct(lin_part0, time_end, 1, false);
auto lin_rep_seq_0 = chrono_types::make_shared<ChFunction_Repeat>();
lin_rep_seq_0->Set_fa(lin_seq_0);
lin_rep_seq_0->Set_window_length(lin_period);
lin_rep_seq_0->Set_window_start(0.0);
lin_rep_seq_0->Set_window_phase(lin_period);
// Let the linearmotor use this motion function:
motor_lin_0->SetSpeedFunction(lin_rep_seq_0);
/////////////////////////////////////////////////////
// Create and initialize the linear motor
auto motor_lin_1 = chrono_types::make_shared<ChLinkMotorLinearSpeed>();
motor_lin_1->Initialize(
cone_body, // body A (slave)
guide_body, // body B (master)
Xb); // motor frame, in abs. coords
// Add the motor to the system
sys_cone_2.AddLink(motor_lin_1);
// Create a ChFunction to be used for the ChLinkMotorLinearSpeed
// lin_rep_seq is a step function for the linear velocity control
// auto lin_part1 = chrono_types::make_shared<ChFunction_Const>();
// lin_part1->Set_yconst(-lin_vel);
// auto lin_seq_1 = chrono_types::make_shared<ChFunction_Sequence>();
// lin_seq_1->InsertFunct(lin_part1, params.time_end, 1, false);
// auto lin_rep_seq_1 = chrono_types::make_shared<ChFunction_Repeat>();
// lin_rep_seq_1->Set_fa(lin_seq_1);
// lin_rep_seq_1->Set_window_length(lin_period);
// lin_rep_seq_1->Set_window_start(0.0);
// lin_rep_seq_1->Set_window_phase(lin_period);
// Let the linearmotor use this motion function:
motor_lin_1->SetSpeedFunction(lin_rep_seq_0);
int currframe = 0;
unsigned int curr_step = 0;
std::vector<float> cone_fx(total_frames);
std::vector<float> cone_fy(total_frames);
std::vector<float> cone_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
gpu_sys.ApplyMeshMotion(1, cone_body->GetPos(), cone_body->GetRot(),
cone_body->GetPos_dt(), cone_body->GetWvel_par());
gpu_sys.ApplyMeshMotion(0, tip_body->GetPos(), tip_body->GetRot(),
tip_body->GetPos_dt(), tip_body->GetWvel_par());
//// calculate forces on the cone
ChVector<> cone_force;
ChVector<> cone_torque;
gpu_sys.CollectMeshContactForces(0, cone_force, cone_torque);
tip_body->Empty_forces_accumulators();
tip_body->Accumulate_force(cone_force, tip_body->GetPos(), false);
tip_body->Accumulate_torque(cone_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 = cone_body->GetPos();
std::cout << "Position of cone in z direction" << pos <<
std::endl;
const ChVector<>& pos_2 = tip_body->GetPos();
std::cout << "Position of tip in z direction" << pos_2 <<
std::endl;
cone_fx[currframe] = cone_force.x();
cone_fy[currframe] = cone_force.y();
cone_fz[currframe] = cone_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_cone.DoStepDynamics(step_size);
sys_cone_2.DoStepDynamics(step_size);
}
std::vector<std::pair<std::string, std::vector<float>>> vals = {{"Fx",
cone_fx}, {"Fy", cone_fy}, {"Fz", cone_fz}};
write_csv("cone_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));
// One thing we can do is to move the Big Box Domain by (0, 0, Z/2) using
SetBDCenter, so the
// coordinate range we are now working with is (-X/2,-Y/2,0) to
(X/2,Y/2,Z), instead of (-X/2,-Y/2,-Z/2) to (X/2, Y/2, Z/2).
// gpu_sys.SetBDCenter(ChVector<float>(0, 0, box_Z / 2));
SetupGPUSystem(gpu_sys);
gpu_sys.EnableMeshCollision(true);
gpu_sys.Initialize();
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);
}
// void Fraction calcs Mohammad:
double Volume_Of_SettledParticles = (chrono::CH_C_PI*box_X*box_Y/4)*
(gpu_sys.GetMaxParticleZ () - gpu_sys.GetMinParticleZ ());
double Total_Volume_Of_Particles = (1.334 * chrono::CH_C_PI *
std::pow(sphere_radius,3))*gpu_sys.GetNumParticles () ;
double voids_volume = Volume_Of_SettledParticles -
Total_Volume_Of_Particles;
double porosity = Total_Volume_Of_Particles/Volume_Of_SettledParticles;
double vr1 = (1.0f - porosity)/porosity;
std::cout << "Void_Ratio1 " << vr1 << std::endl;
double pr = vr1/(1.0f+vr1);
double vr2 = voids_volume/Total_Volume_Of_Particles;
std::cout << "Void_Ratio2 " << vr2 << std::endl;
std::cout << "MAX_z " << gpu_sys.GetMaxParticleZ () << std::endl;
std::cout << "Min_z " << gpu_sys.GetMinParticleZ () << std::endl;
std::cout << "NOP " << gpu_sys.GetNumParticles () << std::endl;
std::cout << "porosity " << pr << std::endl;
// This is settling phase, so output a checkpoint file
gpu_sys.WriteCheckpointFile(checkpoint_file);
return 0;
}
