New question #446267 on Yade:
https://answers.launchpad.net/yade/+question/446267
I am trying to follow the internal energy of a simple two-sphere contact and
compare it with the external work. However for the case of a combination
between normal and shear stress, there seems to be a discrepancy between the
internal energy and external work. One sphere is fixed, while the other one is
displacement-controlled and I'm using a cohesive-frictional contact law, which
allows me to model a cohesive bond between the spheres. From what I am able to
tell, the discrepancy between the internal energy and external work in my
simulation seems to origin from two separate issues:
1.) An overestimation of the bond braking energy: Since following the
bond-breaking energy is not possible, I wrote my own code in order to monitor
it. The code is based on simply following the elastic energy difference before
and after the bond breakage. It seems to work fine in most scenarios, but in
this particular kinematics of two spheres, it seems to overestimate the bond
breaking energy. Or perhaps the error is somewhere else? What I find
particularly strange is the shear dissipation energy jump at bond breaking.
Since the simulation is displacement controlled, that is not what I would have
expected.
2.) A slowly growing discrepancy between the internal energy and external work:
This issue is independent of the cohesion - it appears even if the cohesive
bond in not created (if the lines 133-138 are commented). It also seems to
appear with other contact laws - I have reproduced it with a combination of
FrictMat and CundallStrack contact law. I know, that the imposed kinematics are
somewhat unnatural, due to being displacement controlled and produce big shear
displacements, but I suppose i should still be able to follow the contact
energy?
I deliberately chose a very simple scenario in order to demonstrate the issue
I'm experiencing. If there is an error in my code, I will be very glad if you
point it out to me. I'm attaching my code to this message:
''' Two-sphere contact '''
import numpy as np
import sys, os, os.path
import time
import operator
import timeit
import matplotlib.pyplot as plt
import yade
from yade import plot, export
##########Parameters##########
density = 9700
damping = 0
velocity = -1e-2
frictionAngle = 0.2
cohesion = 1e6
young = 1e7 *10
####Material######
Cohesion_material = yade.CohFrictMat(young=young,poisson=0.3,density=density,
frictionAngle=frictionAngle,
normalCohesion=cohesion,
shearCohesion=cohesion,
momentRotationLaw=False,
etaRoll=-1,alphaKr=0,alphaKtw=0,label='cohesive'
)
yade.O.materials.append(Cohesion_material)
data, energys, bond_energy = [], [], []
broken_bonds_dissipation = 0
external_work = 0
up,Fp = [0,0,0],[0,0,0]
def
plot_energy(previous_normEnergy,previous_shearEnergy,previous_broken_bonds_list,broken_bonds_dissipation,external_work,up,Fp):
#displacement and forces on the top sphere
u
=(yade.O.bodies[id_top_clump].state.pos-yade.O.bodies[id_top_clump].state.refPos)
F = yade.O.forces.f(id_top_clump)
du = u-up
#energy tracking
try:
Ek = O.energy['kinetic']
except KeyError:
Ek = 0
if damping == 0:
Edamp = 0
else:
Edamp = yade.O.energy['nonviscDamp']
Een = myLaw2s[0].normElastEnergy()
Ees = myLaw2s[0].shearElastEnergy()
Eplast = myLaw2s[0].plasticDissipation()
Etot = yade.O.energy.total()
#calculate the energy of a particular cohesive bond
normEnergy, shearEnergy = {},{}
for i in cohesive_interactions:
if i.phys:
normEnergy[i] = 0.5*(i.phys.normalForce.squaredNorm()/i.phys.kn)
shearEnergy[i] = 0.5*(i.phys.shearForce.squaredNorm()/i.phys.ks)
else:
normEnergy[i] = 0
shearEnergy[i] = 0
broken_bonds = 0
broken_bonds_list = {}
for i in cohesive_interactions:
try:
if i.phys.cohesionBroken == False:
broken_bonds_list[i] = 0
continue
else:
broken_bonds = broken_bonds + 1
broken_bonds_list[i] = 1
except AttributeError:
broken_bonds = broken_bonds + 1
broken_bonds_list[i] = 1
for i in cohesive_interactions:
if broken_bonds_list[i] == 1 and previous_broken_bonds_list[i] == 0:
broken_bonds_dissipation = broken_bonds_dissipation +
(previous_normEnergy[i]-normEnergy[i]) +
(previous_shearEnergy[i]-shearEnergy[i])
previous_normEnergy = normEnergy.copy()
previous_shearEnergy = shearEnergy.copy()
previous_broken_bonds_list = broken_bonds_list.copy()
#data to plot
external_work = external_work + Fp[0]*du[0]
internal_energy = Ek + Een + Ees + Eplast + Edamp + broken_bonds_dissipation
ux,uy,uz = u
plot.addData(ux=ux,internal_energy=internal_energy,external_work=-external_work,Eplast=Eplast)
#export energy parameters into the next iteration
up=u
Fp=F
return
previous_normEnergy,previous_shearEnergy,previous_broken_bonds_list,broken_bonds_dissipation,external_work,up,Fp
#Defining interaction geometry, physics and law
myIgeoms =
[yade.Ig2_Sphere_Sphere_ScGeom6D(),yade.Ig2_Facet_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom6D()]
myIphyss = [yade.Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=False)]
myLaw2s =
[yade.Law2_ScGeom6D_CohFrictPhys_CohesionMoment(useIncrementalForm=True,
always_use_moment_law=False)]
yade.O.engines=[
yade.ForceResetter(),
yade.InsertionSortCollider([yade.Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Box_Aabb()]),
yade.InteractionLoop(myIgeoms,myIphyss,myLaw2s),
yade.NewtonIntegrator(damping=damping),
yade.PyRunner(command='previous_normEnergy,previous_shearEnergy,previous_broken_bonds_list,broken_bonds_dissipation,external_work,up,Fp=plot_energy(previous_normEnergy,previous_shearEnergy,previous_broken_bonds_list,broken_bonds_dissipation,external_work,up,Fp)',iterPeriod=1),
]
#sphere construction
b1 = yade.utils.sphere(center =
(-0.5,0,0),radius=1,material='cohesive',wire=False)
b2 = yade.utils.sphere(center =
(0,0,1.95),radius=1,material='cohesive',wire=False)
id_top_clump = yade.O.bodies.append(b1)
id_bottom_clump = yade.O.bodies.append(b2)
#cohesive interactions
intr = yade.utils.createInteraction(id_top_clump,id_bottom_clump)
intr.phys.unp =
-(yade.O.bodies[b1.id].state.pos-yade.O.bodies[b2.id].state.pos).norm()+\
yade.O.bodies[b1.id].shape.radius +
yade.O.bodies[b2.id].shape.radius
intr.phys.cohesionBroken = False
intr.phys.normalAdhesion = 1000000
intr.phys.shearAdhesion = 1000000
cohesive_interactions = []
for i in yade.O.interactions:
cohesive_interactions.append(i)
previous_normEnergy, previous_shearEnergy, previous_broken_bonds_list = {}, {},
{}
for i in cohesive_interactions:
previous_shearEnergy[i] = 0
previous_normEnergy[i] = 0
previous_broken_bonds_list[i] = 0
#kinematic constraints
yade.O.bodies[id_bottom_clump].state.blockedDOFs = 'xyzXYZ'
yade.O.bodies[id_top_clump].state.blockedDOFs = 'xyzXYZ'
yade.O.bodies[id_top_clump].state.vel = [-velocity,0,0]
yade.O.trackEnergy = True
myLaw2s[0].traceEnergy = True
yade.utils.setRefSe3() #Set reference position and orientation as the current
ones.
plot.plots={'ux':('internal_energy','external_work','Eplast')}
plot.plot()
yade.O.dt=0.05*PWaveTimeStep()
O.run()
--
You received this question notification because your team yade-users is
an answer contact for Yade.
_______________________________________________
Mailing list: https://launchpad.net/~yade-users
Post to : [email protected]
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
