How to quickly reduce unbalancedForce to approximate 0?

Dear all,

I am using the PotentialBlock function to do a regular dodecahedron stacking experiment, but I found that after many time steps unbalancedForce still fluctuates around 1, and the dodecahedron is still shaking. Is there any way to quickly reduce the unbalancedForce or let the polyhedron tend to be stationary and stable.

Jie

Question information

Language:
English Edit question
Status:
Solved
For:
Assignee:
No assignee Edit question
Solved by:
Vasileios Angelidakis
Solved:
2020-05-18
Last query:
2020-05-18
2020-05-12

This question was reopened

 Vasileios Angelidakis (vsangelidakis) said on 2020-05-08: #1

Hi Jie,

Depends on what stage of your analysis you are.

If you just want to stabilise your sample/particle during initialisation, you can use local damping (in the NewtonIntegrator), which is un-physical but will do the job.

If you want to keep some damping during the main stage of your main simulation, I would advise to not use local damping, but use a small amount of viscousDamping in the KnKsPBLaw instead.

If now you want to stabilise your system during initialisation momentarily, you can invoke "calm()" for some timesteps, to reset the velocity of the particles.

Hope this helps,
Vasileios

 weijie (amandajoe) said on 2020-05-09: #2

Hi Vasileios, and thank you again.

What is the value range of viscousDamping?

Best regards,
Jie

 Vasileios Angelidakis (vsangelidakis) said on 2020-05-09: #3

Hi Jie,

The viscousDamping parameter takes values from 0 to 1.0 (it is a percentage) and currently it only applies damping on the normalForce (the shearForce remains undamped for the time being).

Regards,
Vasileios

 weijie (amandajoe) said on 2020-05-10: #4

Thanks Vasileios Angelidakis, that solved my question.

 weijie (amandajoe) said on 2020-05-11: #5

Hi Vasileios, and thank you again.

If I use the calm () function to reset the velocity of the particles, the program appears as follows: shearForce: -nan -nan -nan, normalForce: -0.65154 1.59526 -0.0581212, viscousNormal: -0 0 -0, viscousShear: 0 0 0, geom normal: -0.377888 0.925237 -0.0337099, effective_phi: 1718.87, shearIncrement: 0 0 0, cs: 4.30304, incidentVs: 0 0 0, id1: 865, id2: 1562, debugShear: -0 0 0, phys-> mobilizedShear : -nan.

Some particles disappear, and the non-disappearing particles will then pass through the surface of the container. How to solve this problem?

Best regards,
Jie

 Vasileios Angelidakis (vsangelidakis) said on 2020-05-11: #6

Hi Jie,

This is strange. Calm [1] only resets the translational and angular velocity (and the angular momentum), it does not reset forces.
The warning you get comes from [2], which tells you that you have no contact forces (hence the particles pass through boundaries) but I cannot understand why this happens. "Calm" alone cannot cause this behaviour. In order to be more helpful, I need a MWS to reproduce this.

All the best,
Vasileios

 weijie (amandajoe) said on 2020-05-12: #7

Hi Vasileios, and thank you again.

I shortened the code, leaving the part where the main problem occurred. At 21,000 steps, I used the calm () function and the above problem occurred. If I don't use the calm () function, I find that after many steps, the entire system is unstable, unbalancedForce fluctuates around 2, and there is no tendency to decrease, which is obviously inconsistent with the actual reality. Can you help by the way to see what is causing this?

Best regards,
Jie

Following is my codeļ¼
#################

import numpy as np
import math
import random
import os
import errno
try:
os.mkdir('./vtk/')
except OSError as exc:
if exc.errno != errno.EEXIST:
raise
pass

O.engines=[
ForceResetter(),
InteractionLoop(
[Ig2_PB_PB_ScGeom(twoDimension=False, unitWidth2D=1.0, calContactArea=True)],
[Ip2_FrictMat_FrictMat_KnKsPBPhys(kn_i=1e8, ks_i=1e8, Knormal=1e8, Kshear=1e8, useFaceProperties=False, viscousDamping=0.2)],
[Law2_SCG_KnKsPBPhys_KnKsPBLaw(label='law',neverErase=False, allowViscousAttraction=True, traceEnergy=False)]
),
#GlobalStiffnessTimeStepper(),
NewtonIntegrator(damping=0.0,exactAsphericalRot=True,gravity=[0,-9.81,0]),
PotentialBlockVTKRecorder(fileName='./vtk/cubePBscaled',iterPeriod=300000,twoDimension=False,sampleX=50,sampleY=50,sampleZ=50,maxDimension=0.2,label='vtkRecorder')
]

powderDensity = 3000
O.materials.append(FrictMat(young=-1,poisson=-1,frictionAngle=30.0,density=powderDensity,label='frictionless'))
#######################

meanSize = 0.005
wallThickness = 0.5*meanSize
lengthOfBase = 0.250
heightOfBase = 0.600

############################################

sp=pack.SpherePack()
mn,mx=Vector3(-0.5*(lengthOfBase-wallThickness),0.5*meanSize,-0.5*(lengthOfBase-wallThickness)),Vector3(0.5*(lengthOfBase-wallThickness),0.5*heightOfBase,0.5*(lengthOfBase-wallThickness))#
R=sqrt(3.0)*0.009654
sp.makeCloud(mn,mx,R,0,-1,False)

r=0.01*meanSize
aa=[0.5257311121191336, -5.836787854364519e-17, 0.5257311121191336, -5.836787854364519e-17, 0.85065080835204, -0.5257311121191336, 0.85065080835204, -0.5257311121191336, -0.85065080835204, 5.836787854364519e-17, 5.836787854364519e-17, -0.8506508083520399]
bb=[0.85065080835204, -0.5257311121191336, -0.85065080835204, 0.5257311121191336, 5.836787854364519e-17, 0.85065080835204, 5.836787854364519e-17, -0.85065080835204, 5.836787854364519e-17, 0.5257311121191336, -0.5257311121191336, -0.0]
cc=[-5.836787854364519e-17, 0.85065080835204, 5.836787854364519e-17, 0.85065080835204, -0.5257311121191336, 5.836787854364519e-17, 0.5257311121191336, 5.836787854364519e-17, 0.5257311121191336, -0.85065080835204, -0.85065080835204, -0.5257311121191336]
dd=[0.01074978698202965, 0.010749786982029651, 0.01074978698202965, 0.01074978698202965, 0.01074978698202965, 0.01074978698202965, 0.010749786982029651, 0.010749786982029651, 0.01074978698202965, 0.01074978698202965, 0.010749786982029651, 0.01074978698202965]
for s in sp:
b=Body()
b.aspherical=True
wire=False
color=Vector3(random.random(),random.random(),random.random())
highlight=False
b.shape=PotentialBlock(k=0.0, r=r, R=0.0, a=aa, b=bb, c=cc, d=dd, isBoundary=False, color=color,
wire=wire, highlight=highlight, minAabb=Vector3(R,R,R), maxAabb=Vector3(R,R,R), AabbMinMax=True, fixedNormal=False, id=len(O.bodies))
utils._commonBodySetup(b, b.shape.volume, b.shape.inertia, material='frictionless',pos=s[0], fixed=False)
b.state.pos = s[0] #s[0] stores center
b.state.ori = Quaternion((random.random(),random.random(),random.random()),random.random()) #s[2]
O.bodies.append(b)
#######################################################

r=0.1*wallThickness
bbb=Body()
wire=False
color=[0,0.5,1]
highlight=False
bbb.shape=PotentialBlock(k=0.0, r=r, R=0.0, a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1], d=[2.5*lengthOfBase-r,2.5*lengthOfBase-r,0.5*wallThickness-r,0.5*wallThickness-r,lengthOfBase/2.0-r,lengthOfBase/2.0-r], id=len(O.bodies), isBoundary=True, color=color, wire=wire, highlight=highlight, AabbMinMax=True, minAabb=1.05*Vector3(lengthOfBase*2.5,0.5*wallThickness,lengthOfBase/2.0), maxAabb=1.05*Vector3(lengthOfBase*2.5,0.5*wallThickness,lengthOfBase/2.0), fixedNormal=False)
utils._commonBodySetup(bbb, bbb.shape.volume, bbb.shape.inertia, material='frictionless', pos=[0,0,0], fixed=True)
bbb.state.pos = [2*lengthOfBase,0,0]
O.bodies.append(bbb)

#
bA=Body()
wire=False
color=[0,0.5,1]
highlight=False
bA.shape=PotentialBlock(k=0.0, r=r, R=0.0, a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1], d=[0.5*wallThickness-r,0.5*wallThickness-r,0.5*heightOfBase-r,0.5*heightOfBase-r,0.5*lengthOfBase-r,0.5*lengthOfBase-r], id=len(O.bodies), isBoundary=True, color=color, wire=wire, highlight=highlight, AabbMinMax=True, minAabb=1.05*Vector3(0.4*wallThickness,0.5*heightOfBase,0.5*lengthOfBase), maxAabb=1.05*Vector3(0.4*wallThickness,0.5*heightOfBase,0.5*lengthOfBase), fixedNormal=False)
utils._commonBodySetup(bA, bA.shape.volume, bA.shape.inertia, material='frictionless', pos=[0,0,0], fixed=True)
bA.state.pos = [0.5*lengthOfBase,0.5*heightOfBase,0]
O.bodies.append(bA)

bB=Body()
wire=False
color=[0,0.5,1]
highlight=False
bB.shape=PotentialBlock(k=0.0, r=r, R=0.0, a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1], d=[0.5*wallThickness-r,0.5*wallThickness-r,0.5*heightOfBase-r,0.5*heightOfBase-r,0.5*lengthOfBase-r,0.5*lengthOfBase-r], id=len(O.bodies), isBoundary=True, color=color, wire=wire, highlight=highlight, AabbMinMax=True, minAabb=1.05*Vector3(0.4*wallThickness,0.5*heightOfBase,0.5*lengthOfBase), maxAabb=1.05*Vector3(0.4*wallThickness,0.5*heightOfBase,0.5*lengthOfBase), fixedNormal=False)
utils._commonBodySetup(bB, bB.shape.volume, bB.shape.inertia, material='frictionless', pos=[0,0,0], fixed=True)
bB.state.pos = [-0.5*lengthOfBase,0.5*heightOfBase,0]
O.bodies.append(bB)

bC=Body()
wire=True
color=[0,0.5,1]
highlight=False
bC.shape=PotentialBlock(k=0.0, r=r, R=0.0, a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1], d=[2.5*lengthOfBase-r,2.5*lengthOfBase-r,0.5*heightOfBase-r,0.5*heightOfBase-r,0.5*wallThickness-r,0.5*wallThickness-r], id=len(O.bodies), isBoundary=True, color=color, wire=wire, highlight=highlight, AabbMinMax=True, minAabb=1.05*Vector3(2.5*lengthOfBase,0.5*heightOfBase,0.4*wallThickness), maxAabb=1.05*Vector3(2.5*lengthOfBase,0.5*heightOfBase,0.4*wallThickness), fixedNormal=False)
utils._commonBodySetup(bC, bC.shape.volume, bC.shape.inertia, material='frictionless', pos=[0,0,0], fixed=True)
bC.state.pos = [2*lengthOfBase,0.5*heightOfBase,0.5*lengthOfBase]
O.bodies.append(bC)

bD=Body()
wire=False
color=[0,0.5,1]
highlight=False
bD.shape=PotentialBlock(k=0.0, r=r, R=0.0, a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1], d=[2.5*lengthOfBase-r,2.5*lengthOfBase-r,0.5*heightOfBase-r,0.5*heightOfBase-r,0.5*wallThickness-r,0.5*wallThickness-r], id=len(O.bodies), isBoundary=True, color=color, wire=wire, highlight=highlight, AabbMinMax=True, minAabb=1.05*Vector3(2.5*lengthOfBase,0.5*heightOfBase,0.4*wallThickness), maxAabb=1.05*Vector3(2.5*lengthOfBase,0.5*heightOfBase,0.4*wallThickness), fixedNormal=False)
utils._commonBodySetup(bD, bD.shape.volume, bD.shape.inertia, material='frictionless', pos=[0,0,0], fixed=True)
bD.state.pos = [2*lengthOfBase,0.5*heightOfBase,-0.5*lengthOfBase]
O.bodies.append(bD)
##########################################

global wallThickness
global meanSize
uf=utils.unbalancedForce()
if isnan(uf):
uf = 1.0
KE = utils.kineticEnergy()

for b in O.bodies:
if b.state.pos[1] < -5.0*meanSize and len(b.state.blockedDOFs)==0: #i.e. fixed==False
O.bodies.erase(b.id)

plot.plots={'timeStep1':('unbalancedForce'),'timeStep2':('kineticEn'),}
plot.plot() #Uncomment to view plots
O.engines=O.engines+[PyRunner(iterPeriod=20,command='still()')]
O.dt = 0.2*sqrt(0.3*O.bodies[0].state.mass/3.0e5)

def still():
if O.iter>21000 and O.iter<22000:
for b in O.bodies:
calm()

v=qt.View()

O.saveTmp()

 weijie (amandajoe) said on 2020-05-12: #8

You can change the step size of the calm () function in my code to make it appear earlier, and the same problem will occur.

 Vasileios Angelidakis (vsangelidakis) said on 2020-05-12: #9

Hi Jie,

There are several reasons causing instability:
- Your particles pass through the boundaries because the particles forming the boundaries are very thin. A value of wallThickness=0.05 worked well for me.
- frictionAngle must be given in radians, not degrees
- Some particles were erased because you had a statement to erase them inside the myAddPlotData() function

- Use larger values for "r". If your dodecahedra have an inradius d=0.107, then I would advice using r=0.053 (around half the min(d) value). Please recheck the size of your dodecahedra, because I used d=array(dd)-r in the shape class and I am not sure if you had subtracted "r" from the dd values already.
- calm() resets the velocities of all bodies; no need to put it in a loop

Last, I replaced for you the PotentialBlockVTKRecorder with export.exportPotentialBlocks, so you don't have to worry about using minAabb, maxAabb. Also, the "id" parameter is no longer used in the PotentialBlock shape class.

You can compare the updated script below with your previous script using meld [1]

All the best,
Vasileios

import numpy as np
import math
import random
import os
import errno
try:
os.mkdir('./vtk/')
except OSError as exc:
if exc.errno != errno.EEXIST:
raise
pass

#-------------------------------------------
# Engines
Kn=1e8
Ks=Kn/10

O.engines=[
ForceResetter(),
InteractionLoop(
[Ig2_PB_PB_ScGeom(calContactArea=True)],
[Ip2_FrictMat_FrictMat_KnKsPBPhys(kn_i=Kn, ks_i=Ks, Knormal=Kn, Kshear=Ks, viscousDamping=0.2)],
[Law2_SCG_KnKsPBPhys_KnKsPBLaw(allowViscousAttraction=True, traceEnergy=False)]
),
NewtonIntegrator(damping=0.0,exactAsphericalRot=True,gravity=[0,-9.81,0]),
]

#-------------------------------------------
# Dimensions
meanSize = 0.05
wallThickness = 0.5*meanSize
lengthOfBase = 0.250
heightOfBase = 0.600

#-------------------------------------------
# Make cloud
sp=pack.SpherePack()
mn,mx=Vector3(-0.5*(lengthOfBase-wallThickness),0.5*meanSize,-0.5*(lengthOfBase-wallThickness)),Vector3(0.5*(lengthOfBase-wallThickness),0.5*heightOfBase,0.5*(lengthOfBase-wallThickness))#
R=sqrt(3.0)*0.009654 #You can find exact formulae for the circumradius of dodecahedra
sp.makeCloud(mn,mx,R,0,-1,False)

#r=0.01*meanSize
aa=[0.5257311121191336, -5.836787854364519e-17, 0.5257311121191336, -5.836787854364519e-17, 0.85065080835204, -0.5257311121191336, 0.85065080835204, -0.5257311121191336, -0.85065080835204, 5.836787854364519e-17, 5.836787854364519e-17, -0.8506508083520399]
bb=[0.85065080835204, -0.5257311121191336, -0.85065080835204, 0.5257311121191336, 5.836787854364519e-17, 0.85065080835204, 5.836787854364519e-17, -0.85065080835204, 5.836787854364519e-17, 0.5257311121191336, -0.5257311121191336, -0.0]
cc=[-5.836787854364519e-17, 0.85065080835204, 5.836787854364519e-17, 0.85065080835204, -0.5257311121191336, 5.836787854364519e-17, 0.5257311121191336, 5.836787854364519e-17, 0.5257311121191336, -0.85065080835204, -0.85065080835204, -0.5257311121191336]
dd=[0.01074978698202965, 0.010749786982029651, 0.01074978698202965, 0.01074978698202965, 0.01074978698202965, 0.01074978698202965, 0.010749786982029651, 0.010749786982029651, 0.01074978698202965, 0.01074978698202965, 0.010749786982029651, 0.01074978698202965]

r=min(dd)/2 #Suggested value

for s in sp:
b=Body()
b.aspherical=True
color=Vector3(random.random(),random.random(),random.random())
b.shape=PotentialBlock(k=0.0, r=r, R=0.0, a=aa, b=bb, c=cc, d=np.array(dd)-r, isBoundary=False, color=color, wire=False, AabbMinMax=True, fixedNormal=False)
utils._commonBodySetup(b, b.shape.volume, b.shape.inertia, material='frictionless',pos=s[0], fixed=False)
b.state.pos = s[0] #s[0] stores center
b.state.ori = Quaternion((random.random(),random.random(),random.random()),random.random()) #s[2]
O.bodies.append(b)

#-------------------------------------------
# Bottom plate
color=[0,0.5,1]
r=0.15*wallThickness

bbb=Body()
wire=False
bbb.shape=PotentialBlock(k=0.0, r=r, R=0.0, a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1], d=[2.5*lengthOfBase-r,2.5*lengthOfBase-r,0.5*wallThickness-r,0.5*wallThickness-r,lengthOfBase/2.0-r,lengthOfBase/2.0-r], isBoundary=False, color=color, wire=wire, AabbMinMax=True, fixedNormal=False)
utils._commonBodySetup(bbb, bbb.shape.volume, bbb.shape.inertia, material='frictionless', pos=[0,0,0], fixed=True)
bbb.state.pos = [2*lengthOfBase,0,0]
O.bodies.append(bbb)

bA=Body()
bA.shape=PotentialBlock(k=0.0, r=r, R=0.0, a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1], d=[0.5*wallThickness-r,0.5*wallThickness-r,0.5*heightOfBase-r,0.5*heightOfBase-r,0.5*lengthOfBase-r,0.5*lengthOfBase-r], isBoundary=False, color=color, wire=False, AabbMinMax=True, fixedNormal=False)
utils._commonBodySetup(bA, bA.shape.volume, bA.shape.inertia, material='frictionless', pos=[0,0,0], fixed=True)
bA.state.pos = [0.5*lengthOfBase,0.5*heightOfBase,0]
O.bodies.append(bA)

bB=Body()
bB.shape=PotentialBlock(k=0.0, r=r, R=0.0, a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1], d=[0.5*wallThickness-r,0.5*wallThickness-r,0.5*heightOfBase-r,0.5*heightOfBase-r,0.5*lengthOfBase-r,0.5*lengthOfBase-r], isBoundary=False, color=color, wire=False, AabbMinMax=True, fixedNormal=False)
utils._commonBodySetup(bB, bB.shape.volume, bB.shape.inertia, material='frictionless', pos=[0,0,0], fixed=True)
bB.state.pos = [-0.5*lengthOfBase,0.5*heightOfBase,0]
O.bodies.append(bB)

bC=Body()
bC.shape=PotentialBlock(k=0.0, r=r, R=0.0, a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1], d=[2.5*lengthOfBase-r,2.5*lengthOfBase-r,0.5*heightOfBase-r,0.5*heightOfBase-r,0.5*wallThickness-r,0.5*wallThickness-r], isBoundary=False, color=color, wire=True, AabbMinMax=True, fixedNormal=False)
utils._commonBodySetup(bC, bC.shape.volume, bC.shape.inertia, material='frictionless', pos=[0,0,0], fixed=True)
bC.state.pos = [2*lengthOfBase,0.5*heightOfBase,0.5*lengthOfBase]
O.bodies.append(bC)

bD=Body()
bD.shape=PotentialBlock(k=0.0, r=r, R=0.0, a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1], d=[2.5*lengthOfBase-r,2.5*lengthOfBase-r,0.5*heightOfBase-r,0.5*heightOfBase-r,0.5*wallThickness-r,0.5*wallThickness-r], isBoundary=False, color=color, wire=False, AabbMinMax=True, fixedNormal=False)
utils._commonBodySetup(bD, bD.shape.volume, bD.shape.inertia, material='frictionless', pos=[0,0,0], fixed=True)
bD.state.pos = [2*lengthOfBase,0.5*heightOfBase,-0.5*lengthOfBase]
O.bodies.append(bD)

#-------------------------------------------
# Plot results
global wallThickness
global meanSize
uf=utils.unbalancedForce()
if isnan(uf):
uf = 1.0
KE = utils.kineticEnergy()

# for b in O.bodies:
# if b.state.pos[1] < -5.0*meanSize and len(b.state.blockedDOFs)==0: #i.e. fixed==False
# O.bodies.erase(b.id)

plot.plots={'timeStep1':('unbalancedForce'),'timeStep2':('kineticEn'),}
plot.plot() #Uncomment to view plots
O.engines=O.engines+[PyRunner(iterPeriod=200,command='still()',firstIterRun=4000)] #enable this engine after 4000 iterations
O.dt = 0.2*sqrt(0.3*O.bodies[0].state.mass/Kn)*100

#-------------------------------------------
# Invoke calm()
def still():
KE = utils.kineticEnergy()
if KE<0.01:
calm()

#-------------------------------------------
# exportPotentialBlocks
vtkExporter_PotentialBlocks = export.VTKExporter('vtk/pb')

def vtkExport():
vtkExporter_PotentialBlocks.exportPotentialBlocks(what=dict(n='b.id'))

O.engines=O.engines+[PyRunner(iterPeriod=500,command='vtkExport()')]

v=qt.View()

O.saveTmp()

 Vasileios Angelidakis (vsangelidakis) said on 2020-05-12: #10

A correction: I used meanSize = 0.05. wallThickness is half that.

Cheers,
Vasileios

 weijie (amandajoe) said on 2020-05-18: #11

Thanks Vasileios Angelidakis, that solved my question.

 Bruno Chareyre (bruno-chareyre) said on 2020-05-18: #12

A short comment for future reference: in my experience tricks like "calm()" are needed only when there is another problem in parameters or boundary conditions (e.g. too dynamic).
When the parameters are not doomed calm() only makes the whole thing slower (which is logical: if a system is unbalanced, freezing it will not make it converge faster to equilibrium). It is thus more efficient to track the initial problem when possible.
Bruno