Hello Ruochun,
This is a DEME-related question.
I am trying to compress a cylindrical bed of particles similar to what is
done in the cone penetration demo. However, for some reason, the plane that
is used to compress the bed does not seem to be moving or to be interacting
with the bed material. I have tried many things but nothing seems to be
working. I have attached my code for your reference.
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/58b3a741-7955-4cdb-a4d7-18e1be50ae39n%40googlegroups.com.
//Cone Penetraion test: Mohammad
#include <core/ApiVersion.h>
#include <core/utils/ThreadManager.h>
#include <DEM/API.h>
#include <DEM/HostSideHelpers.hpp>
#include <DEM/utils/Samplers.hpp>
#include <cstdio>
#include <chrono>
#include <filesystem>
using namespace deme;
using namespace std::filesystem;
///////////////////////////////////////////////////////////////////////////////////////////////////////
int main() {
DEMSolver DEMSim;
DEMSim.SetVerbosity(INFO);
DEMSim.SetOutputFormat(OUTPUT_FORMAT::CSV);
DEMSim.SetOutputContent(OUTPUT_CONTENT::ABSV);
DEMSim.SetMeshOutputFormat(MESH_FORMAT::VTK);
DEMSim.SetContactOutputContent(OWNER | FORCE | POINT);
// Material Properties for cone and particles:
// E, nu, CoR, mu, Crr, CED...
auto mat_type_cone = DEMSim.LoadMaterial({{"E", 2e5}, {"nu", 0.27}, {"CoR",
0.88}, {"mu", 0.5}, {"Crr", 1.00}});
auto mat_type_particle = DEMSim.LoadMaterial({{"E", 1e9}, {"nu", 0.32},
{"CoR", 0.36}, {"mu", 0.5}, {"Crr", 1.00}});
// Cone & cylindrical bin Parameters:
double time_settle =0.4;
double time_compress = 0.55;
double time_decompress = 0.65;
//float cone_speed = 0.1;//m/s
float step_size = 5e-7;//s
double world_size = 4;
double soil_bin_diameter = 0.1; //m
//double cone_diameter = 0.01;//m
//float shaft_mass = 47.436 ;
//float tip_mass = 1.8137 ;
//double starting_height = -0.02;
float particle_radius = 0.002; //m
float sphere_density = 920;//kg/m3
// 1- Particles Initiation:
auto template_sphere = DEMSim.LoadSphereType(particle_radius *
particle_radius * particle_radius * sphere_density * 4 / 3 * 3.14159,
particle_radius,
mat_type_particle);
std::vector<std::shared_ptr<DEMClumpTemplate>> clump_types;
clump_types.push_back(template_sphere);
// Generate initial clumps for piling
int num_template = 1;
// Generate initial clumps for piling
float spacing = 2.05*particle_radius;
float fill_width = soil_bin_diameter / 2 - 4*particle_radius;
float fill_height = 0.4;
float fill_bottom = -0.05 + fill_width + spacing;
PDSampler sampler(spacing);
// Use a PDSampler-based clump generation process
std::vector<std::shared_ptr<DEMClumpTemplate>> input_pile_template_type;
std::vector<float3> input_pile_xyz;
float layer_z = 0;
while (layer_z < fill_height) {
float3 sample_center = make_float3(0, 0, fill_bottom + layer_z +
spacing / 2);
auto layer_xyz = sampler.SampleCylinderZ(sample_center, fill_width, 0);
unsigned int num_clumps = layer_xyz.size();
// Select from available clump types
for (unsigned int i = 0; i < num_clumps; i++) {
input_pile_template_type.push_back(clump_types.at(i %
num_template));
}
input_pile_xyz.insert(input_pile_xyz.end(), layer_xyz.begin(),
layer_xyz.end());
layer_z += spacing;
}
// Calling AddClumps a to add clumps to the system
auto the_pile = DEMSim.AddClumps(input_pile_template_type, input_pile_xyz);
size_t n_particles = input_pile_xyz.size();
std::cout << "Added clump type " << n_particles << std::endl;
// Create a inspector to find out stuff
auto max_z_finder = DEMSim.CreateInspector("clump_max_z");
auto min_z_finder = DEMSim.CreateInspector("clump_min_z");
float max_z;
float min_z;
// 2- Domain Initiation:
DEMSim.InstructBoxDomainDimension(world_size, world_size, world_size);
// No need to add simulation `world' boundaries, b/c we'll add a
cylinderical container manually
DEMSim.InstructBoxDomainBoundingBC("none", mat_type_particle);
// Now add a cylinderical boundary along with a bottom plane
auto walls = DEMSim.AddExternalObject();
walls->AddCylinder(make_float3(0), make_float3(0, 0, 1), soil_bin_diameter
/ 2., mat_type_particle, 0);
walls->AddPlane(make_float3(0, 0, -0.052), make_float3(0, 0, 1),
mat_type_particle);
// 3- Compressor:
// Now add a plane to compress the sample
auto compressor = DEMSim.AddExternalObject();
compressor->AddPlane(make_float3(0, 0, 0.5), make_float3(0, 0, -1),
mat_type_particle);
compressor->SetFamily(10);
DEMSim.DisableContactBetweenFamilies(0, 10);
DEMSim.SetFamilyFixed(10);
auto compressor_tracker = DEMSim.Track(compressor);
DEMSim.SetCoordSysOrigin("center");
// 3- Simulation Parameters & initilization:
DEMSim.SetInitTimeStep(step_size);
DEMSim.SetGravitationalAcceleration(make_float3(0, 0, -9.81));
DEMSim.SetCDUpdateFreq(6);
DEMSim.SetMaxVelocity(25.);
DEMSim.SetExpandSafetyParam(2.);
DEMSim.SetInitBinSize(2* particle_radius);
DEMSim.SetIntegrator(TIME_INTEGRATOR::CENTERED_DIFFERENCE);
DEMSim.Initialize();
DEMSim.UpdateSimParams();
// 4-Runing Simulation:
path out_dir = current_path();
out_dir += "/DemoOutput_Pile";
create_directory(out_dir);
float curr_time_of_sim = 0; // keeps track of simulation time
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_decompress * 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) {
char filename[100];
sprintf(filename, "%s/step%06d.csv", out_dir.c_str(), currframe++);
DEMSim.WriteSphereFile(std::string(filename));
std::cout << "Output frame " << currframe + 1 << " of " << total_frames
<< std::endl;
}
DEMSim.DoDynamics(step_size);
curr_time_of_sim += step_size;
}
auto total_volume_finder = DEMSim.CreateInspector("clump_volume");
max_z = max_z_finder->GetValue();
min_z = min_z_finder->GetValue();
float porosity;
float total_volume_of_spheres = 4 / 3 * 3.14159
*particle_radius*particle_radius*particle_radius*n_particles;
float total_volume = 3.14159 * (soil_bin_diameter * soil_bin_diameter / 4)
* (max_z - min_z);
double eta = total_volume_of_spheres / total_volume;
double voidRatio = (1.0f - eta)/eta;
porosity = voidRatio/(1.0f+voidRatio);
std::cout << "Max Z is: " << max_z << std::endl;
std::cout << "Min Z is: " << min_z << std::endl;
std::cout << "Porosity Before Compression: " << porosity << std::endl;
std::cout << "Void Ratio Before Compression: " << voidRatio << std::endl;
std::cout << "curr_time_of_sim1: " << curr_time_of_sim << std::endl;
DEMSim.EnableContactBetweenFamilies(0, 10);
//Compresse
std::cout << "here1: " << std::endl;
double compressor_vel = 1;
float terrain_max_z = max_z_finder->GetValue();
compressor_tracker->SetPos(make_float3(0, 0, terrain_max_z));
curr_step = 0;
for (double t = curr_time_of_sim; t < (double)time_compress; t +=
step_size, curr_step++) {
float matter_volume = 4 / 3 * 3.14159
*particle_radius*particle_radius*particle_radius*n_particles;
float total_volume2 = 3.14159 * (soil_bin_diameter * soil_bin_diameter
/ 4) * (terrain_max_z - min_z);
double eta = matter_volume / total_volume2;
double voidRatio = (1.0f - eta)/eta;
porosity = voidRatio/(1.0f+voidRatio);
float3 Position_of_compressor = compressor_tracker->Pos();
if (curr_step % out_steps == 0) {
char filename[100];
sprintf(filename, "%s/step%06d.csv", out_dir.c_str(), currframe++);
DEMSim.WriteSphereFile(std::string(filename));
std::cout << "Compression porosity: " << porosity << std::endl;
std::cout << "Position_of_compressor " << Position_of_compressor.z
<< std::endl;
std::cout << "Void Ratio After Compression: " << voidRatio <<
std::endl;
std::cout << "Output frame " << currframe + 1 << " of " <<
total_frames << std::endl;
}
terrain_max_z -= compressor_vel * step_size;
compressor_tracker->SetPos(make_float3(0, 0, terrain_max_z));
DEMSim.DoDynamics(step_size);
curr_step++;
curr_time_of_sim += step_size;
}
std::cout << "curr_time_of_sim2: " << curr_time_of_sim << std::endl;
double decompressor_vel = 0.1;
terrain_max_z = max_z_finder->GetValue();
compressor_tracker->SetPos(make_float3(0, 0, terrain_max_z));
curr_step = 0;
for (double t = curr_time_of_sim; t < (double)time_decompress; t +=
step_size, curr_step++) {
float matter_volume = 4 / 3 * 3.14159
*particle_radius*particle_radius*particle_radius*n_particles;
float total_volume2 = 3.14159 * (soil_bin_diameter * soil_bin_diameter
/ 4) * (terrain_max_z - min_z);
double eta = matter_volume / total_volume2;
double voidRatio = (1.0f - eta)/eta;
porosity = voidRatio/(1.0f+voidRatio);
float3 Position_of_compressor = compressor_tracker->Pos();
if (curr_step % out_steps == 0) {
char filename[100];
sprintf(filename, "%s/step%06d.csv", out_dir.c_str(), currframe++);
DEMSim.WriteSphereFile(std::string(filename));
std::cout << "DeCompression porosity: " << porosity << std::endl;
std::cout << "Position_of_compressor " << Position_of_compressor.z
<< std::endl;
std::cout << "Void Ratio After DeCompression: " << voidRatio <<
std::endl;
std::cout << "Output frame " << currframe + 1 << " of " <<
total_frames << std::endl;
}
terrain_max_z += decompressor_vel * step_size;
compressor_tracker->SetPos(make_float3(0, 0, terrain_max_z));
DEMSim.DoDynamics(step_size);
curr_step++;
curr_time_of_sim += step_size;
}
std::cout << "curr_time_of_sim3: " << curr_time_of_sim << std::endl;
//Remove compressor
DEMSim.DoDynamicsThenSync(step_size);
DEMSim.DisableContactBetweenFamilies(0, 3);
DEMSim.DoDynamicsThenSync(step_size);
terrain_max_z = max_z_finder->GetValue();
std::cout << "Max Z is: " << terrain_max_z << std::endl;
//DEMSim.ShowTimingStats();
//DEMSim.ClearTimingStats();
std::cout << "Exiting..." << std::endl;
return 0;
}
