Help with multiple materials

Submitted by msandli on Fri, 09/21/2012 - 03:19

I'm trying to run a simulation that uses two different atom types for the particles and for the walls. My code was based on a code that took its material properties from the meshGran tutorial in LIGGGHTS - ie, small young's modulus, only one atom type. I could get that code to run to completion (this simulation would fill up my hopper with particles and then empty them out in 150000 time steps, roughly 5 minutes).

Now, I've put 2 atom types in the simulation, and put realistic values for the material properties. For some reason, the simulation will only run to approximately 500 timesteps, and then just seems to freeze. Will someone take a glance at my code and see if there's anything obvious? I've tried to keep everything labeled consistently, but I'm not sure if I've used the group command wrong, given something an improper group ID, or if the simulation is now so complex that it will simply take longer than I'm used to getting past that first 500 timesteps.

Please note that some commands are commented out - these are leftover from previous versions of this code, and help give me references of what things looked like when it worked.

Any help is greatly appreciated!

cstoltz | Fri, 09/21/2012 - 04:35

At a quick glance, your timestep is way too large to run in a stable fashion. Use the fix_check_timestep_gran command to check against the Rayleigh time. You're at something like 5-7*T_R, when a stable condition is usually more like 0.2*T_R.

Check when you expect the first contact to occur. I'm guessing it might be around the 500 timestep point.

Regards,
Chris

msandli | Fri, 09/21/2012 - 21:33

Chris,

I was able to get my timestep below .2*T_R (got rid of warning messege given by that fix). However, this still doesn't seem to fix it. I've looked at a series of dumped jpgs given by the code, and the simulation stops before the initial dump of particles is even out of the insertion zone. Using the timestep I've provided, the simulation stops at step 54500, but the first contact shouldn't be until around 500000 (or, using the old timestep, simulation stops at step 500, but first contact isn't until around 5000).

I'll edit my first post and add the updated code. I also changed some material properties because I was using shear modulus instead of young's modulus in the fix globals, but that didn't seem to do anything.

Any thoughts are appreciated. I'm running 2.1 if it makes a difference.

ckloss's picture

ckloss | Sun, 09/23/2012 - 12:34

Hi,

could you please perform the following steps: (1) apply your fix check/timestep/gran to "all" and (2) do a git pull update to 2.1.2 (3) ensure the insertion extrusion volume does not overlap any wall?
If the problem persists, I'll have a closer look!

Cheers, Christoph

msandli | Wed, 09/26/2012 - 17:03

Christoph,

I have done your 3 suggestions, and my simulation still pauses at the same time step. I am getting a new warning message that I haven't been getting in older versions: WARNING: Mesh contains highly skewed element (multi_node_mesh_parallel_I.h:513), and there is a line in my output file saying that I have 16 mesh elements with high aspect ratio, but otherwise, my simulation behaves exactly as it did with 2.1.

I will edit my first post and upload my current input script (liggghts_test_3), an output file (standard_out_file...), and an .stl dump (dump30000.txt). I don't think there's any overlapping between my extrusion volume (which is extended in the positive Z direction), and any part of the mesh that I want to use as a real wall, or even my region boundaries.

Your help is greatly appreciated!

ckloss's picture

ckloss | Fri, 09/28/2012 - 14:15

You didn't attach all the stl files. Please, if you post cases for debugging, attach all files needed to run them

Christoph

msandli | Fri, 09/28/2012 - 15:39

Sorry about that. I've added the .stl files (as .txt files) to my original post. hopper_open_2 is my hopper, hopper_insert_2 is my insertion face, and valve is my "valve" - basically a rectangle at the bottom of the hopper. Liggghts_test_3 is my current input script. Again, I'm just trying to dump some particles into a hopper, then open the "valve" and let them drain out.

Thanks again!

ckloss's picture

ckloss | Tue, 10/02/2012 - 21:04

Hi,

this was actually really a bug - will be fixed in the next release! In the meantime, you could (as a quickfix) replace the lines

if(! (mask[iAtom] & groupbit))
continue;

by
if(! (mask[iAtom] & groupbit))
{
if(bins) iAtom = bins[iAtom];
else iAtom = -1;
continue;
}

in fix_neighlist_mesh.cpp

Cheers, Christoph

msandli | Wed, 10/03/2012 - 17:33

Christoph,

Changed code as described. Appears to have no effect on my simulations - they still stop at the same timestep. This code appears at line 196, correct? Should I delete where this added code appears a few lines later (which I didn't do, but was just curious. Whenever I have to replace some code in something I've written with a line(s) that appears elsewhere, it usually means the whole thing is in the wrong place)?

Otherwise, I will wait until the next release.

Thank you so much for the help!

cstoltz | Thu, 10/04/2012 - 13:17

I tried running your code and it ran fine on my machine. It didn't stall. I did notice that the simulation is still unstable and many particles are lost. If you lower the modulus to something more reasonable from a simulation standpoint, e.g. O(10^7), it runs fine with no particles lost. Not sure why it would stall on your machine though.

Regards,
Chris

msandli | Thu, 10/04/2012 - 16:27

Chris,

My simulation uses glass and steel as the two materials, and they have Young's modulus in the range of GPa. Assuming the unit used for Young's modulus in the fix command is Pa (for si units), then my material attributes should be accurate. I'm comparing the results of LIGGGHTS with both the results of another DEM program, and a small lab scale experiment, so I need those parameters as close to "reality" as I can get them.

Even when I had two materials with relatively low Young's modulus (10^6), my program would still freeze at the same spot. One material is fine, seemingly without regard to material properties.

cstoltz | Thu, 10/04/2012 - 18:34

Agreed that the Young's modulus should be in the range of GPa for the true materials. That said, for your simulation (which I'm assuming will be just emptying a hopper to understand flow rates, segregation, etc.), I think that you'll find that the results are not highly sensitive to the Young's modulus provided you don't use a value that is too low. I've found in the past that for hopper simulations, which are dominated by friction, using a modulus of O(10^7) is usually more than sufficient, and will allow you to run much larger simulations in a reasonable amount of time. For other simulations, such as a high speed mixer, getting the modulus right can be of greater importance.

I should note that your x_max domain coordinate is coincident with the wall of your geometry. Just to be on the safe side, I would expand your region slightly to ensure that your wall stays within the domain.

Regards,
Chris

msandli | Mon, 10/01/2012 - 18:03

I've been trying to figure out where the problem is, and I'm still stumped. I took my "simple" simulation - same idea, but with a single material with young's mod around 10^6 - and made the material more realistic (young's mod 10^10), and the simulation ran fine. But whenever I try to use two materials, even if I make them have identical properties, my simulation still freezes with a group of particles hanging in midair after falling a bit from the insertion region.

So I suspect I've done some incorrect coding in how I set up the two different groups, but everything in my script looks right to me, and the fact that the simulation inserts a group of particles and then actually runs for a bit before freezing (without any interaction between the groups) is very puzzling to me.

Any ideas for me to try? Again, all of my code should be appended to my first post.