From 0d8472188ca2f6eb8704a3aeb3f996701de98dae Mon Sep 17 00:00:00 2001 From: DanNegrut Date: Thu, 3 Sep 2026 18:54:23 -0500 Subject: [PATCH] Add DEMdemo_PlateSinkage: a two-stage plate pressure-sinkage test A rigid circular plate is pressed at constant speed into a granular bed and the pressure it carries is recorded against sinkage, which is the classical bevameter experiment terramechanics uses to characterize soil. The demo's structure is the part meant for reuse. Stage 1 prepares the bed cohesionless and saves it with WriteClumpFile; bed preparation is typically around half the cost of a run, and a saved bed serves every later test that wants the same packing. Stage 2 rebuilds a fresh solver, reloads the bed by clump-type name (positions and orientations both, since restoring clump centers without rotations silently changes the packing's interlocking), switches to the stock ForceModelWithCohesion, re-equilibrates, and presses. Preparing cohesionless and testing cohesive is deliberate: cohesion acting during preparation props the packing open and yields a looser bed than asked for. Uses only stock assets: the 3_clump template with the same mass and MOI numbers as DEMdemo_ConePenetration, the cyl_r1_h2 mesh scaled into a disc, and the shipped cohesion force model. Output is a pressure-versus-sinkage CSV plus a console table, with the sinkage datum taken at first contact (mean plate pressure reaching 1 kPa). Co-Authored-By: Claude Opus 5 (1M context) --- src/demo/CMakeLists.txt | 1 + src/demo/DEMdemo_PlateSinkage.cpp | 274 ++++++++++++++++++++++++++++++ 2 files changed, 275 insertions(+) create mode 100644 src/demo/DEMdemo_PlateSinkage.cpp diff --git a/src/demo/CMakeLists.txt b/src/demo/CMakeLists.txt index bffb73e5..80e70e18 100644 --- a/src/demo/CMakeLists.txt +++ b/src/demo/CMakeLists.txt @@ -35,6 +35,7 @@ SET(DEMOS DEMdemo_WheelDPSimplified DEMdemo_WheelSlopeSlip DEMdemo_Indentation + DEMdemo_PlateSinkage DEMdemo_SolarSystem DEMdemo_Electrostatic DEMdemo_FlexibleMesh diff --git a/src/demo/DEMdemo_PlateSinkage.cpp b/src/demo/DEMdemo_PlateSinkage.cpp new file mode 100644 index 00000000..2563d754 --- /dev/null +++ b/src/demo/DEMdemo_PlateSinkage.cpp @@ -0,0 +1,274 @@ +// Copyright (c) 2021, SBEL GPU Development Team +// Copyright (c) 2021, University of Wisconsin - Madison +// +// SPDX-License-Identifier: BSD-3-Clause + +// ============================================================================= +// A plate pressure-sinkage test, the classical bevameter experiment: a rigid +// circular plate is pressed at constant speed into a granular bed and the +// pressure it carries is recorded against how deep it has sunk. The output is +// the pressure-sinkage curve, which is how terramechanics characterizes a soil. +// +// The demo runs in two stages, and the split is the part worth copying: +// +// Stage 1 prepares the bed cohesionless (pour, settle) and saves it with +// WriteClumpFile. Bed preparation is typically around half the cost of a +// run, and a saved bed can be reused by every later test that wants the +// same packing. +// +// Stage 2 rebuilds a fresh solver, reloads the bed (positions AND +// orientations: these are clumped grains, and restoring centers without +// rotations silently changes the interlocking), switches to a force model +// WITH cohesion, lets the bed re-equilibrate, and only then presses the +// plate. Preparing cohesionless and loading cohesive is deliberate: if +// cohesion acts while the bed is being prepared, it props the packing open +// and you get a looser, weaker bed than the density you asked for. +// +// The cohesion model is the stock ForceModelWithCohesion.cu that ships with +// DEME: a constant attractive acceleration between contacting grains, set as +// the pairwise material property "Cohesion". +// ============================================================================= + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +using namespace deme; +using namespace std::filesystem; + +const double math_PI = 3.14159; + +int main() { + // --------------------------------------------------------------------- + // Everything a user might want to change, in one place. + // --------------------------------------------------------------------- + const float plate_diam = 0.35; // plate diameter, m + const float plate_thickness = 0.02; + const float plate_speed = 0.05; // press speed, m/s; slow it down for quasi-static work + const float target_sinkage = 0.06; // stop once the plate is this deep, m + + const double bin_diameter = 1.05; // soil bin, 3 plate diameters across + const float bottom = -0.5; // bin floor, in world coordinates + const float fill_height = 0.55; // loose pour column; settles to roughly 0.28 m of bed + + const float terrain_density = 2.6e3; // grain material density, kg/m^3 + const float cohesion = 200.; // stage-2 grain-grain cohesion, m/s^2 (an acceleration: + // the model applies force = Cohesion * effective mass) + const float step_size = 1e-5; + const double clump_scale = 0.01; // scales the stock 3-sphere template to ~22 mm grains, + // about 16 grain diameters under the plate + + path out_dir = current_path(); + out_dir /= "DemoOutput_PlateSinkage"; + create_directory(out_dir); + const std::string bed_file = (out_dir / "bed.csv").string(); + const std::string curve_file = (out_dir / "pressure_sinkage.csv").string(); + const std::string template_name = "bed_grain"; + + float bed_surface = 0.; // measured at the end of stage 1, reused in stage 2 + + // --------------------------------------------------------------------- + // Stage 1: prepare the bed cohesionless and save it. + // --------------------------------------------------------------------- + { + DEMSolver DEMSim; + DEMSim.SetVerbosity(INFO); + DEMSim.SetOutputFormat(OUTPUT_FORMAT::CSV); + // QUAT matters: the saved file must carry orientations or the reload + // in stage 2 cannot restore the packing's interlocking. + DEMSim.SetOutputContent(OUTPUT_CONTENT::XYZ | OUTPUT_CONTENT::QUAT | OUTPUT_CONTENT::VEL); + + auto mat_type_terrain = DEMSim.LoadMaterial({{"E", 1e9}, {"nu", 0.3}, {"CoR", 0.4}, {"mu", 0.5}, {"Crr", 0.01}}); + + const double world_size = 2; + DEMSim.InstructBoxDomainDimension(world_size, world_size, world_size); + DEMSim.InstructBoxDomainBoundingBC("top_open", mat_type_terrain); + // A cylindrical bin, so the circular plate sees the same wall distance everywhere. + auto walls = DEMSim.AddExternalObject(); + walls->AddCylinder(make_float3(0), make_float3(0, 0, 1), bin_diameter / 2., mat_type_terrain, 0); + walls->AddPlane(make_float3(0, 0, bottom), make_float3(0, 0, 1), mat_type_terrain); + + // The stock 3-sphere clump template, same numbers as DEMdemo_ConePenetration. + double clump_vol = 5.5886717; + float mass = terrain_density * clump_vol; + float3 MOI = make_float3(2.928, 2.6029, 3.9908) * terrain_density; + std::shared_ptr my_template = + DEMSim.LoadClumpType(mass, MOI, GetDEMEDataFile("clumps/3_clump.csv"), mat_type_terrain); + my_template->SetVolume(clump_vol); + my_template->Scale(clump_scale); + // The restart readers in stage 2 look clumps up by this name. + my_template->AssignName(template_name); + + // Loose-pour the fill column, then let it settle. + HCPSampler sampler(clump_scale * 3.); + float3 fill_center = make_float3(0, 0, bottom + 0.03 + fill_height / 2); + auto input_xyz = sampler.SampleCylinderZ(fill_center, bin_diameter / 2. - clump_scale * 2., fill_height / 2); + DEMSim.AddClumps(my_template, input_xyz); + std::cout << "Stage 1: poured " << input_xyz.size() << " grains" << std::endl; + + auto max_z_finder = DEMSim.CreateInspector("clump_max_z"); + auto total_mass_finder = DEMSim.CreateInspector("clump_mass"); + + DEMSim.SetInitTimeStep(step_size); + DEMSim.SetGravitationalAcceleration(make_float3(0, 0, -9.81)); + DEMSim.Initialize(); + + DEMSim.DoDynamicsThenSync(0.8); + + bed_surface = max_z_finder->GetValue(); + float bed_depth = bed_surface - bottom; + float bulk_density = total_mass_finder->GetValue() / (math_PI * bin_diameter * bin_diameter / 4. * bed_depth); + std::cout << "Stage 1: settled bed is " << bed_depth << " m deep, bulk density " << bulk_density << " kg/m^3" + << std::endl; + + DEMSim.WriteClumpFile(bed_file); + std::cout << "Stage 1: bed saved to " << bed_file << std::endl; + } + + // --------------------------------------------------------------------- + // Stage 2: reload the bed, switch cohesion on, press the plate. + // --------------------------------------------------------------------- + { + DEMSolver DEMSim; + DEMSim.SetVerbosity(INFO); + DEMSim.SetOutputFormat(OUTPUT_FORMAT::CSV); + DEMSim.SetOutputContent(OUTPUT_CONTENT::XYZ | OUTPUT_CONTENT::ABSV); + + // The stock cohesion force model. "Cohesion" is a pairwise material + // property; grain-grain gets the real value, plate-grain gets zero so + // the plate face does not adhere (set it nonzero for a sticky soil). + auto my_force_model = DEMSim.ReadContactForceModel("ForceModelWithCohesion.cu"); + my_force_model->SetMustHaveMatProp({"E", "nu", "CoR", "mu", "Crr", "Cohesion"}); + my_force_model->SetMustPairwiseMatProp({"CoR", "mu", "Crr", "Cohesion"}); + my_force_model->SetPerContactWildcards({"delta_time", "delta_tan_x", "delta_tan_y", "delta_tan_z"}); + + auto mat_type_terrain = DEMSim.LoadMaterial( + {{"E", 1e9}, {"nu", 0.3}, {"CoR", 0.4}, {"mu", 0.5}, {"Crr", 0.01}, {"Cohesion", cohesion}}); + auto mat_type_plate = + DEMSim.LoadMaterial({{"E", 1e9}, {"nu", 0.3}, {"CoR", 0.4}, {"mu", 0.5}, {"Crr", 0.01}, {"Cohesion", 0.}}); + DEMSim.SetMaterialPropertyPair("Cohesion", mat_type_plate, mat_type_terrain, 0.); + + const double world_size = 2; + DEMSim.InstructBoxDomainDimension(world_size, world_size, world_size); + DEMSim.InstructBoxDomainBoundingBC("top_open", mat_type_terrain); + auto walls = DEMSim.AddExternalObject(); + walls->AddCylinder(make_float3(0), make_float3(0, 0, 1), bin_diameter / 2., mat_type_terrain, 0); + walls->AddPlane(make_float3(0, 0, bottom), make_float3(0, 0, 1), mat_type_terrain); + + // Same template, same name, so the saved clumps map back onto it. + double clump_vol = 5.5886717; + float mass = terrain_density * clump_vol; + float3 MOI = make_float3(2.928, 2.6029, 3.9908) * terrain_density; + std::shared_ptr my_template = + DEMSim.LoadClumpType(mass, MOI, GetDEMEDataFile("clumps/3_clump.csv"), mat_type_terrain); + my_template->SetVolume(clump_vol); + my_template->Scale(clump_scale); + my_template->AssignName(template_name); + + // Reload the saved bed: positions AND orientations. The readers are + // keyed by clump-type name and return an empty vector for an unknown + // key, so guard against that rather than pressing a plate into vacuum. + auto xyz_map = DEMSim.ReadClumpXyzFromCsv(bed_file); + auto quat_map = DEMSim.ReadClumpQuatFromCsv(bed_file); + auto& in_xyz = xyz_map[template_name]; + auto& in_quat = quat_map[template_name]; + if (in_xyz.empty() || in_quat.size() != in_xyz.size()) { + std::cerr << "Reloading " << bed_file << " failed: got " << in_xyz.size() << " positions and " + << in_quat.size() << " orientations for clump type '" << template_name << "'." << std::endl; + return 1; + } + DEMClumpBatch bed_batch(in_xyz.size()); + bed_batch.SetTypes(std::vector>(in_xyz.size(), my_template)); + bed_batch.SetPos(in_xyz); + bed_batch.SetOriQ(in_quat); + bed_batch.SetFamily(0); + DEMSim.AddClumps(bed_batch); + std::cout << "Stage 2: reloaded " << in_xyz.size() << " grains from " << bed_file << std::endl; + + // The plate: the stock unit cylinder mesh (radius 1, height 2), + // scaled into a disc. Parked just above the stage-1 bed surface. + auto plate = DEMSim.AddWavefrontMeshObject(GetDEMEDataFile("mesh/cyl_r1_h2.obj"), mat_type_plate); + plate->Scale(make_float3(plate_diam / 2., plate_diam / 2., plate_thickness / 2.)); + float plate_mass = 7.8e3 * math_PI * (plate_diam / 2.) * (plate_diam / 2.) * plate_thickness; + plate->SetMass(plate_mass); + plate->SetMOI(make_float3(plate_mass * (3 * plate_diam * plate_diam / 4. + plate_thickness * plate_thickness) / 12., + plate_mass * (3 * plate_diam * plate_diam / 4. + plate_thickness * plate_thickness) / 12., + plate_mass * plate_diam * plate_diam / 8.)); + float park_z = bed_surface + plate_thickness / 2. + 0.01; + plate->SetInitPos(make_float3(0, 0, park_z)); + plate->SetFamily(2); + auto plate_tracker = DEMSim.Track(plate); + + // Family 2 holds still while the bed re-equilibrates under cohesion; + // family 1 descends at the press speed. The switch starts the test. + DEMSim.SetFamilyFixed(2); + DEMSim.SetFamilyPrescribedLinVel(1, "0", "0", "-" + to_string_with_precision(plate_speed)); + + auto max_z_finder = DEMSim.CreateInspector("clump_max_z"); + + DEMSim.SetInitTimeStep(step_size); + DEMSim.SetGravitationalAcceleration(make_float3(0, 0, -9.81)); + DEMSim.Initialize(); + + // Cohesion just switched on and tangential contact history is not + // carried across the two stages, so give the bed a moment to grip. + DEMSim.DoDynamicsThenSync(0.3); + bed_surface = max_z_finder->GetValue(); + std::cout << "Stage 2: bed re-equilibrated, surface at z = " << bed_surface << std::endl; + + DEMSim.ChangeFamily(2, 1); + + // March down, recording pressure against sinkage. The sinkage datum + // is first contact, taken as the depth where the mean plate pressure + // first reaches 1 kPa; before that the plate is falling through air. + const float plate_area = math_PI * (plate_diam / 2.) * (plate_diam / 2.); + const float contact_threshold = 1e3; // Pa + const float sample_every = 2e-3; // s + std::ofstream curve(curve_file); + curve << "sinkage_m,pressure_Pa" << std::endl; + double contact_z = 0.; + bool in_contact = false; + float next_print = 0.; + std::cout << "Stage 2: pressing at " << plate_speed << " m/s" << std::endl; + std::cout << " sinkage (mm) pressure (kPa)" << std::endl; + + while (true) { + DEMSim.DoDynamics(sample_every); + float3 plate_pos = plate_tracker->Pos(); + // ContactAcc returns the contact acceleration on the owner, so + // multiply by the mass we assigned to get the contact force. + float3 plate_force = plate_tracker->ContactAcc() * plate_mass; + float pressure = plate_force.z / plate_area; + + if (!in_contact) { + if (pressure >= contact_threshold) { + in_contact = true; + contact_z = plate_pos.z; + } + continue; + } + float sinkage = contact_z - plate_pos.z; + curve << sinkage << "," << pressure << std::endl; + if (sinkage * 1000. >= next_print) { + std::printf(" %6.1f %9.2f\n", sinkage * 1000., pressure / 1000.); + next_print += 5.; + } + if (sinkage >= target_sinkage) + break; + } + DEMSim.DoDynamicsThenSync(0.); + std::cout << "Stage 2: curve written to " << curve_file << std::endl; + DEMSim.ShowTimingStats(); + } + + std::cout << "DEMdemo_PlateSinkage exiting..." << std::endl; + return 0; +}