yade-users team mailing list archive

[Nicola <question693083@xxxxxxxxxxxxxxxxxxxxx>]

Question #693083 on Yade changed:
https://answers.launchpad.net/yade/+question/693083

Description changed to:
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'

# CANTILIEVER BENDING
#block=O.bodies.append(sphere((O.bodies[nodesIds[-1]].state.pos[0],0,0.3),radius=4*r,material='mat1'))
#joint=O.bodies.append(sphere((O.bodies[nodesIds[0]].state.pos[0],0,-.01),radius=r,material='mat1'))

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

#O.bodies[joint2].state.blockedDOFs='xyzXYZ'
#O.bodies[joint].state.blockedDOFs='xyzXYZ'
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
#                moment = abs(O.interactions[3,4].phys.fragile), 
#                force = abs(O.interactions[3,4].phys.normalForce),
#                force = abs(O.forces.f(nodesIds[3])[2]),
                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)
#               force2 = abs(O.forces.f(nodesIds[2])[2]),
                )

#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()

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.