yade-dev team mailing list archive

[luc scholtes <lscholtes63@xxxxxxxxx>]

Hi guys,

After few investigations about CohesiveFrictionalPM in latest version, I
must admit that I still do not have any idea where to look for, so if one of
you could have any suggestion, it would be very kind.

1- CohesiveFrictionalPM seems to work well when it is used in basic
simulations (say two elements, sphere/sphere or sphere/box) and without the
moment law (the moment law produces the same "explosive behaviour" when used
with CohesiveFrictionalContactLaw, so it is probably due to some problems
with quaternion computation in Eigen as it worked well with wm3... Unless
there was something wrong with wm3...)

2- I don't understand why, but when CohesiveFrictionalPM is used in a
TriaxialTest or a uniaxial test (without the moment law, though), everything
explodes. Since there have been a lot of changes during the last weeks (I
use rev2185 without any problem), I don't really know where to check
first... I have attached a script that I use for Triaxial Testing. I tested
it with CohesiveFrictionalPM as well as with CohesiveFrictionalContactLaw.
RESULTS -> if the moment law is desactivated, it works for
CohesiveFrictionalContactLaw but not for CohesiveFrictionalPM.

The first thing I have noticed is that, if the timestep is fixed with
O.dt=-0.2*utils.PWaveTimeStep() (see the - before
0.2*utils.PWaveTimeStep()), GlobalStiffnessTimeStepper() does not update the
timestep for CohesiveFrictionalPM (it does for
CohesiveFrictionalContactLaw), and, the walls just go through the assembly.

Then, fixing the timestep with O.dt=0.2*utils.PWaveTimeStep() (without the
-), I could see that after few iterations, everything explodes again. Having
a look into the TriaxialRecorder file, I could see that the kinematic
energy, as well as the unbalanced force increase dramatically.

The only (!!!!) difference between the two laws is that, one uses the
penetration depth as the interacting distance, whereas the other initializes
the interacting distance to the penetration depth when an interaction is
detected, and I saw this: https://bugs.launchpad.net/yade/+bug/399963. Do
you think it could be related?

Please... help

  Luc

# -*- coding: utf-8 -*-
# encoding: utf-8
from yade import ymport, export, utils

"""
A basic script to allow triaxial testing of a given pre existing assembly

"""

##### definition of the material 
## CohesiveFrictionalContactLaw
#def sphereMat(): return CohFrictMat(young=115e8,poisson=.2,frictionAngle=radians(40),density=1400)
#wallMat=O.materials.append(CohFrictMat(young=115e8,poisson=.2,frictionAngle=radians(0),density=1400))

## CohesiveFrictionalPM
def sphereMat(): return CFpmMat(type=1,young=115e8,poisson=.2,frictionAngle=radians(40),density=1400)
wallMat=O.materials.append(CFpmMat(type=0,young=115e8,poisson=.2,frictionAngle=radians(0),density=1400))

##### insert spheres from existing packing into the scene
O.bodies.append(ymport.text('400.spheres',scale=1.,material=sphereMat))

##### create walls around the packing
walls=utils.aabbWalls(oversizeFactor=1.5,material=wallMat)
wallIds=O.bodies.append(walls)

##### triaxial compressionEngine parameters
triax=TriaxialCompressionEngine(
	wall_bottom_id=wallIds[2],
	wall_top_id=wallIds[3],
	wall_left_id=wallIds[0],
	wall_right_id=wallIds[1],
	wall_back_id=wallIds[4],
	wall_front_id=wallIds[5],
	StabilityCriterion=0.01,
	translationAxis=(0,1,0),
	internalCompaction=False,
	sigmaIsoCompaction=2e5,
	autoCompressionActivation=True,
	sigmaLateralConfinement=2e5,

	wallDamping=0.,
	strainRate=0.1,
	max_vel=0.1,
	epsilonMax=0.015,
	label='compressionEngine')

##### define the interaction radius for the first time step of the simulation (create bonds)
interactionRadius=-1;

##### engines definition
O.engines=[
	ForceResetter(),
	BoundDispatcher([Bo1_Sphere_Aabb(aabbEnlargeFactor=interactionRadius,label='is2aabb'),Bo1_Box_Aabb()]), 
	InsertionSortCollider(),

	#### CohesiveFrictionalContactLaw
	#InteractionDispatchers(
		#[Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=interactionRadius,label='ss2d3dg'),Ig2_Box_Sphere_ScGeom()],
		#[Ip2_2xCohFrictMat_CohFrictPhys(setCohesionNow=True,normalCohesion=115e5,shearCohesion=115e6)],
		#[Law2_ScGeom_CohFrictPhys_ElasticPlastic(momentRotationLaw=False,shear_creep=0,twist_creep=0)]
	#),

	#### CohesiveFrictionalPM
	InteractionDispatchers(
		[Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=interactionRadius,label='ss2d3dg'),Ig2_Box_Sphere_ScGeom()], 
		[Ip2_CFpmMat_CFpmMat_CFpmPhys(cohesiveTresholdIteration=1,useAlphaBeta=True,Alpha=0.2,Beta=0.,eta=0.,tensileStrength=0.,cohesion=0.,strengthSoftening=0.,label='interactionPhys')],
		[Law2_ScGeom_CFpmPhys_CohesiveFrictionalPM(label='interactionLaw')]
	),

	GlobalStiffnessTimeStepper(),
	triax,
	TriaxialStateRecorder(file='test_bzr',iterPeriod=1,label='recorder'),
	NewtonIntegrator(damping=.4)
	
]

##### timeStep definition
O.dt=-.2*utils.PWaveTimeStep() # initial timestep to avoid initial explosion -> sign (-) should enable GlobalStiffnessTimeStepper(), but does not work for CohesiveFrictionalPM!!!
print O.dt

##### opening YADE interface
from yade import qt
v=qt.Controller()
v=qt.View()

##### to manage interaction detection factor during the first timestep
O.step();
## initializes the interaction detection factor to the default value
ss2d3dg.interactionDetectionFactor=-1.
is2aabb.aabbEnlargeFactor=-1.

##### save initial simualtion
O.save('TriaxialTest.xml');

O.wait()

Attachment: 400.spheres
Description: Binary data