How to relate to each other the stresses calculated respectively from force on the wall and contact force

Dear All,

I added left, right, front and back walls to the script of rot_compression.py in the Tutorial Version 2.3,
and calculated the stresses via two methods ( one was the same as in the tutorial using forces on the wall and the position of wall, the other was via interparticle contact force and branchvector).
However, I found that the stresses from the two methods differ from each other.
So, What's the relation between the stresses from the two methods.

Related codes are as following

1) modified rot_compression.py
#scripts from rot_compress.py in Tutorial version 2.3
#import the appropriate ESyS-Particle modules:
from esys.lsm import *
from esys.lsm.util import *
from esys.lsm.geometry import *
from WallLoader import WallLoaderRunnable
from math import *

print pi

#create a simulation container object:
# N.B. there must be at least two sub-divisions
# in the X-direction for periodic boundaries
sim = LsmMpi (numWorkerProcesses=2, mpiDimList=[1,2,0])
sim.initNeighbourSearch (
 particleType = "RotSphere",
 gridSpacing = 5.0000,
 verletDist = 0.08000
)

#specify the number of timesteps and timestep increment:
sim.setNumTimeSteps(250000)
sim.setTimeStepSize(1.0000e-06)

#specify the spatial domain and direction of periodic boundaries:
domain = BoundingBox ( Vec3 (-20,-20,-20), Vec3 (20,20,20) )
sim.setSpatialDomain (domain)

#construct a prism of spherical particles:
geoRandomBlock = RandomBoxPacker (
 minRadius = 0.400,
 maxRadius = 2.000,
 cubicPackRadius = 2.2000,
 maxInsertFails = 5000,
 bBox = BoundingBox(
 Vec3(-5.0, 0.0,-5.0),
 Vec3(5.0, 20.0, 5.0)
 ),
 circDimList = [False, False, False],
 tolerance = 1.0e-5
)
geoRandomBlock.generate()
geoRandomBlock_particles = geoRandomBlock.getSimpleSphereCollection()

sim.createParticles(geoRandomBlock_particles)

#bond particles together with bondTag = 1:
sim.createConnections(
 ConnectionFinder(
  maxDist = 0.005,
  bondTag = 1,
  pList = geoRandomBlock_particles
 )
)

#create a wall on the left of the model:
sim.createWall (
 name = "left_wall",
 posn = Vec3(-5.0000, 0.0000, -5.0000),
 normal = Vec3(1.0000, 0.0000, 0.0000)
)

#create a wall at the bottom of the model:
sim.createWall (
 name = "bottom_wall",
 posn = Vec3(-5.0000, 0.0000, -5.0000),
 normal = Vec3(0.0000, 1.0000, 0.0000)
)

#create a wall at the back of the model:
sim.createWall (
 name = "back_wall",
 posn = Vec3(-5.0000, 0.0000, -5.0000),
 normal = Vec3(0.0000, 0.0000, 1.0000)
)

#create a wall at the right of the model:
sim.createWall (
 name = "right_wall",
 posn = Vec3(20.0000, 20.0000, 20.0000),
 normal = Vec3(-1.0000, 0.0000, 0.0000)
)

#create a wall at the top of the model:
sim.createWall (
 name = "top_wall",
 posn = Vec3(20.0000, 20.0000, 20.0000),
 normal = Vec3(0.0000, -1.0000, 0.0000)
)

#create a wall at the top of the model:
sim.createWall (
 name = "front_wall",
 posn = Vec3(20.0000, 20.0000, 20.0000),
 normal = Vec3(0.0000, 0.0000, -1.0000)
)

#create rotational elastic-brittle bonds between particles:
pp_bonds = sim.createInteractionGroup (
 BrittleBeamPrms(
  name="pp_bonds",
  youngsModulus=100000.0,
  poissonsRatio=0.25,
  cohesion=100.0, #100.0
  tanAngle=1.0, #1.0
  tag=1
 )
)

#initialise frictional interactions for unbonded particles:
sim.createInteractionGroup (
 FrictionPrms(
  name="friction",
  youngsModulus=100000.0,
  poissonsRatio=0.25,
  dynamicMu=0.4, # 0.4
  staticMu=0.6 # 0.6
 )
)

#create an exclusion between bonded and frictional interactions:
sim.createExclusion (
 interactionName1 = "pp_bonds",
 interactionName2 = "friction"
)

#specify elastic repulsion from the bottom wall:
sim.createInteractionGroup (
 NRotElasticWallPrms (
  name = "bottom_wall_repel",
  wallName = "bottom_wall",
  normalK = 100000.0
 )
)

#specify elastic repulsion from the top wall:
sim.createInteractionGroup (
 NRotElasticWallPrms (
  name = "top_wall_repel",
  wallName = "top_wall",
  normalK = 100000.0
 )
)

#specify elastic repulsion from the left wall:
sim.createInteractionGroup (
 NRotElasticWallPrms (
  name = "left_wall_repel",
  wallName = "left_wall",
  normalK = 100000.0
 )
)

#specify elastic repulsion from the right wall:
sim.createInteractionGroup (
 NRotElasticWallPrms (
  name = "right_wall_repel",
  wallName = "right_wall",
  normalK = 100000.0
 )
)

#specify elastic repulsion from the back wall:
sim.createInteractionGroup (
 NRotElasticWallPrms (
  name = "back_wall_repel",
  wallName = "back_wall",
  normalK = 100000.0
 )
)

#specify elastic repulsion from the top wall:
sim.createInteractionGroup (
 NRotElasticWallPrms (
  name = "back_wall_repel",
  wallName = "back_wall",
  normalK = 100000.0
 )
)

#add translational viscous damping:
sim.createInteractionGroup (
 LinDampingPrms(
  name="damping1",
  viscosity=0.002,
  maxIterations=50
 )
)

#add rotational viscous damping:
sim.createInteractionGroup (
 RotDampingPrms(
  name="damping2",
  viscosity=0.002,
  maxIterations=50
 )
)

#add a wall loader to move the top wall:
wall_loader1 = WallLoaderRunnable(
 LsmMpi = sim,
 wallName = "top_wall",
 vPlate = Vec3 (0.0, -0.125, 0.0),
 startTime = 0,
 rampTime = 50000
)
sim.addPreTimeStepRunnable (wall_loader1)

#add a wall loader to move the bottom wall:
wall_loader2 = WallLoaderRunnable(
 LsmMpi = sim,
 wallName = "bottom_wall",
 vPlate = Vec3 (0.0, 0.125, 0.0),
 startTime = 0,
 rampTime = 50000
)
sim.addPreTimeStepRunnable (wall_loader2)

#create a FieldSaver to store number of bonds:
sim.createFieldSaver (
 InteractionScalarFieldSaverPrms(
  interactionName="pp_bonds",
  fieldName="count",
  fileName="nbonds.dat",
  fileFormat="SUM",
  beginTimeStep=0,
  endTimeStep=250000,
  timeStepIncr=1000
 )
)

#create a FieldSaver to store the total kinetic energy of the particles:
sim.createFieldSaver (
 ParticleScalarFieldSaverPrms(
  fieldName="e_kin",
  fileName="ekin.dat",
  fileFormat="SUM",
  beginTimeStep=0,
  endTimeStep=250000,
  timeStepIncr=1000
 )
)

#create a FieldSaver to store potential energy stored in bonds:
sim.createFieldSaver (
 InteractionScalarFieldSaverPrms(
  interactionName="pp_bonds",
  fieldName="potential_energy",
  fileName="epot.dat",
  fileFormat="SUM",
  beginTimeStep=0,
  endTimeStep=250000,
  timeStepIncr=1000
 )
)

#create a FieldSaver to wall positions:
posn_saver = WallVectorFieldSaverPrms(
 wallName=["bottom_wall", "top_wall","left_wall","right_wall","back_wall", "front_wall"],
 fieldName="Position",
 fileName="out_Position.dat",
 fileFormat="RAW_SERIES",
 beginTimeStep=0,
 endTimeStep=250000,
 timeStepIncr=1000
)
sim.createFieldSaver(posn_saver)

#create a FieldSaver to wall forces:
force_saver = WallVectorFieldSaverPrms(
 wallName=["bottom_wall", "top_wall","left_wall","right_wall","back_wall", "front_wall"],
 fieldName="Force",
 fileName="out_wallForce.dat",
 fileFormat="RAW_SERIES",
 beginTimeStep=0,
 endTimeStep=250000,
 timeStepIncr=1000
)
sim.createFieldSaver(force_saver)

#create a FieldSaver to contact force and particle positions
bond_contact_force_saver = InteractionVectorFieldSaverPrms (
   interactionName = "pp_bonds",
   fieldName = "force",
   fileName = "bond_contactForces",
   fileFormat = "RAW_WITH_POS_ID",
   beginTimeStep = 0,
   endTimeStep = 250000,
   timeStepIncr = 1000
)
sim.createFieldSaver(bond_contact_force_saver)

#create a FieldSaver to contact force and particle positions
friction_contact_force_saver = InteractionVectorFieldSaverPrms (
   interactionName = "friction",
   fieldName = "force",
   fileName = "friction_contactForces",
   fileFormat = "RAW_WITH_POS_ID",
   beginTimeStep = 0,
   endTimeStep = 250000,
   timeStepIncr = 1000
)
sim.createFieldSaver(friction_contact_force_saver)

#execute the simulation:
sim.run()

2) code for calculating the stress via forces on wall and the position of wall
from math import *

posnfile = open("out_Position.dat","r")
posn = posnfile.readlines()
posnfile.close()

forcefile = open("out_wallForce.dat","r")
force = forcefile.readlines()
forcefile.close()

for i in range (len(posn)):
 Y_bottom = float(posn[i].split()[1])
 Y_top = float(posn[i].split()[4])
 F_bottom = float(force[i].split()[1])
 F_top = float(force[i].split()[4])

 stress = (F_bottom - F_top)/200.0
 strain = 1.0 - (Y_top - Y_bottom)/20.0
 print i*1000,strain,stress
3) code for calculating stress via interparticle contact forces and branchvector
from math import *

posnfile = open("out_Position.dat","r")
posn = posnfile.readlines()
posnfile.close()

for j in range(251):

 #calculating volume
 X_left = float(posn[j].split()[6])
 X_right = float(posn[j].split()[9])
 Y_bottom = float(posn[j].split()[1])
 Y_top = float(posn[j].split()[4])
 Z_back = float(posn[j].split()[14])
 Z_front = float(posn[j].split()[17])
 volume = (Z_front-Z_back)*(Y_top-Y_bottom)*(X_right-X_left)

 #set zero
 stress = [[0,0,0],[0,0,0],[0,0,0]]
 force = [0,0,0]
 branchvector = [0,0,0]

 #calculate the stress from bond_contactForces
 bondcontactforcefilename = "bond_contactForces.{0:0d}.dat".format(j)
 bondcontactforcefile = open(bondcontactforcefilename,"r")
 bondcontactforce = bondcontactforcefile.readlines()
 bondcontactforcefile.close()

 for i in range (len(bondcontactforce)):

  force[0] = float(bondcontactforce[i].split()[11])
  force[1] = float(bondcontactforce[i].split()[12])
  force[2] = float(bondcontactforce[i].split()[13])

  branchvector[0] = float(bondcontactforce[i].split()[5])-float(bondcontactforce[i].split()[2])
  branchvector[1] = float(bondcontactforce[i].split()[6])-float(bondcontactforce[i].split()[3])
  branchvector[2] = float(bondcontactforce[i].split()[7])-float(bondcontactforce[i].split()[4])

  stress[0][0] += force[0]*branchvector[0]/volume
  stress[0][1] += force[0]*branchvector[1]/volume
  stress[0][2] += force[0]*branchvector[2]/volume

  stress[1][0] += force[1]*branchvector[0]/volume
  stress[1][1] += force[1]*branchvector[1]/volume
  stress[1][2] += force[1]*branchvector[2]/volume

  stress[2][0] += force[2]*branchvector[0]/volume
  stress[2][1] += force[2]*branchvector[1]/volume
  stress[2][2] += force[2]*branchvector[2]/volume

 #calculate the stress from friction_contactForces
 frictioncontactforcefilename = "friction_contactForces.{0:0d}.dat".format(j)
 frictioncontactforcefile = open(frictioncontactforcefilename,"r")
 frictioncontactforce = frictioncontactforcefile.readlines()
 frictioncontactforcefile.close()

 for i in range (len(frictioncontactforce)):

  force[0] = float(frictioncontactforce[i].split()[11])
  force[1] = float(frictioncontactforce[i].split()[12])
  force[2] = float(frictioncontactforce[i].split()[13])

  branchvector[0] = float(frictioncontactforce[i].split()[5])-float(frictioncontactforce[i].split()[2])
  branchvector[1] = float(frictioncontactforce[i].split()[6])-float(frictioncontactforce[i].split()[3])
  branchvector[2] = float(frictioncontactforce[i].split()[7])-float(frictioncontactforce[i].split()[4])

  stress[0][0] += force[0]*branchvector[0]
  stress[0][1] += force[0]*branchvector[1]
  stress[0][2] += force[0]*branchvector[2]

  stress[1][0] += force[1]*branchvector[0]
  stress[1][1] += force[1]*branchvector[1]
  stress[1][2] += force[1]*branchvector[2]

  stress[2][0] += force[2]*branchvector[0]
  stress[2][1] += force[2]*branchvector[1]
  stress[2][2] += force[2]*branchvector[2]

 #output stress in TimeStep j*timeStepIncr

 print j*1000,volume,stress[0][0],stress[0][1],stress[0][2],stress[1][0],stress[1][1],stress[1][2],stress[2][0],stress[2][1],stress[2][2]

Best regards,
Jiadun