yade-users team mailing list archive

[Saeed <question681991@xxxxxxxxxxxxxxxxxxxxx>]

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

Hi There,
When I run the code below, I will be given the Error when arrives to > unsat.initialization()
I guessed the grains out of the box, which I asked in the code to remove that grain, but the grains were not placed outside the box.
Also, using the lower timestep but the problem remains.

Code :

## encoding: utf-8
import matplotlib; matplotlib.rc('axes',grid=True)
import matplotlib.pyplot as plt
import pylab
from yade import pack
import numpy as np
from numpy import *
import math

## Define Parameters
num_spheres=5000			
compFricDegree=30	
confiningS=-12500
RhoW=1000	
RhoS=2700
youngModulus=75e6
poissonRatio=0.3

## creat a packing with a specific particle size distribution (PSD)
psdSizes=[0.0016,0.002,0.005,0.006,0.0075,0.01,0.015,0.025,0.035,0.05,0.075] #(mm)
psdCumm=[0.0,0.075,0.21,0.23,0.28,0.37,0.48,0.70,0.85,0.89,1.0] #cumulative
psdSizesArray=np.array(psdSizes)
psdSizesMeter=psdSizesArray*0.001 #Convert the size of particles to meter
sp=pack.SpherePack()
mn,mx=Vector3(0,0,0),Vector3(.0012,0.0012,0.0012) #initial box = 1.2*1.2*1.2 mm
sp.makeCloud(minCorner=mn,maxCorner=mx,num=num_spheres,psdSizes=psdSizesMeter,psdCumm=psdCumm,distributeMass=False,seed=True)
sp.psd(bins=150,mass=False)
sp.toSimulation

## create material, which will be used as default
O.materials.append(FrictMat(young=youngModulus,poisson=poissonRatio,frictionAngle=radians(compFricDegree),density=RhoS,label='spheres'))
O.materials.append(FrictMat(young=youngModulus,poisson=poissonRatio,frictionAngle=0,density=0,label='frictionless'))

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

O.bodies.append([utils.sphere(center,rad,material='spheres') for center,rad in sp])

triax=TriaxialStressController(
	internalCompaction=False,	#external (walls) compaction , turn internal compaction off to keep particles sizes constant
	goal1=confiningS,	
	goal2=confiningS,		
	goal3=confiningS,		
	label="triax"
)

newton=NewtonIntegrator(damping=0.4) 

## Engines
O.engines=[
	ForceResetter(),
	InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
	InteractionLoop(
		[Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
		[Ip2_FrictMat_FrictMat_FrictPhys()],
		[Law2_ScGeom_FrictPhys_CundallStrack()]
	),
	GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
	triax,
	newton
]

O.dt=.1*PWaveTimeStep()
O.dynDt=False

while 1:
  O.run(1000,True)
  unb=unbalancedForce()
  e=(triax.porosity)/(1-triax.porosity)
  print 'unb: ',unb,' e: ',e,' Sigma: ',-triax.stress(3)[1]
  if unb<0.15 and abs(confiningS-triax.meanStress)/(-confiningS)<0.015 and e<=1.05:
    break

triax.depth0=triax.depth
triax.height0=triax.height
triax.width0=triax.width
O.save('Pack60.yade') #save the packing, which can be reloaded later.

## Plot and Save data
from yade import plot
O.engines=O.engines+[PyRunner(iterPeriod=100,command='history()',dead=1,label='recorder')]
X=0
def history():
  	plot.addData(e11=-triax.strain[0], e22=-triax.strain[1], e33=-triax.strain[2],
		    s11=-triax.stress(0)[0],s22down=-triax.stress(2)[1],s22up=-triax.stress(3)[1],s33=-triax.stress(4)[2],
		    sEfficient=abs(triax.stress(3)[1])+X*(unsat.bndCondValue[3]-unsat.bndCondValue[2]),
		    pc=unsat.bndCondValue[3]-unsat.bndCondValue[2],
		    Sr=unsat.getSaturation(False),
		    n=triax.porosity,
		    Vs=triax.particlesVolume,
		    V=triax.boxVolume,
		    Vv=triax.porosity*triax.boxVolume,	
		    Vw=unsat.getSaturation(False)*triax.porosity*triax.boxVolume,
		    Mw=RhoW*unsat.getSaturation(False)*triax.porosity*triax.boxVolume,
		    Ms=RhoS*triax.particlesVolume,
		    WaterContent=((RhoW*unsat.getSaturation(False)*triax.porosity*triax.boxVolume)/(RhoS*triax.particlesVolume))*100,
		    eV=-triax.volumetricStrain,		
		    gPP=unsat.getPorePressure((0,0,0)),
		    up=unsat.porosity,
		    VoidRatio=(triax.porosity)/(1-triax.porosity),
		    meanS=triax.meanStress,
		    xx=triax.width,
		    yy=triax.height,
		    zz=triax.depth,				    
		    i=O.iter
		    )

##########################
##     Start test      ###
##########################
triax.stressMask=2
triax.goal1=triax.goal3=0
goalTop=triax.stress(3)[1]
triax.goal2=goalTop
triax.wall_bottom_activated=False

## Instantiate a two-phase engine
unsat=TwoPhaseFlowEngine()
meanDiameter=(O.bodies[-1].shape.radius + O.bodies[6].shape.radius)/2.0

## set boundary conditions, the drainage is controlled by decreasing Water pressure and keeping air pressure constant
unsat.bndCondIsPressure=[0,0,1,1,0,0]
unsat.bndCondValue=[0,0,0,0,0,0]
unsat.isPhaseTrapped=False 
unsat.defTolerance=0.05 		
unsat.meshUpdateInterval=200 
unsat.useSolver=3 	
unsat.permeabilityFactor=1 	
unsat.updateTriangulation=True 		
unsat.surfaceTension = 0.0728 
unsat.initialization()			

## Drainage
for pg in arange(1.0e-5,15.0,0.125):
  unsat.bndCondValue=[0,0,(-1.0)*pg*unsat.surfaceTension/meanDiameter,0,0,0]
  unsat.invasion()
  unsat.computeCapillaryForce()
  for b in O.bodies:
    O.forces.setPermF(b.id, unsat.fluidForce(b.id))
  while 1:
    O.run(1000,True)
    unb=unbalancedForce()
    Moisture=((RhoW*unsat.getSaturation(False)*triax.porosity*triax.boxVolume)/(RhoS*triax.particlesVolume))*100
    print 'Drainage unb: ',unb,' e: ',(triax.porosity)/(1-triax.porosity),' w: ',Moisture,' PC: ',-unsat.bndCondValue[2]
    if unb<0.5:
      break
  if Moisture<18.5:
      break
--------------------------------------------------------------------.
error:
WARNING! rmin>rmax. rmin=2.10136e-05 ,rmax=1e-10
WARNING! rmin>rmax. rmin=1.75633e-05 ,rmax=1e-10
WARNING! rmin>rmax. rmin=2.02999e-05 ,rmax=1e-10
WARNING! rmin>rmax. rmin=1.50292e-05 ,rmax=1e-10
WARNING! rmin>rmax. rmin=2.16551e-05 ,rmax=1e-10
WARNING! rmin>rmax. rmin=1.67435e-05 ,rmax=1e-10
WARNING! rmin>rmax. rmin=2.45035e-05 ,rmax=1e-10
WARNING! rmin>rmax. rmin=2.10136e-05 ,rmax=1e-10
WARNING! rmin>rmax. rmin=2.27934e-05 ,rmax=1e-10
WARNING! rmin>rmax. rmin=2.02674e-05 ,rmax=1e-10
WARNING! rmin>rmax. rmin=1.51404e-05 ,rmax=1e-10

-------------------------------------------------------------------------


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