Energy tracking of a two-sphere contact

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()