About uniaxial compression test simulation by CohFrictMat in Yade
Purpose: I am trying to simulate uniaxial compression of sea ice using the CohFrictMat material. The python file is as follows. Among them, I used some papers and my own experience to set parameter values, such as normalCohesion and shearCohesion. Problems: However, the simulation results are quite different from the experiment, and the sea ice specimen in the simulation is completely collapsed. In theory, under uniaxial compression of conventional materials, there will be 45 degree oblique section failure phenomenon, but the simulation has no inclined section failure phenomenon or trend. Problem wanted to solve: How can I improve the python file to make the simulation more realistic?
#######PYTHON FILE############ (a_zuobiao.csv, a_top.mesh and a_boot.mesh are not included in this post,I can provide it if needed
from yade import pack
from numpy import genfromtxt
from yade import utils
from yade import export
from yade import ymport
from yade import wrapper
from yade import plot
import csv
import numpy as np
import gts, os.path, locale
############################################
######### 0.定义变量
############################################
rho_ice = 920
gravity = 9.8
length_beam = 1.0
h1_beam = 0.1
h2_beam = 0.1
loading_rate = 0.02
############################################
######### 1.定义材料参数
############################################
#1.1 冰
matice = O.materials.append(CohFrictMat(isCohesive=True, frictionAngle=radians(15), density=rho_ice, poisson=0.33, young=1e9, normalCohesion=1e6, shearCohesion=1e6,label="cohmat"))
Matice=O.materials[matice]
#1.2 钢材
matsteel=O.materials.append(FrictMat(young=205e9,poisson=0.3,density=7580,frictionAngle=radians(30)))
Matsteel=O.materials[matsteel]
#1.3 生成冰梁离散元模型
coordinates=numpy.genfromtxt("a_zuobiao.csv",delimiter=",")
pack_ice=[sphere((x,y,z),radius=0.0035,material=Matice,color=(1,1,0)) for x,y,z in coordinates]
PACK_ICE=O.bodies.append(pack_ice)
#1.4 上加载压头
top=ymport.gmsh('a_top.mesh',color=(0.,0.42,0.42))
TOP=O.bodies.append(top)
#1.5 下加载压头
bott=ymport.gmsh('a_bott.mesh',color=(0.,0.42,0.42))
BOTT=O.bodies.append(bott)
############################################
######### 2.定义全局Engines
############################################
O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=1.1), Bo1_Facet_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom6D(interactionDetectionFactor=1.1), Ig2_Facet_Sphere_ScGeom()],
[Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True), Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_ScGeom6D_CohFrictPhys_CohesionMoment(),Law2_ScGeom_FrictPhys_CundallStrack()]
),
NewtonIntegrator(gravity=(0, 0, -gravity), damping=0.2),
TranslationEngine(translationAxis=[0,0,-1],velocity=loading_rate, ids=TOP,label='line_velocity'),
PyRunner(command='compressionStrength2displace()',iterPeriod=1000),
VTKRecorder(iterPeriod=2000,recorders=['all'],fileName='out-'),
ForceRecorder(iterPeriod=1000,ids=TOP,file='force.csv',label='force')
]
############################################
######### 3.自定义函数
############################################
#3.1 将加载载荷转换为压缩强度(Mpa-mm)
def compressionStrength2displace():
top_force = sum(O.forces.f(b)[2] for b in TOP)
conversion_coeff = 1/(h1_beam*h2_beam)
comp_strengh = top_force*conversion_coeff/1000000
displace = O.dt*O.iter*loading_rate*1000
with open('a_compStr2dis.txt','a') as file0:
print('%s'%displace,'%s'%comp_strengh,file=file0)
return
############################################
######### 4.控制仿真参数
############################################
#timestep
O.dt = 0.1*PWaveTimeStep() #0.00022886885410853175
print('delta_T',O.dt)
# save the simulation, so that it can be reloaded later, for experimentation
#now continue simulation
O.run(2000000,True)