Questions about the load plate test for sandy soil

Hello everyone,
    I am trying to calibrate the contact parameters of sand with the plate load tests in article [1].
    However, due to the large size of the subsequent model I want to build, the number of particles may be too high if using the real sand particle size (maximum around 2mm), so I choose to enlarge the particle size to 3~15mm.
    After that, I have built the model according to the whole size in [1], where the script after sample making is shown in（I am sorry, I have no idea how to upload the dat file after making the sample）.The script of the plate load tests model is as follows：

**************************************************************************************************************************
from yade import pack,qt,plot,utils,polyhedra_utils,ymport,export,pack,timing
from yade import *
import numpy
from pprint import pprint
import random
from random import uniform
from random import randint
from math import *

steel = PolyhedraMat() #the material of wall the parameter is identical with the model in chaning.
steel.density = 7800 #kg/m^3
steel.young = 8E8 #inital steel was ,try
steel.poisson = 0.21
steel.frictionAngle = 0.3 #rad
O.materials.append(steel)

steelp = PolyhedraMat() #the material of plate
steel.density = 2000 #kg/m^3 in order to ignore the gravity of plate
steel.young = 20E9 #inital steel was ,try
steel.poisson = 0.21
steel.frictionAngle = 0.3 #rad
O.materials.append(steel)

spheresand = CohFrictMat()
spheresand.isCohesive = False
spheresand.alphaKr = 1.7 # 0.7
spheresand.alphaKtw = 0
spheresand.momentRotationLaw = True
spheresand.density = 2600
spheresand.young = E7 #in October 3th
spheresand.possion = 0.33
spheresand.frictionAngle = 0.977 #rad 0.77
O.materials.append(spheresand)

'''spheresand = FrictMat()
spheresand.density = 2600
spheresand.young = 8E7 #in10.3
spheresand.possion = 0.33
spheresand.frictionAngle = 0.68 #rad
O.materials.append(spheresand)'''

global a
a=1
#length unit: m; using stander unit

O.bodies.append(utils.wall((-0.19,0,0),0,sense=1,material=steel))
O.bodies.append(utils.wall((0.19,0,0),0,sense=-1,material=steel))
O.bodies.append(utils.wall((0,-0.24,0),1,sense=1,material=steel))
O.bodies.append(utils.wall((0,0.24,0),1,sense=-1,material=steel))
O.bodies.append(utils.wall((0,0,0),2,sense=1,material=steel))

'''sp=pack.SpherePack()
sp.makeCloud(minCorner=(-0.19,-0.24,0),maxCorner=(0.19,0.24,0.1),porosity=0.38,psdSizes=[0.0025,0.0033,0.0039,0.0052,0.0063,0.0085,0.0127,0.0143,0.0169,0.0184],psdCumm=[0,0.066,0.141,0.272,0.359,0.559,0.838,0.924,0.986,1])
sp.toSimulation(material=spheresand)
'''
O.bodies.append(ymport.textExt('/tmp/nd0.dat',format='x_y_z_r',material=spheresand))
#O.bodies.append(ymport.textExt('/tmp/newsx3.dat',format='x_y_z_r',shift=Vector3(0,0,0.0496),material=spheresand))

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Polyhedra_Aabb(),Bo1_Wall_Aabb(),Bo1_Facet_Aabb(),Bo1_Sphere_Aabb()]),
 InteractionLoop(
  [Ig2_Wall_Polyhedra_PolyhedraGeom(),Ig2_Polyhedra_Polyhedra_PolyhedraGeom(),Ig2_Facet_Polyhedra_PolyhedraGeom(),Ig2_Sphere_Sphere_ScGeom(),Ig2_Sphere_Polyhedra_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
  [Ip2_PolyhedraMat_PolyhedraMat_PolyhedraPhys(),Ip2_FrictMat_PolyhedraMat_FrictPhys(),Ip2_FrictMat_FrictMat_FrictPhys(),Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(),],
  [Law2_PolyhedraGeom_PolyhedraPhys_Volumetric(),Law2_ScGeom_FrictPhys_CundallStrack(),Law2_ScGeom6D_CohFrictPhys_CohesionMoment()], # contact law -- apply forces #Bo1_Cylinder_Aabb()
 ),
 NewtonIntegrator(damping=0.3,gravity=(0,0,-9.8)),
 PyRunner(command='cycle()',iterPeriod=300,label='step'),
 PyRunner(command='Compact()',iterPeriod=1,label='compact'),
 PyRunner(command='save()',iterPeriod=100,label='save')
]

def cycle():
 for b in O.bodies:
  if isinstance(b.shape,Sphere):
   if b.state.pos[0] < -0.2 or b.state.pos[0] > 0.2:
    O.bodies.erase(b.id)
   elif b.state.pos[1] < -0.26 or b.state.pos[1] > 0.26:
    O.bodies.erase(b.id)
   elif b.state.pos[2] < -0.2 or b.state.pos[2] > 0.079:
    O.bodies.erase(b.id)
   else:
    pass
 #for b in O.bodies:
  #if isinstance(b.shape,Sphere):
   #b.state.vel=(0,0,0)
   #b.state.angVel=(0,0,0)
   #b.state.angMom=(0,0,0)

def Compact():
 global a
 if a == 1:
  ldpltheight=max([b.state.pos[2] for b in O.bodies if isinstance(b.shape,Sphere)])
  print ("pos=%.5f"%(ldpltheight))
  if O.iter < 20: return
  ldpltheight=0.0745
  ldplt=polyhedra_utils.polyhedra(steel,v=((0.045,0.05,ldpltheight),(0.045,-0.05,ldpltheight),(-0.045,-0.05,ldpltheight),(-0.045,0.05,ldpltheight),(0.045,0.05,ldpltheight+0.02),(0.045,-0.05,ldpltheight+0.02),(-0.045,-0.05,ldpltheight+0.02),(-0.045,0.05,ldpltheight+0.02)),fixed=False, color=(0.75,0.65,0.65))
  O.bodies.append(ldplt)
  a=a+1
 elif a == 2:
  Lastnum=O.bodies[-1].id
  O.bodies[Lastnum].state.blockedDOFs='xyXYZ'
  O.forces.setPermF(Lastnum,(0,0,-450))
  a=a+1
 elif a == 3:
  Lastnum=O.bodies[-1].id
  plateF=O.forces.f(Lastnum)[2]
  print ("Load= %.5f"%(plateF))
  h=O.bodies[-1].state.pos[2]
  print ("h=%.5f"%(h))

def save():
 lpnum=O.bodies[-1].id
 Force=O.forces.f(lpnum)[2]
 Displ=O.bodies[lpnum].state.pos[2]
 plot.addData(Fz=Force,Dis=Displ)
 plot.saveDataTxt('Sphere.txt.bz2',vars=('Dis','Fz'))
 export.textExt('/tmp/SphereCompact.dat')

O.dt=1E-5
**************************************************************************************************************************
    According to my preliminary calculations: under a smaller load (50kPa), the loading plate sank into the sand, which implied that the sand strength is too small. However, in the model, the friction coefficient of particles have raised to a surprisingly high level.
   Now, my questions are as follows:
    1. For the simulation, for the sand particles contact model, I have chosen Cohfrictmat to to simulate the anti-rotation effect between unbonded sand particles, but I am not quite sure if the relevant settings are correct. （May I ask how to check whether the anti-rotational moment is generated between the sand particles）
    2. For the calculation, is there an fatal error in the enlargement of the particles? According to the relevant literature, I found that after the model and particle size greater than 30, it seems that the particle size has less effect on the strength, but here the size ratio of loading plate and particle does not meet this requirement. (For this, I doubled the loading plate size, but again, under a 50 kPa load, the loading plate sank into the sandy soil.)

   Any help and advice would be greatly appreciated!

[1] Behavior of Geocell-Reinforced Sand Under a Vertical Load. https://journals.sagepub.com/doi/abs/10.3141/2045-11