This is the file that I am currently using. I have my calculations for the
void ratio at the very end if you need to look at it.
Thank you so much
On Sunday, May 22, 2022 at 3:44:57 PM UTC-7 Mohammad Wasfi wrote:
> Hello there,
>
> I am using the GPU module to initiate a bed of particles and let it
> settle. I have been able to do that successfully. Now, I am trying to
> control the void volume of the bed using the void ratio value. To do so, I
> need to determine the maximum and minimum x and y positions of the
> particles to find the most accurate cross-sectional area of the bed. I was
> wondering of a way to do this similarly to using the getMaxParticleZ. If
> that's not possible, is there a Chrono function under the GPU module that
> can directly find the void ratio?
>
> Finally, do you have any suggestions for controlling the void volume of a
> bed of particles?
>
>
> Thank you so much,
>
--
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/f8c19276-bbae-464b-a9d8-5cc94177cc00n%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"
#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
double sphere_radius = 0.1; // sphere radius (6 mm glass beads same as
experiments )
double sphere_density = 0.92; // density of sphere particle
double time_settle = 0.75;
double time_end = 0.5;
// Define Big domain size
double box_X = 16;
double box_Y = 16;
double box_Z = 200;
int run_mode_id = 0;
// Gravity
double grav_X = 0.0f;
double grav_Y = 0.0f;
double grav_Z = -981.0f;
double ISE = 100000;
float cone_density = 8; //g/cm^3
float cone_mass = (1. / 3.) * (float)CH_C_PI * 5.0 * 5.0 * 5.0 * cone_density;
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.8; // sphere restitution
double cor_w = 0.8; // mesh & wall restitution (value from Thoesen et al.,
2019)
double cor_m = 0.8; // mesh restitution
double youngs_modulus_p = 9.7e10; // 40Mpa = 4e7Pa = 4e8 g/(cms^2)
double youngs_modulus_m = 1.9e12;
double youngs_modulus_w = 9.7e10;
double poisson_ratio_m = 0.265;
double poisson_ratio_p = 0.35;
double poisson_ratio_w = 0.35;
double mu_s2s = 0.5; // static friction coefficient sphere to sphere
double mu_s2w = 0.5; // static friction coefficient sphere to wall & mesh
(value from Thoesen et al., 2019)
double mu_s2m = 0.04;
float roll_s2s = 0.001;
float roll_s2w = 0.01;
float roll_s2m = 0.01;
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.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
double fill_bottom_below = -95;
double fill_top_below = -65;
chrono::utils::PDSampler<float> sampler(2.4f * sphere_radius);
// leave a 1cm margin at edges of sampling
ChVector<> hdims(box_X / 2 - 2.0, box_Y / 2 - 2.0, 0);
ChVector<> center(0, 0, fill_bottom_below + 2.0 * sphere_radius);
std::vector<ChVector<float>> body_points;
// Shift up for bottom of box
center.z() += 3 * sphere_radius;
// fill the lower half
while (center.z() < fill_top_below) {
std::cout << "Create layer at " << center.z() << std::endl;
auto points = sampler.SampleBox(center, hdims);
body_points.insert(body_points.end(), points.begin(), points.end());
center.z() += 2.05 * sphere_radius;
}
//Send sinusoidal wave through particle bed in x direction after 0.5
seconds?
std::function<double3(float)> pos_func_wave = [¶ms](float t ) {
double3 pos = {0, 0, 0};
double t0 = 0.5;
double freq = CH_C_PI / 4;
if (t > t0) {
pos.x = 0.1 * params.box_X * std::sin((t - t0) * freq);
}
return pos;
}
gpu_sys.setBDWallsMotionFunction(pos_func_wave);
std::cout << "Created " << body_points.size() << "spheres" << std::endl;
gpu_sys.SetParticles(body_points);
}
void SetupGPUParams(ChSystemGpuMesh& gpu_sys) {
double cor_p = 0.8; // sphere restitution
double cor_w = 0.8; // mesh & wall restitution (value from Thoesen et al.,
2019)
double cor_m = 0.8; // mesh restitution
double youngs_modulus_p = 9.7e10; // 40Mpa = 4e7Pa = 4e8 g/(cms^2)
double youngs_modulus_m = 1.9e12;
double youngs_modulus_w = 9.7e10;
double poisson_ratio_m = 0.265;
double poisson_ratio_p = 0.35;
double poisson_ratio_w = 0.35;
double mu_s2s = 0.5; // static friction coefficient sphere to sphere
double mu_s2w = 0.5; // static friction coefficient sphere to wall & mesh
(value from Thoesen et al., 2019)
double mu_s2m = 0.04;
float roll_s2s = 0.001;
float roll_s2w = 0.01;
float roll_s2m = 0.01;
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.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("../Part1.obj", ChVector<float>(0, 0, -78.5173),
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;
sys_cone.SetContactForceModel(ChSystemSMC::ContactForceModel::Hooke);
sys_cone.SetTimestepperType(ChTimestepper::Type::EULER_IMPLICIT);
sys_cone.Set_G_acc(ChVector<>(0, 0, -980));
float cone_inertia_x = (3. / 10.) * cone_mass * 5.0 * 5.0;
float cone_inertia_y = cone_mass * ((3. / 20.) * 5.0 * 5.0 + (3. / 80.)
* 5.0 * 5.0);
float cone_inertia_z = cone_inertia_y;
ChVector<> ball_initial_pos(0, 0, 0);
std::shared_ptr<ChBody> cone_body(sys_cone.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.AddBody(cone_body);
// Setup detailed output options
float out_fps = 350;
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> 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(0, cone_body->GetPos(), cone_body->GetRot(),
cone_body->GetPos_dt(), cone_body->GetWvel_par());
//// calculate forces on the cone
ChVector<> cone_force;
ChVector<> cone_torque;
gpu_sys.CollectMeshContactForces(0, cone_force, cone_torque);
cone_body->Empty_forces_accumulators();
cone_body->Accumulate_force(cone_force, cone_body->GetPos(), false);
cone_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];
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);
}
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);
}
double height = gpu_sys.GetMaxParticleZ();
std::cout << height << std::endl;
// void Fraction calcs Mohammad:
double Volume_Of_SettledParticles = (box_X*box_Y)* (gpu_sys.GetMaxParticleZ
() - gpu_sys.GetMinParticleZ ());
double Total_Volume_Of_Particles = (1.334 * 3.14 *
pow(sphere_radius,3))*gpu_sys.GetNumParticles () ;
double voids_volume = Volume_Of_SettledParticles -
Total_Volume_Of_Particles;
double vr2 = voids_volume/Total_Volume_Of_Particles;
std::cout << "Void_Ratio " << vr2 << std::endl;
// This is settling phase, so output a checkpoint file
gpu_sys.WriteCheckpointFile(checkpoint_file);
return 0;
}
