Steel bar 3-point bending test

I'm trying to simulate a 3-point bending test in order to calibrate/analyze CohFrictMat parameters with the purpose of studying a steel behaviour. The Euler-Bernoulli beam theory is used to calculate the deflection in the middle point.
"Fitted" material parameters should give a bar displacement as close as possible to that value, however the simulated deflection is way smaller than what theory suggests.
I've tried changing all parameters, also in a wide range, but it didn't help.
Someone have any suggestion? Thanks. (Below is attached the code)
-------------------------------------------------------------------------------------------------------

# -*- coding: utf-8 -*-
"""
Created on Fri Sep 18 10:00:32 2020

@author: antonio
"""

#from builtins import range
from yade import plot, pack
from yade.gridpfacet import *
from numpy import linspace

#### Parameters ####
L=.25 # length of the bar
n=7 # number of nodes used to generate the bar
r=0.008/2. # radius of the bar element
I=3.14*(r**4/4) # moment of inertia

#### Engines ####
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),
                           Bo1_GridConnection_Aabb()]),
 InteractionLoop(
  [Ig2_GridNode_GridNode_GridNodeGeom6D(),
         Ig2_Sphere_Sphere_ScGeom(),
         Ig2_Sphere_GridConnection_ScGridCoGeom()],
  [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True,setCohesionOnNewContacts=False)],
  [Law2_ScGeom6D_CohFrictPhys_CohesionMoment(useIncrementalForm=True,always_use_moment_law=True),
         Law2_GridCoGridCoGeom_FrictPhys_CundallStrack(),
         Law2_ScGridCoGeom_CohFrictPhys_CundallStrack()]
 ),
 NewtonIntegrator(gravity=(0.,0.,0.),damping=0.2,label='newton'),
# GlobalStiffnessTimeStepper(timestepSafetyCoefficient=0.1,label='ts'),
 PyRunner(initRun=True,iterPeriod=50,command='addData()',label='adddata'),
 PyRunner(initRun=False,iterPeriod=50,command='checkStop()',label='checkstop')
    ]

#### Material ####
b_young=2.2e9
b_normalCohesion=1.6e9
b_shearCohesion=1.6e9

b_alphaKr=1e8 #Dimensionless rolling stiffness
b_alphaKtw=1e8 #Dimensionless twisting stiffness
b_etaRoll=1 #Dimensionless bending strength
b_etaTwist=1 #Dimensionless twisting strength

O.materials.append(CohFrictMat(young=b_young,poisson=0.3,density=7.85e3,frictionAngle=radians(0.),normalCohesion=b_normalCohesion,shearCohesion=b_shearCohesion,
                               alphaKr=b_alphaKr,alphaKtw=b_alphaKtw,etaRoll=b_etaRoll,etaTwist=b_etaTwist,
                               momentRotationLaw=True,label='mat1'))
#### Nodes and connections ####
nodesIds=[]
for i in linspace(0,L,n):
   nodesIds.append(O.bodies.append(gridNode([i,0,0],r,wire=False,fixed=False,material='mat1',color=[1,0,0])))

for i,j in zip( nodesIds[:-1], nodesIds[+1:]):
   O.bodies.append(gridConnection(i,j,r,color=[1,1,1],material='mat1'))

#### Boundary conditions ####
O.bodies[nodesIds[0]].state.blockedDOFs='xyzXYZ'
O.bodies[nodesIds[-1]].state.blockedDOFs='xyzXYZ'

# DOUBLE-SUPPORTED BEAM BENDING
block=O.bodies.append(sphere((O.bodies[nodesIds[3]].state.pos[0],0,.2),radius=4*r,material='mat1'))

O.bodies[block].state.blockedDOFs='xyzXYZ'

init_pos=O.bodies[nodesIds[3]].state.pos[2]

pushing_vel=-.5
O.bodies[block].state.vel=(0,0,pushing_vel)

#### Set a time step ####
O.dt=1e-06
if O.dt>0.5*PWaveTimeStep():
    print('dt may be too high')

#temporary saving of data inside the simulation
def addData():
    plot.addData(iter=O.iter,
                th_disp=O.time*pushing_vel,
                disp=abs(O.bodies[nodesIds[3]].state.pos[2]-init_pos), #bar deflection
                fblock= abs(O.forces.f(block)[2]),
                fmax=abs(((O.forces.f(block)[2])*L**3)/(48.*b_young*I)) #theoretical deflection = (F*L^3)/(48*E*J)
                )

#real-time plots
plot.plots={'iter':('disp','fmax'),'disp':'fblock'}

#end condition for the simulation
def checkStop():
    if abs(O.bodies[block].state.pos[2])>.5:
        O.pause()
        #to save data in an external *.txt file
        plot.saveDataTxt('dataOutputBending.txt')
        print ('the test is ended')

from yade import qt
plot.plot()
qt.View()
qt.Controller()