plastic dissipation calculation does't match with the result from yade

Hello,

I'm trying to compute the external work, elastic energy and plastic dissipation of a specimen along the biaxial loading test. The final goal is to analyze energy dissipation in specific area, so I want to compare my results (energies in the whole speicmens) with yade at first. The external work and elastic energy from my code are closed to the results by yade. However, the plastic dissipation is about 0.6~0.7 times smaller than the result in yade. In fact, the code for plastic dissipation part is copied from issue #1899 (closed) which I have referred to and modified to avoid repeatly write/read .txt files.

The part of code related to energy calculation is shown as below. It could be have some mistakes and I need someone help me to fix it.

Thanks in advance!

# calculate energy at each time step
def calEnergy():
    global strain0
    global externalWork, elasticStrainEnergy, plasticStrainEnergy, plasticStrainEnergy1
    global current_ExtWork, current_ElasticEnergy, current_PlasticEnergy

    currentElasticEnergy = 0
    currentPlasticEnergy = 0

    t0_Eel = elasticStrainEnergy

    # incremental external work
    disx = triax1.width
    disy = triax1.height
    disz = triax1.depth
    vol = disx * disy * disz

    s11 = -triax1.stress(triax1.wall_right_id)[0]
    s22 = -triax1.stress(triax1.wall_top_id)[1]
    s33 = -triax1.stress(triax1.wall_front_id)[2]

    strain = -triax1.strain
    de11 = strain[0]-strain0[0]
    de22 = strain[1]-strain0[1]
    de33 = strain[2]-strain0[2]

    dW_ext = vol * (s11 * de11 + s22 * de22 + s33 * de33)  # external power
    currentWork = dW_ext  # * O.dt
    externalWork += currentWork
    current_ExtWork = currentWork

    strain0 = strain

    # elastic and plastic energy
    Ft_dict = {}
    for k in O.interactions:
        if k.isReal and isinstance(O.bodies[k.id1].shape, Sphere) and isinstance(O.bodies[k.id2].shape, Sphere):
            phys = k.phys
            normalForce = phys.normalForce.norm()
            shearForce = phys.shearForce.norm()
            shearF = phys.shearForce
            kn = phys.kn
            ks = phys.ks

            # elastic_energy = 0.5 * physics.kn * k.geom.penetrationDepth ** 2 + 0.5 * physics.ks * physics.shearForce.norm() ** 2
            elastic_energy = 0.5 * (normalForce ** 2 / kn + shearForce ** 2 / ks)
            currentElasticEnergy += elastic_energy

            Ft_dict[tuple((k.id1, k.id2))] = shearF
    elasticStrainEnergy = currentElasticEnergy
    current_ElasticEnergy = currentElasticEnergy

    O.step()

    for k in O.interactions:
        if k.isReal and isinstance(O.bodies[k.id1].shape, Sphere) and isinstance(O.bodies[k.id2].shape, Sphere):
            ids = tuple((k.id1, k.id2))
            Ftt = k.phys.shearForce
            du = k.geom.shearInc
            kt = k.phys.ks
            Fn = k.phys.normalForce
            maxFs = Fn.norm() * k.phys.tangensOfFrictionAngle
            Ft = Vector3(0, 0, 0)
            if ids in Ft_dict.keys():
                Ft = Ft_dict[ids]
            Ftc = Ft - kt * du
            if Ftc.norm() > maxFs:
                # Ftt=Ftt/Ftt.norm()*maxFs
                due = -(Ftt - Ft) / kt
                # e+=Ftt.norm()*(du.norm()-due.norm())
                currentPlasticEnergy += Ftt.norm() * (du - due).norm()
    plasticStrainEnergy += currentPlasticEnergy
    current_PlasticEnergy = currentPlasticEnergy

and call this function by:

while record_check == 0:
    if steps < 100:
        O.step()
        steps += 1
        calEnergy()
    else:
        print(externalWork, elasticStrainEnergy - elastic_E0, plasticStrainEnergy, plasticStrainEnergy1)
        e22 = -triax1.strain[1] - strain00[1]
        if (e22 - ini_e22b) >= 1e-4:
            with open('energy.txt', 'a') as f:
                f.write('%d %.5f %.5f %.5f %.5f %.5f %.5f %.5f %.5f\n' % (O.iter, -triax1.strain[1] - strain00[1], current_ExtWork, externalWork, elasticStrainEnergy - elastic_E0, plasticStrainEnergy, plasticStrainEnergy1, -triax1.externalWork - external_W0, elasticStrainEnergy + plasticStrainEnergy - elastic_E0))
            addPlotData()
            ini_e22b = e22
            if (e22 - ini_e22a) >= 0.0005:
                O.save('./states/save' + '{:.4f}'.format(e22) + 'yade.bz2')
                ini_e22a = e22
        steps = 0