yade-users team mailing list archive

[ehsan benabbas <question688575@xxxxxxxxxxxxxxxxxxxxx>]

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

    Status: Answered => Open

ehsan benabbas is still having a problem:
Minimal scripts:


print ('************** START **************')
import numpy as np
import time
import datetime, os
start_time=time.time()
from datetime import datetime
import math
from yade import qt, export, utils
from yade import pack

##################################
######### DEFINING VARIABLES #########

print ('============ DEFINING VARIABLES ============')

nRead=readParamsFromTable(
 num_spheres=20000,
 compFricDegree = 29,
 key='_triax_',
 unknownOk=True
)

from yade.params import table

num_spheres=table.num_spheres
key=table.key
targetPorosity = 0.4
compFricDegree = table.compFricDegree
finalFricDegree = 29
IP=100 # iteration period to record data and stuff
micro_record_iterPeriod=IP
ORN=3000 # O.Run Number of iterations
micro_record_enable_normal_branch=True
micro_record_float_type = np.float32
damp=0.2
thick=0
stabilityThreshold=0.01
PCPC=0.0001 # Precision of Confining Pressure Convergence
r_min=0.1*1e-3 # m
d_min=2*r_min # m
r_max=0.3*1e-3 # m
d_max=2*r_max # m
r_avr=(r_min+r_max)/2 # m
d_avr=2*r_avr # m
r_fuz=(r_max/r_avr)-1 # m
Kn=10e8*(d_avr) ### FIXME
Kt=10e8*(d_avr) ### FIXME
young=Kn/r_avr # 2 (E r1 E r2 / E r1 + E r2) >>> E = Kn/r_avr
poisson=Kn/Kt # Kt/Kn
Ls=0.02 # m length of specimen ### FIXME
L_REV=7*(d_avr) # m
if Ls < L_REV:
    sys.exit("*** ERROR! The specimen's dimension is too samll! ***")
elif Ls==L_REV:
    print ("*** This is the minimum specimen's dimension you can take! ***")
else:
    print ("*** The specimen's dimension is good enough! ***")
mn,mx=Vector3(0,0,0),Vector3(Ls,Ls,Ls)
Vt=-1*1e-3 # m/s # negative sign describes the compression direction
strainRate=Vt/Ls # %/sec
target_strain=0.25 ### FIXME %
print ("The target strain has been set to:", target_strain)
sigmaIso=-5e5 # Pa ### FIXME
particleDensity=2000 #kg/m3

##################################################
################# DEFINING MATERIALS #################

print ('============ DEFINING MATERIALS ============')
O.materials.append(FrictMat(young=young,poisson=poisson,frictionAngle=radians(compFricDegree),density=particleDensity,label='spheres'))
O.materials.append(FrictMat(young=young,poisson=poisson,frictionAngle=0,density=0,label='walls'))

################################################
################# DEFINING PACKING #################

print ('============ DEFINING PACKING ============')
walls=aabbWalls([mn,mx],thickness=thick,material='walls')
for w in walls:w.shape.radius=0
wallIds=O.bodies.append(walls)
sp=pack.SpherePack()
clumps=False
sp.makeCloud(mn,mx,r_avr,r_fuz,num_spheres,False, 0.95,seed=1)
O.bodies.append([sphere(center,rad,material='spheres') for center,rad in sp])

from yade import export
os.mkdir('3axresults')
os.chdir('/home/ehsan/Desktop/3axresults')
export.text('InitialPackingData')

####################################################
################# DEFINING TRIAXIAL TEST #################

print ('============ DEFINING TRIAXIAL TEST ============')
triax=TriaxialStressController(
 maxMultiplier=1.+2e4/young,
 finalMaxMultiplier=1.+2e3/young,
 thickness = thick,
 stressMask = 7,
 internalCompaction=True,
)

##################################################
################# DEFINING FUNCTIONS #################

print ('============ DEFINING FUNCTIONS ============')
from yade import plot
def history():
    plot.addData(
  e11 = -triax.strain[0],
  e22 = -triax.strain[1],
  e33 = -triax.strain[2],
  ev = -triax.strain[0]-triax.strain[1]-triax.strain[2],
  s11 = -triax.stress(triax.wall_right_id)[0],
  s22 = -triax.stress(triax.wall_top_id)[1],
  s33 = -triax.stress(triax.wall_front_id)[2],
  i = O.iter,
        t = O.time, # virtual (yade) time --- time of simulation
        fab = utils.fabricTensor()[0])

################################################
################# DEFINING ENGINES #################

print ('============ DEFINING ENGINES ============')
newton=NewtonIntegrator(damping=damp)

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,
    PyRunner(iterPeriod=IP,command='history()',label='macro_recorder'),
 TriaxialStateRecorder(iterPeriod=IP,file='WallStresses'+table.key),
 newton
]

Gl1_Sphere.stripes=True
if nRead==0: yade.qt.Controller(), yade.qt.View()

##########################################################
################# APPLYING CONFINING PRESSURE #################

print ('============ APPLYING CONFINING PRESSURE ============')

triax.goal1 = sigmaIso
triax.goal2 = sigmaIso
triax.goal3 = sigmaIso

while 1:
    O.run(ORN,True)
    unb = unbalancedForce()
    meanS=(triax.stress(triax.wall_right_id)[0]+triax.stress(triax.wall_top_id)[1]+triax.stress(triax.wall_front_id)[2])/3
    ConfStressRatio=abs(sigmaIso-triax.meanStress)/abs(sigmaIso)
    print ('~~~~~~~~~~~~~ Phase_01: Converging to Isotropic Compression, 50kPa ~~~~~~~~~~~~~')
    print ('mean stress engine:',triax.meanStress,' mean stress (Calculated):',meanS, ' ConfSratio:',ConfStressRatio,' step:', O.iter/ORN, ' Time:',O.time, ' TimeStep',O.dt)
    print ('porosity:',triax.porosity, ' void ratio:',triax.porosity/(1-triax.porosity))
    if unb<stabilityThreshold and ConfStressRatio<PCPC:
        break

export.text('FinalPhase01PackingData')
e22Check=-triax.strain[1] ###%%%***
print ('Axial Strain',e22Check)
print ('Mean stress engine: ',triax.meanStress)
print ('Mean stress (Calculated): ',meanS)
print ('################## Isotropic phase is finished and saved successfully ##################')

########################################################
################# REACHING TARGET POROSITY ##################

print ('============ REACHING TARGET POROSITY ============')
import sys
while triax.porosity>targetPorosity:
 compFricDegree = 0.95*compFricDegree
 setContactFriction(radians(compFricDegree))
 print ('\r Friction: ',compFricDegree,' porosity:',triax.porosity,'step= ',O.iter/ORN,' Time:',O.time, ' TimeStep',O.dt)
 sys.stdout.flush()
 O.run(ORN,True)

print ('################## Target porosity is reached and compacted
state saved successfully ##################')

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

print ('============ APPLYING DEVIATORIC LOADING ============')
triax.internalCompaction=False
setContactFriction(radians(finalFricDegree))

triax.wall_top_activated=True
triax.wall_bottom_activated=False
triax.wall_left_activated=True
triax.wall_right_activated=True
triax.wall_back_activated=True
triax.wall_front_activated=True

triax.stressControl_1=True
triax.stressControl_2=False
triax.stressControl_3=True

triax.strainRate2=strainRate

triax.stressMask = 5
triax.goal1 = sigmaIso
triax.goal2 = strainRate
triax.goal3 = sigmaIso

newton.damping=0.1 ### FIXME
os.mkdir('Phase02PackingData')
os.chdir('/home/ehsan/Desktop/3axresults/Phase02PackingData')

while 1:
    O.run(ORN,True)
    export.text(str(O.iter/ORN))
    record_micro_data(str(O.iter/ORN))
    unb=unbalancedForce()
    axialS=triax.stress(triax.wall_top_id)[1]
    eps2=triax.strain[1]
    eps2cal=(triax.height-Ls)/Ls
    print ('~~~~~~~~~~~~~ Phase_02: Converging to Deviatoric Compression, Strain Rate ~~~~~~~~~~~~~')
    print ('sigma2: ',axialS, ' q: ', axialS-sigmaIso,' step= ', O.iter/ORN,' Time:',O.time, ' TimeStep',O.dt)
    print ('Strain2 Calculated: ',eps2cal, ' Strain2 engine: ',eps2,'axial deformation (%)', (eps2-e22Check)*100)
    print('A:', abs(eps2-e22Check), ' B:',target_strain)
    if abs(eps2-e22Check)>=target_strain: #triax.sigma2
        break

print ('################## Deviatoric phase is finished and saved
successfully ##################')

######################################
################# END #################

print ('************** END **************')
O.realtime
print ('Analysis has been taken for',O.realtime, 'seconds or', O.realtime/60, 'minutes')
print('Real time of run:',(time.time() - start_time), 'seconds or',(time.time() - start_time)/60, 'minutes')

#######################################################
################# RECORD Macro DATA and Plot #################

print ('============ RECORD AND PLOT DATA ============')

O.run(ORN,True)

plot.plots={'e22':('ev',),'j':('s22',),'e22 ':('s22',),' j':('s11','s33',)}
plot.labels={'s11':'$\sigma_{11}$' , 's22':'$\sigma_{22}$' , 's33':'$\sigma_{33}$' , 'e22':'$\epsilon_{22}$' , 'ev':'$\epsilon_{V}$' , 'i':'Time Step'}
 plot.plot()

os.chdir('/home/ehsan/Desktop/3axresults')
plot.saveDataTxt('Macro_results')

O.pause()

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