yade-users team mailing list archive

[Anay Joshi <question661250@xxxxxxxxxxxxxxxxxxxxx>]

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

    Status: Answered => Open

Anay Joshi is still having a problem:
Hi Luc,
Thanks a lot for this descriptive information. I have a few more questions
to ask,
1. I read your paper, ' A DEM model for soft and hard rocks: Role of grain
interlocking on strength' in which you have modeled the granite. However, I
am modeling a concrete cylinder in my work. I have read in several
literatures that they used the bonded particle model (Cundall,2004). My
question is, can I use JCFpm for modeling concrete?

2. I am running the simulation on the remote servers and I would like to
know that is there any way to see the time-dependent animation of the whole
process of my simulation? Also, can you suggest me how to see the
streamlines of the DE with the yade?


Thanks,
Anay

On Wed, Jan 24, 2018 at 9:37 PM, Luc Scholtès <
question661250@xxxxxxxxxxxxxxxxxxxxx> wrote:

> Your question #661250 on Yade changed:
> https://answers.launchpad.net/yade/+question/661250
>
>     Status: Open => Answered
>
> Luc Scholtès proposed the following answer:
> Wow... so many lines.
>
> I would suggest that you use the following first batch of lines to get
> used to yade first. If you copy paste these lines in a file and name it
> "uniaxialCompression.py", you should then be able to run it with "yade
> uniaxialCompression.py". Basically, it simulates a uniaxial compression
> with the JCFPM model. Then, if you want to plot some data from the
> simulation, you can copy the second batch of lines, copy paste them in a
> file name "plotUniaxialCompression.py" and then launch it with "python
> plotUniaxialCompression.py".
>
> Then, once you know how to run and exploit this simulation, you'll be
> able to modify the script with the contact law that suits you the most.
> If not, please ask a simple question with minimum lines and clear
> (short) title.
>
> Luc
>
> --- uniaxial compression:
>
> # -*- coding: utf-8 -*-
> from yade import pack, plot
>
> ################# SIMULATIONS DEFINED HERE
>
> #### packing (previously constructed)
> OUT='compressionTest_JCFPM'
>
> #### Simulation Control
> rate=-0.01 #deformation rate
> iterMax=10000 # maximum number of iterations
> saveVTK=2000 # saving output files for paraview
>
> #### Material microproperties
> intR=1.1 # allows near neighbour interaction (can be adjusted for every
> packing)
> DENS=2500 # could be adapted to match material density:
> dens_DEM=dens_rock*(V_rock/V_particles)=dens_rock*1/(1-poro_DEM) ->
> packing porosity as to be computed?
> YOUNG=20e9
> FRICT=7
> ALPHA=0.1
> TENS=1e6
> COH=1e6
>
> #### material definition
> def sphereMat(): return JCFpmMat(type=1,density=DENS,
> young=YOUNG,poisson=ALPHA,frictionAngle=radians(FRICT),
> tensileStrength=TENS,cohesion=COH)
>
> #### create the specimen
> pred=pack.inCylinder((0,0,0),(0,1,0),0.25)
> O.bodies.append(pack.regularHexa(pred,radius=0.025,
> gap=0.,material=sphereMat))
>
> R=0
> Rmax=0
> nbSpheres=0.
> for o in O.bodies:
>  if isinstance(o.shape,Sphere):
>    nbSpheres+=1
>    R+=o.shape.radius
>    if o.shape.radius>Rmax:
>      Rmax=o.shape.radius
> Rmean=R/nbSpheres
>
> print 'nbSpheres=',nbSpheres,' | Rmean=',Rmean
>
> #### boundary condition (see utils.uniaxialTestFeatures
> bb=utils.uniaxialTestFeatures()
> negIds,posIds,longerAxis,crossSectionArea=bb['negIds'],
> bb['posIds'],bb['axis'],bb['area']
>
> ################# ENGINES DEFINED HERE
>
> O.engines=[
>         ForceResetter(),
>         InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=
> intR,label='Saabb')]),
>         InteractionLoop(
>                 [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=
> intR,label='SSgeom')],
>                 [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(
> cohesiveTresholdIteration=1,label='interactionPhys')],
>                 [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(
> recordCracks=True,Key=OUT,label='interactionLaw')]
>         ),
>         UniaxialStrainer(strainRate=rate,axis=longerAxis,
> asymmetry=0,posIds=posIds,negIds=negIds,crossSectionArea=crossSectionArea,
> blockDisplacements=1,blockRotations=1,setSpeeds=0,
> stopStrain=0.1,dead=1,label='strainer'),
>         GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.5,
> defaultDt=utils.PWaveTimeStep()),
>         NewtonIntegrator(damping=0.4,label='newton'),
>         PyRunner(iterPeriod=int(100),initRun=True,command='
> recorder()',label='data'),
>         VTKRecorder(iterPeriod=int(saveVTK),initRun=True,
> fileName=OUT+'-',recorders=['spheres','jcfpm','cracks'],
> Key=OUT,label='vtk')
> ]
>
> ################# RECORDER DEFINED HERE
>
> def recorder():
>     yade.plot.addData({'i':O.iter,
>                        'eps':strainer.strain,
>                        'sigma':strainer.avgStress,
>                        'tc':interactionLaw.nbTensCracks,
>                        'sc':interactionLaw.nbShearCracks,
>                        'te':interactionLaw.totalTensCracksE,
>                        'se':interactionLaw.totalShearCracksE,
>                        'unbF':utils.unbalancedForce()})
>     plot.saveDataTxt(OUT)
>
> # if you want to plot during simulation
> plot.plots={'i':('sigma')}
> #plot.plot()
>
> ################# PREPROCESSING
>
> #### manage interaction detection factor during the first timestep and
> then set default interaction range ((cf. A DEM model for soft and hard
> rock, Scholtes & Donze, JMPS 2013))
> O.step();
> ### initializes the interaction detection factor
> SSgeom.interactionDetectionFactor=-1.
> Saabb.aabbEnlargeFactor=-1.
>
> #### coordination number verification
> numSSlinks=0
> numCohesivelinks=0
> for i in O.interactions:
>     if isinstance(O.bodies[i.id1].shape,Sphere) and
> isinstance(O.bodies[i.id2].shape,Sphere):
>       numSSlinks+=1
>     if i.phys.isCohesive :
>       numCohesivelinks+=1
>
> print "nblinks=", numSSlinks, " | nbCohesivelinks=", numCohesivelinks, "
> || Kcohesive=", 2.0*numCohesivelinks/nbSpheres
>
> ################# SIMULATION REALLY STARTS HERE
> strainer.dead=0
> O.run(iterMax)
>
> --- plot uniaxial test data:
>
> # -*- coding: utf-8 -*-
> from pylab import *
>
> ### processing function
> def store(var,textFile):
>     data=loadtxt(textFile,skiprows=1)
>     it=[]
>     eps=[]
>     sig=[]
>     tc=[]
>     sc=[]
>     te=[]
>     se=[]
>     ubf=[]
>     for i in range(0,len(data)):
>       it.append(float(data[i,1]))
>       eps.append(abs(float(data[i,0])))
>       sig.append(abs(float(data[i,4])))
>       tc.append(float(data[i,5]))
>       sc.append(float(data[i,2]))
>       te.append(float(data[i,6]))
>       se.append(float(data[i,3]))
>       ubf.append(float(data[i,7]))
>     var.append(it)
>     var.append(eps)
>     var.append(sig)
>     var.append(tc)
>     var.append(sc)
>     var.append(te)
>     var.append(se)
>     var.append(ubf)
>
> ### data input
> dataFile1='compressionTest_JCFPM'
>
> a1=[]
> store(a1,dataFile1)
>
> rcParams.update({'legend.numpoints':1,'font.size':
> 20,'axes.labelsize':25,'xtick.major.pad':10,'ytick.major.
> pad':10,'legend.fontsize':20})
> lw=2
> ms=10
>
> ### plots
>
> figure(0,figsize=(10,10))
> ax1=subplot(1,1,1)
> xlabel('iteration [-]')
> ax1.plot(a1[0],[x/1e6 for x in a1[2]],'-k',linewidth=lw)
> ylabel(r'$\sigma_1$ [MPa]')
> ax2 = ax1.twinx()
> ax2.plot(a1[0],a1[5],'-r',linewidth=lw)
> ylabel('unbForce [-]')
>
> figure(1,figsize=(10,10))
> ax1=subplot(1,1,1)
> xlabel(r'$\epsilon_1$ [millistrain]')
> plot([x*1e3 for x in a1[1]],[x/1e6 for x in a1[2]],'-k',linewidth=lw)
> ylabel(r'$\sigma_1$ [MPa]')
>
> figure(2,figsize=(10,10))
> ax1=subplot(1,1,1)
> xlabel(r'$\epsilon_1$ [millistrain]')
> ax1.plot([x*1e3 for x in a1[1]],[x/1e6 for x in a1[2]],'-k',linewidth=lw)
> ylabel(r'$\sigma_1$ [MPa]')
> ax2 = ax1.twinx()
> ax2.plot([x*1e3 for x in a1[1]],a1[3],'-b',linewidth=lw)
> ax2.plot([x*1e3 for x in a1[1]],a1[4],'-r',linewidth=lw)
> ylabel('cumulative number of microcracks [-]')
> legend(('tensile','shear'),loc=2)
>
> figure(3,figsize=(10,10))
> ax1=subplot(1,1,1)
> xlabel(r'$\epsilon_1$ [millistrain]')
> ax1.plot([x*1e3 for x in a1[1]],[x/1e6 for x in a1[2]],'-k',linewidth=lw)
> ylabel(r'$\sigma_1$ [MPa]')
> ax2 = ax1.twinx()
> ax2.plot([x*1e3 for x in a1[1]],[x+y for x,y in zip(a1[5],a1[6])],'-r',
> linewidth=lw)
> ylabel('released energy [J]')
>
> ### show or save
> show()
>
> --
> If this answers your question, please go to the following page to let us
> know that it is solved:
> https://answers.launchpad.net/yade/+question/661250/+confirm?answer_id=11
>
> If you still need help, you can reply to this email or go to the
> following page to enter your feedback:
> https://answers.launchpad.net/yade/+question/661250
>
> You received this question notification because you asked the question.
>

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