DEM-Engine geometry motion setting

[simboys]

Hi,

     I want to achieve a composite rotational motion of a geometric body, as follows

Motion 1: Rotate at a speed of w1 rad/s, with the rotation axis being the z-axis
Motion 2: Rotate at a speed of w2 rad/s, with a rotation vector of [1 * cos (w1 * time), 1 * sin (w1 * time), 0], and an action point of [0, 0, 5], as shown in attached figure.

However, the following issues were encountered:
1. After using DEMTracker to track geometry, even if I don't add any motion, the geometry will continue to drop down. How can I fix the geometry?
2. How to achieve the composite rotational motion mentioned above?

Best regards,

Wenxuan Xu

[Ruochun Zhang]

Hi Wenxuan,

1. When you load a mesh or an analytical object, it's by default in family 255, a reserved and automatically fixed family. You probably changed your mesh into a family that does not have any motion prescription, then in that case it follows the "normal physics" like all the particles you load in, including falling under gravity. If you wish to fix that family, just use SetFamilyFixed.

2. This is more complicated. First, we have to understand that in DEME, strictly speaking, objects work with their principal frames of inertia (centered at CoM, having a diagonal inertia matrix). This is a current limitation due to how the inertia matrix-defining APIs are designed: you can only input the principal moments of inertia. But this also means that setting the frame correctly is not necessary if you completely prescribe the motion of that object (a.k.a. this object does not respond to physics), and you don't care about getting a torque reading. However, in your case, I suggest you do so, since you'll prescribe a complex motion, and standardizing the frame to use at least the CoM frame can save you some head-scratching. I suggest that you make it such that in the mesh file you load in, the (0,0,0) point is right at the CoM of this mesh; but if not, you can adjust the CoM and principal axes with InformCentroidPrincipal after loading in the mesh.

After that is taken care of, you need to combine the two rotations into a rotation about this object's own CoM and a linear translation of this object's CoM (because family motion prescription allows for setting these two motions only). Use your MBD knowledge to do that. I don't guarantee the correctness of this, but some back-of-the-envelop calculation seems to put the angular velocity to be (w2×cos(w1×t), w2×sin(w1×t), w1), and the linear velocity to be 5×(w2×cos(w2×t)×sin(w1×t)+sin(w2×t)×w1×cos(w1×t), -w2×cos(w2×t)×cos(w1×t)+sin(w2×t)×w1×sin(w1×t), 0). Again, please verify yourself. After that, prescribe the object's family's motion using it.

But if this is too complicated and you can write your code to calculate the location and quaternion of this object at any given t (which for some people is easier than mathematically deriving the analytical expression for the linear and rotational motions), then you can set this object fixed, then use the tracker methods SetPos and SetOriQ to enforce the location and rotation of this object at each time step. It should achieve the same assuming this object is not moving at a very high velocity.

Thank you,

Ruochun

[simboys]

Hi, Ruochun

       Thank you for your reply, I wiil try it. By the way, can you share me the method how to debug the code step-by-step in linux? In addition, if I use the DEMTracker to set the location and rotation of geometry at each time step, will this affect computational efficiency? And how to set the location of geometry at each time step? From the demo, DoDynamics(frameTime) is many times longer than timestep, and if we set the frameTime by timestep, it will result in a huge data file.

Best regards,

Wenxuan

[simboys]

Hi, Ruochun

       I have tried it, but the geometric motion doesn't seem to achieve the expected result. Can you help me troubleshoot the error? Ths mesh obj file are attached.

#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(STEP_METRIC);

// For general use cases, you want to set the verbosity to INFO: It's also a bit faster than STEP_METRIC.

// DEMSim.SetVerbosity(INFO);

DEMSim.SetOutputFormat(OUTPUT_FORMAT::CSV);

DEMSim.SetOutputContent(OUTPUT_CONTENT::ABSV);

DEMSim.SetMeshOutputFormat(MESH_FORMAT::VTK);

// If you don't need individual force information, then this option makes the solver run a bit faster.

DEMSim.SetNoForceRecord();

// E, nu, CoR, mu, Crr...

auto mat_type_mixer = DEMSim.LoadMaterial({{"E", 1e8}, {"nu", 0.3}, {"CoR", 0.6}, {"mu", 0.5}, {"Crr", 0.0}});

auto mat_type_granular = DEMSim.LoadMaterial({{"E", 1e8}, {"nu", 0.3}, {"CoR", 0.6}, {"mu", 0.2}, {"Crr", 0.0}});

// If you don't have this line, then mu between mixer material and granular material will be 0.35 (average of the

// two).

DEMSim.SetMaterialPropertyPair("mu", mat_type_mixer, mat_type_granular, 0.5);

float step_size = 5e-5;

const double world_size = 10;

const float chamber_height = world_size / 3.;

const float fill_height = chamber_height;

const float chamber_bottom = -world_size / 2.;

const float fill_bottom = chamber_bottom + chamber_height;

DEMSim.InstructBoxDomainDimension(world_size, world_size, world_size);

DEMSim.InstructBoxDomainBoundingBC("all", mat_type_granular);

auto mixer = DEMSim.AddWavefrontMeshObject((GET_DATA_PATH() / "mesh/Chute_Rotating.obj").string(), mat_type_mixer);

std::cout << "Total num of triangles: " << mixer->GetNumTriangles() << std::endl;

mixer->SetFamily(10);

mixer->InformCentroidPrincipal(make_float3(0, 0, 6.793), make_float4(0, 0, 0, 1));

double w1 = 3.14159;

double w2 = 1.0;

std::string angVel_x = "fabs(1.0 * cos(deme::PI * t))";

std::string angVel_y = "fabs(1.0 * sin(deme::PI * t))";

std::string angVel_z = "deme::PI";

std::string vel_x = "(1.0 * cos(1.0 * t) * sin(deme::PI * t) + sin(1.0 * t) * deme::PI * cos(deme::PI * t))";

std::string vel_y = "(-1.0 * cos(1.0 * t) * cos(deme::PI * t) + sin(1.0 * t) * deme::PI * sin(deme::PI * t))";

std::string vel_z = "0";

DEMSim.SetFamilyPrescribedAngVel(10, angVel_x, angVel_y, angVel_z);

DEMSim.SetFamilyPrescribedLinVel(10, vel_x, vel_y, vel_z);

float granular_rad = 0.005;

// Calculate its mass and MOI

float mass = 2.6e3 * 5.5886717; // in kg or g

float3 MOI = make_float3(2.928, 2.6029, 3.9908) * 2.6e3;

std::shared_ptr<DEMClumpTemplate> template_granular = DEMSim.LoadClumpType(mass, MOI, GetDEMEDataFile("clumps/3_clump.csv"), mat_type_granular);

template_granular->Scale(granular_rad);

// Track the mixer

auto mixer_tracker = DEMSim.Track(mixer);

DEMSim.SetInitTimeStep(step_size);

DEMSim.SetGravitationalAcceleration(make_float3(0, 0, -9.81));

// Initialize the simulation system

DEMSim.Initialize();

path out_dir = current_path();

out_dir += "/DEMdemo_RotatingChute";

create_directory(out_dir);

float sim_end = 2.0;

unsigned int fps = 20;

float frame_time = 1.0 / fps;

std::cout << "Output at " << fps << " FPS" << std::endl;

unsigned int currframe = 0;

// mixer_tracker->SetOriQ(make_float4(0, 0, 0, 1));

for (float t = 0; t < sim_end; t += frame_time) // 0.05迭代一次

{

std::cout << "Position: " << mixer_tracker->GetPos()[0] << " " << mixer_tracker->GetPos()[1] << " " << mixer_tracker->GetPos()[2] << std::endl;

std::vector<float> Qua = mixer_tracker->GetOriQ(); // x, y, z, w

// std::cout << Qua[0] << " " << Qua[1] << " " << Qua[2] << " " << Qua[3] << std::endl;

std::cout << "Frame: " << currframe << std::endl;

char filename[200], meshfilename[200], cnt_filename[200];

sprintf(filename, "%s/DEMdemo_output_%04d.csv", out_dir.c_str(), currframe);

sprintf(meshfilename, "%s/DEMdemo_mesh_%04d.vtk", out_dir.c_str(), currframe);

sprintf(cnt_filename, "%s/Contact_pairs_%04d.csv", out_dir.c_str(), currframe++);

DEMSim.WriteMeshFile(std::string(meshfilename));

DEMSim.DoDynamics(frame_time);

}

std::cout << "DEMdemo_RotatingChute simulating complete!" << std::endl;

return 0;

}

[Ruochun Zhang]

Hi Wenxuan,

Regarding debugging on Linux, I recommend seeking professional assistance. However, I find that VSCode's GUI for stepping through code is sufficient for my needs.

For enforcing motion on an entity step-by-step, use DoDynamics(step_size) or DoStepDynamics() instead of DoDynamics(frame_duration). This way, you can count the number of steps and use if statements to control outputs. The RotatingDrum demo provides a good example of this.

The performance impact should be minimal, as DEME is designed for fine-grain operations like these to enable more flexible co-simulation. The two GPU streams synchronize only with DoDynamicsThenSync(), not DoDynamics(). Also, transferring the explicit position and velocity of one or two entities to GPUs at each step introduces only a small and manageable overhead, too.

Regarding your script, it doesn't fully implement my specified formula. I didn't use fabs, and the linear velocity should have a scaling factor of 5, which seems to be missing. Even with my specified formula, I can't guarantee its correctness (without knowing details like what the coordinate system is, and what the "expected results" is). But more realistically, it's that I cannot do your research for you. You see, this is becoming more of a consulting task than troubleshooting.

I suggest starting simple, using only one of the implemented motions, verifying it works as expected, and then gradually adding components of the second one to approach your desired outcome.

Thank you,

Ruochun