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

Hello,

I wish you and your family perfect health.
I am modeling a triaxial test. The problem that I have is that at the beginning 
of the simulation, some particles come out from the walls. Therefore, it 
destroys the shape of micro-strain contours. In addition, the initial stresses 
which are sensed at the walls are so high (more than 1 MPa). I am using YADE 
2020.01a, and Ubuntu 20.04. My python version is 3.8.5. Attached, you can find 
my script. I know it should not be too long, but this is the problem that I 
have. I was wondering if you could give me a solution.
Thank you very much.

Best Regards,

Alireza

from yade import utils, plot
from yade import pack, qt
from datetime import datetime
import math
import pylab
from numpy import *


qtr=qt.Renderer()
qtr.bgColor=(1,1,1)


#==============================================================
#================= define the materials =======================
#==============================================================

O.materials.append(CohFrictMat(normalCohesion= 1e20, shearCohesion= 1e20, 
isCohesive= True, young=6.81e8, density=1377.5, poisson=0.3, frictionAngle= 
0.31, fragile=False, label='aggregate-48'))

O.materials.append(CohFrictMat(normalCohesion= 1e20, shearCohesion= 1e20, 
isCohesive= True, young=6.81e8, density=1532.2, poisson=0.3, frictionAngle= 
0.31, fragile=False, label='aggregate-814'))

O.materials.append(CohFrictMat(normalCohesion= 1e20, shearCohesion= 1e20, 
isCohesive= False, young=2e11, 
density=0, poisson=0.3, frictionAngle= 0, fragile=False, 
label='wall-frictionless'))


#===============================================================
#=================== define walls ============================
#===============================================================

walls=aabbWalls([(0,0,0),(10e-2,10e-2,10e-2)],thickness=0.001,oversizeFactor=1.5,material='wall-frictionless')
wallIds=O.bodies.append(walls)


#===============================================================
#=================== define packing ============================
#===============================================================


nums=['t','t','t']

mats=['aggregate-48','aggregate-814','aggregate-814']

coke=((1.875e-3,3*1000),(0.9475e-3,3*2734),(0.50125e-3,3*24970))

color=((0,0,1),(0,1,0),(1,0,0))

tolerance=[(2e-4),(1e-4),(5e-5)]

for i in range(len(nums)):

     nums[i]=pack.SpherePack()
     
nums[i].makeCloud((5e-3,5e-3,5e-3),(9.95e-2,9.95e-2,9.95e-2),rMean=coke[i][0],rRelFuzz=tolerance[i],num=coke[i][1])
     O.bodies.append([utils.sphere(c,r,material=mats[i],color=color[i]) for c,r 
in nums[i]])


#===============================================================
#=================== define Engine =============================
#===============================================================

O.engines=[
        ForceResetter(),
        
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb(),Bo1_Facet_Aabb()]),
        
InteractionLoop([Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom6D(),Ig2_Facet_Sphere_ScGeom6D()],
        [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
        
[Law2_ScGeom_FrictPhys_CundallStrack(),Law2_ScGeom6D_CohFrictPhys_CohesionMoment()]
        ),
        TriaxialStateRecorder(iterPeriod=10000,file='WallStresses'),
#        VTKRecorder(fileName='3d-vtk-',recorders=['all'],iterPeriod=1000),
#        qt.SnapshotEngine(fileBase='3d-',iterPeriod=200,label='snapshot'),
        NewtonIntegrator(damping=0.4,gravity=[0,0,-9.81]),
        PyRunner(iterPeriod=10000,command='addPlotData()'),
]

#yade.qt.Controller(), yade.qt.View()

O.dt=5e-8
#O.dt=.2*PWaveTimeStep()

def addPlotData():
        plot.addData(exx=-triax.strain[0], eyy=-triax.strain[1], 
ezz=-triax.strain[2], ev=(triax.strain[0]+triax.strain[1]+triax.strain[2]), 
sxx=-triax.stress(triax.wall_right_id)[0], 
syy=-triax.stress(triax.wall_top_id)[1], 
szz=-triax.stress(triax.wall_front_id)[2], mass=sum(b.state.mass for b in 
O.bodies if isinstance(b.shape,Sphere)), 
q=-(triax.stress(triax.wall_front_id)[2]-triax.stress(triax.wall_right_id)[0]), 
p=triax.meanStress, R=triax.stress(triax.wall_front_id)[2]/sigmaIso, 
por=triax.porosity, i=O.iter,**O.energy,total=O.energy.total()),


#=================================================================
#=================  APPLYING CONFINING PRESSURE   ================
#=================================================================
sigmaIso=-5e4
triax=TriaxialStressController(
        maxMultiplier=1.000,  
        finalMaxMultiplier=1.000,  
        thickness = 0,
        stressMask = 7,
        internalCompaction=False,
        max_vel=0.05,
)


O.engines=O.engines+[triax]

triax.goal1=sigmaIso
triax.goal2=sigmaIso
triax.goal3=sigmaIso
triax.wall_back_activated=True


while 1:
        O.run(100, True)
        unb=unbalancedForce()
        
meanS=(triax.stress(triax.wall_right_id)[0]+triax.stress(triax.wall_top_id)[1]+triax.stress(triax.wall_front_id)[2])/3
        print('mean Stress:',triax.meanStress,'porosity:', 
triax.porosity,'Unbalancing:',unb)
        if  triax.porosity<0.37 and 
abs(sigmaIso-triax.meanStress)/abs(sigmaIso)<0.02:
                print('Isotropic strain1:',-triax.strain[0], 'Isotropic strain 
2:',-triax.strain[1], 'Isotropic strain 3:', -triax.strain[2])
                break
                print("###      Isotropic state saved      ###")


triax.depth0=triax.depth
triax.height0=triax.height
triax.width0=triax.width
O.save('RVE3000-sizeDis3-solid-Iso5e4-Rate005-Isopart.yade')



#====================================================================
#=====================   DEVIATORIC LOADING   =======================
#====================================================================

triax.stressMask = 3
triax.goal1 = sigmaIso
triax.goal2 = sigmaIso
triax.goal3 = -0.05
triax.max_vel=0.05


O.trackEnergy=True
O.saveTmp()


#====================================================================
#====================     Micro Strain    ===========================
#====================================================================


TW=TesselationWrapper()
TW.triangulate()
TW.computeVolumes()
TW.volume(10)
TW.setState(0)
step=0
while 1:
        step +=1
        O.run(1000000,True)
        print('Dev strain1:',-triax.strain[0], 'Dev strain 
2:',-triax.strain[1], 'Dev strain 3:', -triax.strain[2])
        TW.setState(1)
        TW.defToVtk("strain_"+str(step)+".vtk")
        if abs(triax.strain[2]) > 0.5:
#               makeVideo(snapshot.snapshots,'3d.mpeg',fps=10,bps=10000)
                print(O.iter)
                print('time is ', O.time)
                print('porosity:',triax.porosity)
                print('strain 1:',-triax.strain[0], 'strain 
2:',-triax.strain[1], 'strain 3:',-triax.strain[2])
                O.save('RVE3000-sizeDis3-solid-Iso5e4-Rate005-Devpart.yade')
                
plot.saveDataTxt('RVE3000-sizeDis3-solid-Iso5e4-Rate005-Devpart.txt')
                break
                O.pause()

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

_______________________________________________
Mailing list: https://launchpad.net/~yade-users
Post to     : yade-users@lists.launchpad.net
Unsubscribe : https://launchpad.net/~yade-users
More help   : https://help.launchpad.net/ListHelp

Reply via email to