Wrong movement of wall under flow triaxial compression

Hello,
I'm simulating triaxial compression under Z-direction flow pressure.The confining pressure I set is 20MPa and the osmotic pressure is 5MPa in Z-direction(script in the end).
However,compared with the case without osmotic pressure, I found that when the simulation termination condition is reached, the wall in X direction moves about 25mm in the negative direction, the wall in Y direction moves about 15mm in the positive direction, and the wall in Z direction does not move significantly.
My guess is whether this is related to the deviation of the wall caused by the excessive osmotic pressure I set?
Thanks for help.

--------------------------script--------------------
from yade import pack, ymport, plot, utils, export, timing
Wy=-20e6
Sy=-5e6
rate=-0.1
damp=0.4
stabilityThreshold=0.001
key='_triax_base_'
young=3000e9
name='JCFPM_triax'
compFricDegree=30
poisson=0.4

mn,mx=Vector3(0,0,0),Vector3(0.05,0.05,0.05)
O.materials.append(JCFpmMat(type=1,density=2640,young=young,poisson=poisson,tensileStrength=160e6,cohesion=1600e6,frictionAngle=radians(10),label='sphere'))
O.materials.append(JCFpmMat(type=1,frictionAngle=0,density=0,label='wall'))

walls=aabbWalls([mn,mx],thickness=0,material='wall')
wallIds=O.bodies.append(walls)
O.bodies.append(ymport.text('packing-JCFPM.spheres',scale=1,shift=Vector3(0,0,0),material='sphere',color=(0.6,0.5,0.15)))

#(The model here is a cube generated by RandomDensePack,composed of 5000 particles with a side length of 0.05m )

triax=TriaxialStressController(
 maxMultiplier=1.+4e8/young,
 finalMaxMultiplier=1.+4e7/young,
 thickness = 0,
 stressMask = 7,
 internalCompaction=True,
)
newton=NewtonIntegrator(damping=damp)
O.timingEnabled=True

def recorder():
 yade.plot.addData(
 i=O.iter,
 e11=-triax.strain[0],#e22=-triax.strain[1],
 e33=-triax.strain[2],
 s11=-triax.stress(triax.wall_right_id)[0],#0+边界id为right
 s33=-triax.stress(triax.wall_front_id)[2],#2+边界id为front
 numberTc=interactionLaw.nbTensCracks,
 numberSc=interactionLaw.nbShearCracks,
 unb=unbalancedForce(),
 flux1=flow.getBoundaryFlux(4),
 flux2=flow.getBoundaryFlux(5)
)
 plot.saveDataTxt('triax_JCFPM_HM_Wy=20.0Mpa,Sy=5.0MPa')

def stop_condition():
 extremum=max(abs(s) for s in plot.data['s33'])
 s=abs(plot.data['s33'][-1])
 e=abs(plot.data['e33'][-1])
 if e < 0.001:
  return
 if abs(s)/abs(extremum) < 0.9 :
  print('Max stress and strain:',extremum,e)
  presentcohesive_count = 0
  for i in O.interactions:
          if hasattr(i.phys, 'isCohesive'):
               if i.phys.isCohesive == 1:
                   presentcohesive_count+=1
  print('the number of cohesive bond now is:',presentcohesive_count)
  O.pause()

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=1.5,label='is2aabb'),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=1.5,label='ss2sc'),Ig2_Box_Sphere_ScGeom()],
  [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys()],
  [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(recordCracks=True,Key='triax_JCFPM_HM',label='interactionLaw')]
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.7),
 triax,
 FlowEngine(dead=1,label='flow'),
 PyRunner(iterPeriod=int(1000),initRun=True,command='recorder()',label='data',dead=0),
 PyRunner(iterPeriod=1000,command='stop_condition()',dead=0),
 VTKRecorder(iterPeriod=500,initRun=True,fileName='triax/JFFPM-',recorders=['spheres','jcfpm','cracks'],Key='triax',label='vtk',dead=1),
 newton,
]

O.step()
ss2sc.interactionDetectionFactor=-1
is2aabb.aabbEnlargeFactor=-1
cohesiveCount = 0
for i in O.interactions:
 if hasattr(i.phys, 'isCohesive'):
              if i.phys.isCohesive == True:
                     cohesiveCount+=1
print('the origin total number of cohesive bond is:',cohesiveCount)
triax.goal1=triax.goal2=triax.goal3=Wy

while 1:
 O.run(1000,1)
 unb=unbalancedForce()
 print('unbalanced force:',unb,'mean stres:',triax.meanStress)
 if unb<stabilityThreshold and abs(Wy-triax.meanStress)/abs(Wy)<0.001:
  break

yade.timing.reset()
flow.dead=0
flow.debug=False
flow.fluidBulkModulus=2.2e9
flow.permeabilityFactor=5
flow.decoupleForces = False
flow.meshUpdateInterval=1000
flow.useSolver=4
flow.viscosity=0.001
flow.bndCondIsPressure=[0,0,0,0,1,1]
flow.bndCondValue=[0,0,0,0,0,Sy]

triax.internalCompaction=False
flow.emulateAction()
triax.stressMask=3
triax.goal1=Wy
triax.goal2=Wy
triax.goal3=rate

plot.plots={'e33':(('s33','g'),None,('unb','b')),'i':(('numberTc','b'),('numberSc','r'),None,('s33','y'))}
plot.plot(subPlots=False)
O.run()
-------------------------------------------------