yade-users team mailing list archive

[hanzj <question443228@xxxxxxxxxxxxxxxxxxxxx>]

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

hello all
I was doing the cohesive triaxial test.
 first ,I  am going to reach the confined state by internalcompaction,  and then compact it to reach the porosity desired, the last step is to do deviatoric loading. 
my system is ubuntu 14.04, the trunk version is Yade 2016-03-28.git-1f77f38.
Apart from the wall using the frictmat, I am also using two kinds of material,but they are both cohefricmat, but have different parameters.but when i turn the "O.engines[2].physDispatcher.functors[1].setCohesionNow = True" on , there exists segmentation fault, if not, it works. and if i use one cohefricmat , it works no matter setCohesionNow is True or False.
any ideas?

thank you very much

Han Zhijian 

my script as follows:

#!/usr/bin/python
# -*- coding: utf-8 -*-
from __future__ import division

from yade import plot,pack,timing
import time, sys, os, copy 
from yade import ymport, utils, pack, export
#from pylab import *
import math
import numpy as np

############################################
###   DEFINING VARIABLES AND MATERIALS   ###
############################################
# default parameters or from table
readParamsFromTable(noTableOk=True, # unknownOk=True,
	targetPorosity=0.45, #the porosity we want for the packing
	num_spheres=3000,# number of spheres

	young=2.4e6,
	poisson=.2,
	density=4800,
	frictionAngle=radians(30),
	normalCohesion=3.e6,
	shearCohesion=3.e6,
	etaRoll=0.1,

	compFricDegree=30, # contact friction during the confining phase
	finalFricDegree=30, # contact friction during the deviatoric loading
	key='_triax_base_', # put you simulation's name here

	rate=-0.02, # loading rate (strain rate)
	damp=0.7, # damping coefficient
	stabilityThreshold=0.01, # we test unbalancedForce against this value in different loops (see below)
)
from yade.params import table
from yade.params.table import *
mn,mx=Vector3(0,0,0),Vector3(1,1,1) # corners of the initial packing
## create materials for spheres and plates
mat1=CohFrictMat(young=young,poisson=poisson,density=density,frictionAngle=frictionAngle,normalCohesion=normalCohesion,shearCohesion=shearCohesion,momentRotationLaw=True,etaRoll=etaRoll,label='cement')
mat_1=O.materials.append(mat1)
mat2=CohFrictMat(young=young,poisson=poisson,density=density,frictionAngle=frictionAngle,normalCohesion=normalCohesion,shearCohesion=shearCohesion,momentRotationLaw=True,etaRoll=etaRoll,label='glass_sphere')
mat_2=O.materials.append(mat2)
O.materials.append(FrictMat(young=2*young,poisson=.25,frictionAngle=0,density=0,label='frictionlessWalls'))

## create walls around the packing
walls=aabbWalls([mn,mx],thickness=0,material='frictionlessWalls')
wallIds=O.bodies.append(walls)

## use a SpherePack object to generate a random loose particles packing
sp=pack.SpherePack()
sp.makeCloud(mn,mx,-1,0,num_spheres,False, 0.95,seed=1) #"seed" make the "random" generation always the same
O.bodies.append([sphere(center,rad,material='cement') for center,rad in sp])
cement=[]
glass_sphere=[]
for b in O.bodies:
	if random.random() < 0.5:
		b.mat = mat1
		b.shape.color = (1,0,0)
		b.state.mass*=mat1.density/mat1.density
		cement.append(b.id)
	else:
		b.mat = mat2
		b.shape.color = (0,1,1)
		b.state.mass*=mat2.density/mat1.density
		glass_sphere.append(b.id)
############################
###   DEFINING ENGINES   ###
############################

triax=TriaxialStressController(
	maxMultiplier=1.+2.4e5/young, # spheres growing factor (fast growth)
	finalMaxMultiplier=1.+2.4e4/young, # spheres growing factor (slow growth)
	thickness = 0,
	stressMask = 7,
	internalCompaction=True, # If true the confining pressure is generated by growing particles
)

newton=NewtonIntegrator(damping=damp)

O.engines=[
	ForceResetter(),
	InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
	InteractionLoop(
		#box-sphere interactions will be the simple normal-shear law, we use ScGeom for them
		[Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom()],
		#Boxes will be frictional (FrictMat), so the sphere-box physics is FrictMat vs. CohFrictMat, the Ip type will be found via the inheritance tree (CohFrictMat is a FrictMat) and will result in FrictPhys interaction physics
		#and will result in a FrictPhys
		[Ip2_FrictMat_FrictMat_FrictPhys(),Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(label="cohesiveIp")],
		#Finally, two different contact laws for sphere-box and sphere-sphere
		[Law2_ScGeom_FrictPhys_CundallStrack(),Law2_ScGeom6D_CohFrictPhys_CohesionMoment(
			useIncrementalForm=True, #useIncrementalForm is turned on as we want plasticity on the contact moments
			always_use_moment_law=False,  #if we want "rolling" friction even if the contact is not cohesive (or cohesion is broken), we will have to turn this true somewhere
			label='cohesiveLaw')]
	),
	GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
	triax,
	newton
]
triax.goal1=triax.goal2=triax.goal3=-10000
while 1:
  O.run(1000, True)
  ##the global unbalanced force on dynamic bodies, thus excluding boundaries, which are not at equilibrium
  unb=unbalancedForce()
  print 'unbalanced force:',unb,' mean stress: ',triax.meanStress
  if unb<stabilityThreshold and abs(-10000-triax.meanStress)/10000<0.001:
    break
print 'timestep at confinedState:',O.iter

import sys #this is only for the flush() below
while triax.porosity>targetPorosity:
	## we decrease friction value and apply it to all the bodies and contacts
	compFricDegree = 0.95*compFricDegree
	setContactFriction(radians(compFricDegree))
	#print "\r Friction: ",compFricDegree," porosity:",triax.porosity,
	sys.stdout.flush()
	## while we run steps, triax will tend to grow particles as the packing
	## keeps shrinking as a consequence of decreasing friction. Consequently
	## porosity will decrease
	O.run(1000,1)
print 'porosity:',triax.porosity
print 'timestep at compactedState:',O.iter
print "###    Compacted state saved      ###"

##############################
###   DEVIATORIC LOADING   ###
##############################

##We move to deviatoric loading, let us turn internal compaction off to keep particles sizes constant
triax.internalCompaction=False

## Change contact friction (remember that decreasing it would generate instantaneous instabilities)
setContactFriction(radians(finalFricDegree))

##set stress control on x and z, we will impose strain rate on y
triax.stressMask = 5
##now goal2 is the target strain rate
triax.goal2=rate
## we define the lateral stresses during the test, here the same 10kPa as for the initial confinement.
triax.goal1=-10000
triax.goal3=-10000

newton.damping=0.1
O.engines[2].lawDispatcher.functors[1].always_use_moment_law = True
#We assign cohesion to all contacts at the next iteration
O.engines[2].physDispatcher.functors[1].setCohesionNow = True

O.run(1000,True)


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