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

Hi, I'm new in Yade. I wanted to do a triaxial test on a cylindrical specimen. 
I found a code on github and modified the code for my specimen (the original 
code creates a spherical pack and do triaxial test but I created a cylindrical 
sample  of desired porosity and imported the coordinates to this code). I ran 
the code but the cylindrical membrane that was supposed to be created around 
the spherical pack is created away from the pack. Can somebody help me?

Here is the link to the actual code: 
http://bazaar.launchpad.net/~yade-pkg/yade/git-trunk/view/head:/examples/concrete/triax.py
Here is the code:

# -*- encoding=utf-8 -*-
from __future__ import print_function
################################################################################
#
# Triaxial test. Axial strain rate is prescribed and transverse prestress.
# Test is possible on prism or cylinder
# An independent c++ engine may be created from this script in the future.
#
################################################################################
from builtins import range
from yade import ymport, plot
from yade import pack, plot
import os

# default parameters or from table
readParamsFromTable(noTableOk=True,
        # material parameters
        young = 10e7,
        poisson = .3,
        frictionAngle = 0.15,
        sigmaT = 1.5e6,
        
        # prestress
        preStress = -3e6,
        # axial strain rate
        strainRate = -100,

        # assamlby parameters
        rParticle = 0.0005,
        width = 5,
        height = 9.85,
        bcCoeff = 5,

        # facets division
        nw = 48,
        nh = 30,

        # output specifications
        fileName = 'test',
        exportDir = '/tmp',
        runGnuplot = False,
        runInGui = True,
)
from yade.params.table import *

# materials
sphereMat = 
O.materials.append(CohFrictMat(young=young,poisson=poisson,frictionAngle=frictionAngle,alphaKr=0.25,alphaKtw=0,etaRoll=0.005,etaTwist=0,normalCohesion=5e6,shearCohesion=5e6,momentRotationLaw=True,density=2530))

frictMat = O.materials.append(FrictMat(
        young=young,poisson=poisson,frictionAngle=frictionAngle
))

# spheres
sp= ymport.textExt('DensepackRD70-80',format='x_y_z_r',shift= 
Vector3(0,0,0.5*height), scale=1.0, material=sphereMat)
spheres = O.bodies.append(sp)

#pred = pack.inCylinder((0,0,0),(0,0,height),.5*width) if testType=='cyl' else 
pack.inAlignedBox((-.5*width,-.5*width,0),(.5*width,.5*width,height)) if 
testType=='cube' else None
#sp=SpherePack()
#sp = 
pack.randomDensePack(pred,spheresInCell=2000,radius=rParticle,memoizeDb='/tmp/triaxTestOnCylinder.sqlite',returnSpherePack=True)
#spheres=sp.toSimulation(color=(0,1,1))

# bottom and top of specimen. Will have prescribed velocity
bot = [O.bodies[s] for s in spheres if 
O.bodies[s].state.pos[2]<rParticle*bcCoeff]
top = [O.bodies[s] for s in spheres if 
O.bodies[s].state.pos[2]>height-rParticle*bcCoeff]
vel = strainRate*(height-rParticle*2*bcCoeff)
top_limit = 0
top_id = 0
for s in top:
        if s.state.pos[2]>=top_limit:
                top_limit = s.state.pos[2]
                top_id = s.id
bot_limit = height
bot_id = 0
for s in bot:
        if s.state.pos[2]<=bot_limit:
                bot_limit = s.state.pos[2]
                bot_id = s.id
# facets
facets = []
rCyl2 = .5*width / cos(pi/float(nw))
for r in range(nw):
        for h in range(nh):
                v1 = Vector3( rCyl2*cos(2*pi*(r+0)/float(nw)), 
rCyl2*sin(2*pi*(r+0)/float(nw)), height*(h+0)/float(nh) )
                v2 = Vector3( rCyl2*cos(2*pi*(r+1)/float(nw)), 
rCyl2*sin(2*pi*(r+1)/float(nw)), height*(h+0)/float(nh) )
                v3 = Vector3( rCyl2*cos(2*pi*(r+1)/float(nw)), 
rCyl2*sin(2*pi*(r+1)/float(nw)), height*(h+1)/float(nh) )
                v4 = Vector3( rCyl2*cos(2*pi*(r+0)/float(nw)), 
rCyl2*sin(2*pi*(r+0)/float(nw)), height*(h+1)/float(nh) )
                f1 = facet((v1,v2,v3),color=(0,0,1),material=frictMat)
                f2 = facet((v1,v3,v4),color=(0,0,1),material=frictMat)
                facets.extend((f1,f2))

O.bodies.append(facets)
mass = O.bodies[0].state.mass
for f in facets:
        f.state.mass = mass
        f.state.blockedDOFs = 'XYZz'

# plots 
plot.plots = { 'e':('s',), }
def plotAddData():
        f1 = sum(O.forces.f(b.id)[2] for b in top)
        f2 = sum(O.forces.f(b.id)[2] for b in bot)
        f = .5*(f2-f1)
        s = f/(pi*.25*width*width)
        e = (top[0].state.displ()[2] - bot[0].state.displ()[2]) / 
(height-rParticle*2*bcCoeff)
        plot.addData(
                i = O.iter,
                s = s,
                e = e,
        )

# apply prestress to facets
def addForces():
        for f in facets:
                n = f.shape.normal
                a = f.shape.area
                O.forces.addF(f.id,preStress*a*n)

# stop condition and exit of the simulation
def stopIfDamaged(maxEps=5e-3):
        extremum = max(abs(s) for s in plot.data['s'])
        s = abs(plot.data['s'][-1])
        e = abs(plot.data['e'][-1])
        if O.iter < 1000 or s > .5*extremum and e < maxEps:
                return
        f = os.path.join(exportDir,fileName)
        print('gnuplot',plot.saveGnuplot(f,term='png'))
        if runGnuplot:
                import subprocess
                os.chdir(exportDir)
                subprocess.Popen(['gnuplot',f+'.gnuplot']).wait()
        print('Simulation finished')
        O.pause()
        #sys.exit(0) # results in some threading exception

O.dt=.5*utils.PWaveTimeStep()
enlargeFactor=1.5
O.engines=[
        ForceResetter(),
        InsertionSortCollider([
                Bo1_Sphere_Aabb(aabbEnlargeFactor=enlargeFactor,label='bo1s'),
                Bo1_Facet_Aabb()
        ]),
        InteractionLoop(
                [
                        
Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=enlargeFactor,label='ss2d3dg'),
                        Ig2_Facet_Sphere_ScGeom(),
                ],
                [
                        
Ip2_CpmMat_CpmMat_CpmPhys(cohesiveThresholdIter=O.iter+5),
                        Ip2_FrictMat_CpmMat_FrictPhys(),
                        Ip2_FrictMat_FrictMat_FrictPhys(),
                ],
                [
                        Law2_ScGeom_CpmPhys_Cpm(),
                        Law2_ScGeom_FrictPhys_CundallStrack(),
                ],
        ),
        PyRunner(iterPeriod=1,command="addForces()"),
        NewtonIntegrator(damping=.3),
        CpmStateUpdater(iterPeriod=50,label='cpmStateUpdater'),
        PyRunner(command='plotAddData()',iterPeriod=100),
        PyRunner(iterPeriod=50,command='stopIfDamaged()'),
]

# run one step
O.step()

# reset interaction detection enlargement
bo1s.aabbEnlargeFactor=ss2d3dg.interactionDetectionFactor=1.0

# initialize auto-updated plot
if runInGui:
        plot.plot()
        try:
                from yade import qt
                renderer=qt.Renderer()
                # uncomment following line to exagerate displacement
                #renderer.dispScale=(100,100,100)
        except:
                pass

# run
O.run()


-- 
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