yade-users team mailing list archive

[Horation <question707087@xxxxxxxxxxxxxxxxxxxxx>]

New question #707087 on Yade:
https://answers.launchpad.net/yade/+question/707087

I was conducting a triaxial test on the argillite, which contained a matrix and two inclusions, to investigate its small strain stiffness. I tried to use spherical grains and JCF model to simulate this kind of rock. And I seperated the particle into three groups and assigned them different material properties (Young's modulous and Poisson's ratio). But the results were not desired. Because the small strain stiffness of the sample is almost equivalent to the sample containg only one material, which is confusing. Theoretically, the change of materials in the sample will affect the properties of the rock. How can I realize this goal? Could you give me some advice?

from yade import ymport, utils , plot

O.load('isocompact200kpa2.yade.bz2')
O.engines[6].dead=True
#---------------- ENGINES ARE DEFINED HERE
try:
    os.mkdir('post')
except:
    pass

PACKING='121_1k.spheres'
OUT=PACKING+'_1MPa_r0.02'
iterMax=50000 # maximum number of iterations
saveVTK=iterMax/5 # Vtk files record interval

YOUNG2=100e9
YOUNG3=50e9
ALPHA2=0.06
ALPHA3=0.27

bodies2=[b for b in O.bodies if 5.0<b.id<178.0]
bodies3=[b for b in O.bodies if 177.0<b.id<349.0]

for b in bodies2:
 b.mat.young=YOUNG2
 b.mat.poisson=ALPHA2

for b in bodies3:
 b.mat.young=YOUNG3
 b.mat.poisson=ALPHA3

confinement=-7e6
strainRate=-0.0001
triax=O.engines[3]
#### triaxial Engine
triax.stressMask=5
triax.goal1=confinement
triax.goal2=strainRate
triax.goal3=confinement
triax.max_vel=1

#### simulation is defined here (DEM loop, interaction law, servo control, recording, etc...)
O.engines=[
        ForceResetter(),
        InsertionSortCollider([Bo1_Box_Aabb(),Bo1_Sphere_Aabb(aabbEnlargeFactor=-1.,label='Saabb')]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=-1.,label='SSgeom'),Ig2_Box_Sphere_ScGeom()],
  [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1,label='interactionPhys')],
  [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(recordCracks=True,Key='data',label='interactionLaw')]
 ),
        triax,
        GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.5,defaultDt=0.1*utils.PWaveTimeStep()),
        NewtonIntegrator(damping=0.5,label="newton"),
        PyRunner(command='addPlotData()',iterPeriod=10,label='record'),
        VTKRecorder(iterPeriod=int(5000),initRun=True,fileName='./post/p1-',recorders=['spheres','jcfpm','cracks'],Key='data',label='vtk')
]

#### custom recording functions
#tensCks=shearCks=0

# if you want to plot during simulation

#---------------- SIMULATION STARTS HERE

#### manage interaction detection factor during the first timestep and then set default interaction range (intRadius=1)

def addPlotData():

    global ini_e22a
    global ini_e22b

    e22=-triax.strain[1]
    unb = unbalancedForce()
    mStress = -(triax.stress(triax.wall_right_id)[0]+triax.stress(triax.wall_top_id)[1]+triax.stress(triax.wall_front_id)[2])/3.0
    volume = (O.bodies[1].state.pos[0]-O.bodies[0].state.pos[0])*(O.bodies[3].state.pos[1]-O.bodies[2].state.pos[1])*(O.bodies[5].state.pos[2]-O.bodies[4].state.pos[2])
    Porosity=1-sum(4*pi*b.shape.radius**3/3 for b in O.bodies if isinstance(b.shape,Sphere))/volume

    plot.addData(e11=-triax.strain[0], e22=-triax.strain[1], e33=-triax.strain[2],
     ev=-triax.strain[0]-triax.strain[1]-triax.strain[2],
   s11=-triax.stress(triax.wall_right_id)[0],
   s22=-triax.stress(triax.wall_top_id)[1],
   s33=-triax.stress(triax.wall_front_id)[2],
                 ub=unbalancedForce(),
                 dstress=(-triax.stress(triax.wall_top_id)[1])-(-triax.stress(triax.wall_right_id)[0]),
                 disx=triax.width,
                 disy=triax.height,
                 disz=triax.depth,
   i=O.iter,
                 porosity=Porosity,
                 cnumber=avgNumInteractions(),
                 tc=interactionLaw.nbTensCracks,
                 sc=interactionLaw.nbShearCracks,
                 te=interactionLaw.totalTensCracksE,
                 se=interactionLaw.totalShearCracksE,
                 )

    #if (e22 - ini_e22a) > 0.001 and e22 < 0.05:
    # O.save('save'+'{:.4f}'.format(e22)+'yade.bz2')
    # ini_e22a = e22

    #if (e22 - ini_e22b) > 0.00005:
    print ('unbalanced force: %f, mean stress: %f, s11: %f, s22: %f, s33: %f, coordination number: %f, porosity: %f' %(unb,mStress,-triax.stress(triax.wall_right_id)[0],-triax.stress(triax.wall_top_id)[1],-triax.stress(triax.wall_front_id)[2],avgNumInteractions(),Porosity))
        #Output()
        #ini_e22b = e22
    plot.saveDataTxt('loadinglog3.txt.bz2')

plot.plots={'e22':('s22',None,'ev')}
plot.plot()
#O.run(iterMax)

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