yade-users team mailing list archive

[James X <question314841@xxxxxxxxxxxxxxxxxxxxx>]

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

    Status: Answered => Open

James X is still having a problem:
Hello,

Thank you guys so far for your help. I imposed a small stress along the
x and z axis, with a set strain rate, but instead of pulling apart as I
would hope in uniaxial tension, I just see certain elements 'explode'.
Unfortunately, after testing with various strain rates and tensile
stresses, I find that the mesh will first move a bit and then 'break
apart' or do nothing. The behavior i expect to see is that the system
acts in uniaxial tension.

The modified code is (only significant change is in the 'TriaxialStressController' section):
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 = 2
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']


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


strainRateTension=1e3
triax=TriaxialStressController(
    maxMultiplier=1.0002,
    finalMaxMultiplier=1.002,
    thickness = 0,
    stressMask = 5, #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 = -.001,
    goal2 = strainRateTension, # positive is tension, negative is compression
    goal3 = -.001
)


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,
	#PyRunner(iterPeriod=5000,initRun=False,command='print("5000 runs in fracture run, elapsed time is: " + str(time.time() - start) + ", goal2 is: " + str(triax.goal2)) '),
]

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

O.run()

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