yade-users team mailing list archive

[James X <question314841@xxxxxxxxxxxxxxxxxxxxx>]

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

Hi,
I am still learning many of the features in Yade and I think I must be doing something very wrong with TriaxialStressController. I want to simulate a 2D uniaxial tension experiment on this sytem, ideally by holding stress constant in one direction while increasing strain in another direction.  Below is a stripped down 2D simulation that is not working correctly, which does work if I comment out the "triax" section and uncomment the unaxialstrainer line.

I have set the bitmask to zero, as I would for now just like to reproduce the unaxial strain behavior. Once I get that working, I intend to add stress into the system by manipulating bitmask and the goal1 parameter, but for now I just want to see that I can reproduce uniaxial behavior. I tried looking at the session 1 tutorial example, but am still confused about how I can apply it for my problem.

What do I need to do in order to see the uniaxial strain behavior with TriaxialStressController?

--------------------------- BEGIN CODE ---------------------------

import numpy as np
import time

from yade import pack
from yade import geom
from yade import polyhedra_utils
from yade import utils
from yade import linterpolation
from yade import export
from yade import ymport

#other random parameters
start = time.time()

#this function just makes a 2D sphere mesh and exports it as <FILENAME>.sphere
def Make_Packed_file(filename, xSize, ySize, radius, rRelFuzz_value, nbIter=3000):
    #filename = filename #output a msh file
    sp=pack.SpherePack()
    sp.makeCloud(minCorner=(0,0,0),maxCorner=(xSize,ySize*2.5,0),rMean=radius, rRelFuzz = 0)
    sp.toSimulation()
    for b in O.bodies:
        b.state.blockedDOFs = 'zXY'
    box=O.bodies.append(utils.geom.facetBox((xSize/2,1.25*ySize,0),(xSize/2,1.25*ySize,radius),wallMask=63)) # the same box of makeCloud

    O.engines=[
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Wall_Aabb()]),
        InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
        [Ip2_FrictMat_FrictMat_FrictPhys()],
        [Law2_ScGeom_FrictPhys_CundallStrack()]
        ),
        NewtonIntegrator(gravity=(0,-2.0,0),damping=.1),
        PyRunner(iterPeriod=500,initRun=False,command='print("elapsed time is: " + str(time.time() - start) ) '),
        DomainLimiter(lo=(-0.01,-0.01,-0.01),hi=(xSize,ySize,1.1*radius),iterPeriod=nbIter-1)
    ]
    O.dt=utils.PWaveTimeStep()
    O.stopAtIter=nbIter
    O.run()
    O.wait()
    export.text(filename+'.spheres')


###random parameters
radius = 3
smoothContact=True
jointFrict=radians(20)
jointDil=radians(0)
new_scale = 1.0 
xSize = 250;
ySize = 150;
rRelFuzz_value = 0;

#actually make this packed mesh
Make_Packed_file("simplified_packed_2D_mesh", xSize, ySize, radius, rRelFuzz_value)

O.reset() #reset everything, just in case


#defin the material
def mat_1(): return JCFpmMat(type=1,young=30.0e9,poisson=0.3,frictionAngle=radians(30),density=3000,tensileStrength=1.23e8,cohesion=1.23e8,jointNormalStiffness=1e7,jointShearStiffness=1e7,jointCohesion=1e6,jointFrictionAngle=radians(20),jointDilationAngle=0.0) 


O.materials.append((mat_1()))

#import the mesh created in the previous simulation
import_spheres = ymport.text("simplified_packed_2D_mesh.spheres",scale=new_scale,shift=Vector3(0,0,0),material=mat_1)
O.bodies.append(import_spheres)


#add in wall information, from uniax.py
young = 90e9
mn,mx=Vector3(0,0,0),Vector3(xSize*new_scale,1.01*ySize*new_scale,2.00001*radius) # corners of the initial packing
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls'))
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)


#from uniax.py, original needed for original uniaxialstrainer
bb=uniaxialTestFeatures()
negIds,posIds,axis,crossSectionArea=bb['negIds'],bb['posIds'],bb['axis'],bb['area']

strainRateTension=.05
setSpeeds=True


interactionRadius=1.20 # to set initial contacts to larger neighbours and a bit further


triax=TriaxialStressController(
    maxMultiplier=1.0002,
    finalMaxMultiplier=1.002,
    thickness = 0,
    stressMask = 0, #Bitmask, if imposed goal`s are stresses (0 for none, 7 for all, 1 for direction 1, 5 for directions 1 and 3, etc. :ydefault:`7)
    internalCompaction=False,
    goal1 = 0,
    goal2 = strainRateTension, # positive is tension, negative is compression
    goal3 = 0
)



O.engines=[

	ForceResetter(),
	InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=interactionRadius,label='is2aabb'),]),
	InteractionLoop(
		[Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=interactionRadius,label='ss2d3dg')],
		[Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1,label='interactionPhys')],
		[Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(smoothJoint=smoothContact,label='interactionLaw')]
	),
	GlobalStiffnessTimeStepper(timestepSafetyCoefficient=0.8),
	NewtonIntegrator(damping=0.2,gravity=(0.,0,0.)),
	triax,
	#TriaxialStateRecorder(iterPeriod=100,file='WallStresses'+table.key),
	PyRunner(iterPeriod=5000,initRun=False,command='print("5000 runs in fracture run, elapsed time is: " + str(time.time() - start) + ", goal2 is: " + str(triax.goal2)) '),
        #UniaxialStrainer(strainRate=strainRateTension,axis=axis,asymmetry=0,posIds=posIds,negIds=negIds,crossSectionArea=crossSectionArea,blockDisplacements=False,blockRotations=False,setSpeeds=setSpeeds,label='strainer'),
]

#### time step definition (low here to create cohesive links without big changes in the assembly)
O.dt=0.05*utils.PWaveTimeStep()

#### set cohesive links with interaction radius>=1
O.step();

#### initializes now the interaction detection factor to strictly 1
ss2d3dg.interactionDetectionFactor=-1.0
is2aabb.aabbEnlargeFactor=-1.0



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