Hi,
I am trying to use the Chrono gpu to achieve a similar effect as the run a
cone penetration test but have been running into some difficulties. I hope
that somebody here would be able to help me out.
When setting the young's modulus for the particles and boundaries to 1e9, I
get the following error:
- No available contact pair slots for body xxxxx and body xxxxx
However, when I run the simulation with values larger than 1e9, the
simulation runs smoothly.
It seems strange to me that there seems to exist a critical value of
cohesion that causes models to behave in unexpected ways. I am wondering if
anybody here has had success with using low values for young's modulus in
Chrono GPU. I am attaching my file for your reference.
--
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/24445804-3aad-4821-b6f4-52f02b9d6bbfn%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;
}
