Unbalanced force problem with clump

Hello Zhaochen,

could you please provide us with the script you use? It will help us to
find the problem and solve your question.

How many bodies do you have in your simulation? Just the one clump?

Thanks
Jan
Dne 7.12.2012 2:51 "kelaogui" <email address hidden>
napsal(a):

> New question #216150 on Yade:
> https://answers.launchpad.net/yade/+question/216150
>
> Dear sir
> I am a fresh man of Yade. Now I face a problem in triaxial test
> with clump. I use ThreeDTriaxialEngine and move the wall to do the
> confinement. But during the confining process, the unbalanced force is
> rather high, about 10e3 magnitude. My clump is the most simple one which is
> composed by two spheres with equal radius. I want to know where the problem
> comes from?
>
> Thank you
> Zhaochen
>
> --
> You received this question notification because you are a member of
> yade-users, which is an answer contact for Yade.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-users
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>

Yeah,
    Thank you for your suggestion. Here is my code.

from yade import pack
from yade import *

num_spheres=2000
## corners of the initial packing
mn,mx=Vector3(0,0,0),Vector3(1,1,1)
thick = 0.01
compFricDegree = 2
rate=0.2
damp=0.1
stabilityThreshold=0.001
key='_define_a_name_'

## create material #0, which will be used as default
O.materials.append(FrictMat(young=5e6,poisson=0.5,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
O.materials.append(FrictMat(young=5e6,poisson=0.5,frictionAngle=0,density=0,label='walls'))

## create walls around the packing
walls=utils.aabbWalls([mn,mx],thickness=thick,material='walls')
wallIds=O.bodies.append(walls)

sp=pack.SpherePack()
sp.makeCloud(mn,mx,-1,0.3333,num_spheres,False, 0.95)

volume = (mx[0]-mn[0])*(mx[1]-mn[1])*(mx[2]-mn[2])
mean_rad = pow(0.09*volume/num_spheres,0.3333)

clumps=True
if clumps:
 c1=pack.SpherePack([((-0.2*mean_rad,0,0),0.5*mean_rad),((0.2*mean_rad,0,0),0.5*mean_rad)])
 sp.makeClumpCloud((0,0,0),(1,1,1),[c1],periodic=False)
 O.bodies.append([utils.sphere(center,rad,material='spheres') for center,rad in sp])
 standalone,clumps=sp.getClumps()
 for clump in clumps:
  O.bodies.clump(clump)
  for i in clump[1:]: O.bodies[i].shape.color=O.bodies[clump[0]].shape.color
 #sp.toSimulation()
else:
 O.bodies.append([utils.sphere(center,rad,material='spheres') for center,rad in sp])

O.dt=.5*utils.PWaveTimeStep() # initial timestep, to not explode right away
O.usesTimeStepper=True
triax=ThreeDTriaxialEngine(
 #maxMultiplier=1.005,
 #finalMaxMultiplier=1.002,
 thickness = thick,
 stressControl_1 = True,
 stressControl_2 = True,
 stressControl_3 = True,
        strainRate1=0.01,
        strainRate2=0.01,
        strainRate3=0.01,
 ## The stress used for (isotropic) internal compaction
 sigma_iso = 10000,
 ## Independant stress values for anisotropic loadings
 sigma1=10000,
 sigma2=10000,
 sigma3=10000,
 internalCompaction=False,
 Key=key,
)

newton=NewtonIntegrator(damping=damp)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()],verletDist=-mean_rad*0.06),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack(label='cundall')]
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8, defaultDt=4*utils.PWaveTimeStep()),
 triax,
 TriaxialStateRecorder(iterPeriod=100,file='WallStresses'+key),
 newton
]

#Display spheres with 2 colors for seeing rotations better
Gl1_Sphere.stripes=0
yade.qt.Controller(), yade.qt.View()

while 1:
  O.run(100, True)
  #the global unbalanced force on dynamic bodies, thus excluding boundaries, which are not at equilibrium
  unb=unbalancedForce()
  meanS=(triax.stress(triax.wall_right_id)[0]+triax.stress(triax.wall_top_id)[1]+triax.stress(triax.wall_front_id)[2])/3
  sidelen1=O.bodies[1].state.se3[0][0]-O.bodies[0].state.se3[0][0]-thick
  print 'unbalanced force:',unb,' mean stress: ',meanS,'side: ',sidelen1
  if unb<stabilityThreshold and abs(meanS-triax.sigma_iso)/triax.sigma_iso<0.001:
    print "### Isotropic state saved ###"
    break

You script is in fact the triax-tutorial script. :)
I must admit it has this block for clumps which does not work properly, so if you fix it let me know, I'll update it.

What version of yade are you using? For older versions unbalancedForce did not work.

Yeah, it's triax-tutorial script. The version is the newest trunk version. So, where can i find the code for computing unbalanced force?

In latest versions, unbalancedForce is ok for clumps.
So the problem is the clump generation part of triax-tutorial.py, which is broken (I made an update of the triax tutorial today and removed the clump part).
Actually, if you look what happens in the first steps, you will see some particles flying away from the box, this is because initial positions are not consistent with the size of the problem. I don't have a fixed version of this script personaly. You may try and fix it yourself, it should not be really difficult. Let me know if you find a fix please, so that I will fix the original tutorial script.

I am sorry but the script you offered will result in this error:

Traceback (most recent call last):
  File "/usr/local/bin/yade", line 144, in runScript
    execfile(script,globals())
  File "triaxial clump.py", line 35, in <module>
    sp.toSimulation(material='spheres')
  File "/home/dylan/Yade/lib/yade-2012-11-01.git-ce12d2b/py/yade/pack.py", line 112, in SpherePack_toSimulation
    clumpIds.append(O.bodies.clump(clump)) # clump spheres with given ids together, creating the clump object as well
RuntimeError: Clump::updateProperties: NaNs in eigen-decomposition of inertia matrix?!

It seems the material parameter did not pass to the clump? Is it because of the version of trunk?

It's OK now.
It's the problem of my version. Thanks, Chareyre.

Thanks Chareyre, that solved my question.