Hello people. I'm working on a simulation of a tensile test with granular bonds. I'm currently using liggghts 3.5.0 with the bondspackage uploaded by richtie. The thing is that I'm trying to get a tension-deformation diagram for the case, but the force output by the command compute/reduce is just too low. The bonds are working as intended, as they break when they should break (the simulation uses a velocity in one end, while the other end is freezed, so the deformation is directly proportional to the time elapsed), but as i said, the force is orders of magnitude wrong. So the question, is there a problem in the code, or should i compute the force in some other way? what could be that way?
Also I've seen in the example that you can use the "numbond" on thermo output, but it doesnt seem to work in 3.5.0 version. Is there any file that i need to modify to get that?
Thanks
| Attachment | Size |
|---|---|
 simple.txt | 2.82 KB |
mschramm | Wed, 05/10/2017 - 15:55
I believe
I believe this is what you want
compute b1 all property/local batom1x batom1y batom1z batom2x batom2y batom2z batom1 batom2 btype bforceX bforceY bforceZ
as for numbond, is it declared in thermo.cpp?
felipe.merino | Wed, 05/10/2017 - 20:19
Yes
Yes, that is what im looking for, but the args batom1x, batom1y, etc, and bforceX, bforceY and bforceZ, are not recognized by the compute. Do i have to define them before the compute?
mschramm | Wed, 05/10/2017 - 21:17
The Joy of 3.3 to 3.6
Check your compute_property_local.cpp file.
You should see the following lines
} else if (strcmp(arg[iarg],"btype") == 0) {
pack_choice[i] = &ComputePropertyLocal::pack_btype;
if (kindflag != NONE && kindflag != BOND)
error->all(FLERR,
"Compute property/local cannot use these inputs together");
kindflag = BOND;
} else if (strcmp(arg[iarg],"batom1x") == 0)
pack_choice[i] = &ComputePropertyLocal::pack_batom1x;
else if (strcmp(arg[iarg],"batom1y") == 0)
pack_choice[i] = &ComputePropertyLocal::pack_batom1y;
else if (strcmp(arg[iarg],"batom1z") == 0)
pack_choice[i] = &ComputePropertyLocal::pack_batom1z;
else if (strcmp(arg[iarg],"batom2x") == 0)
pack_choice[i] = &ComputePropertyLocal::pack_batom2x;
else if (strcmp(arg[iarg],"batom2y") == 0)
pack_choice[i] = &ComputePropertyLocal::pack_batom2y;
else if (strcmp(arg[iarg],"batom2z") == 0)
pack_choice[i] = &ComputePropertyLocal::pack_batom2z;
else if (strcmp(arg[iarg],"bforceX") == 0)
pack_choice[i] = &ComputePropertyLocal::pack_bforceX;
else if (strcmp(arg[iarg],"bforceY") == 0)
pack_choice[i] = &ComputePropertyLocal::pack_bforceY;
else if (strcmp(arg[iarg],"bforceZ") == 0)
pack_choice[i] = &ComputePropertyLocal::pack_bforceZ;
I had to add these lines to my 3.6 release of liggghts.
felipe.merino | Wed, 05/31/2017 - 17:55
Still the same problem
Ok, so i tried what yo told me. Sadly, i found that the value dumped by compute property/local is the same as the value calculated by the program (the one that shows the regular dump), which is far below the theorical force on the bond. For instance, im getting a force value of 7.84E-06 N when it should be about 31 N. Im using a simulation of a simple bond, with 2 particles and using only axial displacement.
felipe.merino | Wed, 05/31/2017 - 17:57
I have uploaded the input
I have uploaded the input script, in case anyone can look it up.
mschramm | Thu, 06/01/2017 - 20:04
Try This
In the call for bond parameters, ${r} is a constant to the bond radius
rbmin=rb[type]*MIN(radius[i1],radius[i2]); //lamda * min(rA,rB) see P.Cundall, "A bondet particle model for rock"
try using
bond_coeff 1 1.0 ${kn} ${kt} 1 ${sn} ${st} #!!!
felipe.merino | Fri, 06/02/2017 - 17:01
I've noticed that
Yeah, i did read the code and the paper. i have also tried both values, as a constant, around 1.0, and as the radius. With the values i used, considering a radius of the bond equal to the particle, the bond breakage should be at timestep 6324. Using the constant equal to the radius (0.0005 m) im getting the breakage at step 6328. Using a value higher than 0.5, the bond doesnt break. However, in both cases, the force value is not even close to 31 N. I could ignore the force value to get the simulation running, but cannot validate the model, or even get a stress-deformation curve.
mschramm | Fri, 06/02/2017 - 20:26
Running out of Ideas...
I'm running out of things to try.
Yu Guo has done work with this model (Validation and time step determination of discrete element modeling of flexible fibers)
Your time step is below the threshold.
Only thing coming to mind right now is to try more particles. Maybe the fixes at the end are messing with the calculation.
saiya | Sat, 05/02/2020 - 10:24
HOW TO INTRODUCE VISCOSITY
Hello, I want to do a viscous simulation.
But I can't do that simulation.
I think because it is that I can't introduce a new package to use viscous.
Please, tell me how to introduce a package to use viscous.