Make Radius Ratio for the sample and verity of particles radii
Hello,
I made some changes to Bruno's Triaxial code as follows. The problem is, I want to have the ratio of radii as R_max/R_min=3, however, when I get the output of radii at the end of my code, it seems all particles have the same size as 0.6 mm (as the R_avr is even 0.6 mm as well). Is there anything wrong with my code regarding radius defining or I must use another way to make the sample?
Thank you for your help.
#######
######### TRIAXIAL PROBLEM, Y IS THE VERTICAL AXIS, X IS THE RIGHT AXIS, Z IS THE FRONT AXIS #########
#######
print ('************** START **************')
import numpy as np
import time
import datetime, os
import shutil
start_time=
from datetime import datetime
import math
from yade import qt, export, utils
from yade import pack
import yade.timing;
from pprint import pprint
O.timingEnabled
#######
######### DEFINING VARIABLES #########
print ('============ DEFINING VARIABLES ============')
nRead=readParam
num_spheres=5000,
compFricDegree = 10,
key='_triax_',
unknownOk=True
)
from yade.params import table
num_spheres=
key=table.key
targetPorosity = 0.404
compFricDegree = table.compFricD
finalFricDegree = 30
IP=2000 # iteration period to record data and stuff
IPmacro=100
micro_record_
ORN=1000 # O.Run Number of iterations
micro_record_
micro_record_
damp=0.2
thick=0
stabilityThresh
PCPC=0.0001 # Precision of Confining Pressure Convergence
Ls=0.02 # m
Lh=1*Ls
r_min=0.15e-3 # m
d_min=2*r_min # m
r_max=0.45e-3 # m
d_max=2*r_max # m
r_avr=(
d_avr=2*r_avr # m
mn,mx=Vector3(
volume = (mx[0]-
r_fuz=0.2*r_avr # m
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! ***")
young=2*(1e8) # Kn = 2 (E r1 E r2 / E r1 + E r2) >>> E = Kn/r_avr
print ('The Young modulus is:', young)
Kn=young*(r_avr) ### FIXME
print ('The Stiffness is:', Kn)
Kt=Kn ### FIXME
poisson=Kt/Kn # Kt/Kn
cvs=1-targetPor
mean_rad = pow(0.24*
print('
strainRate=-0.5 # %/sec ### FIXME
target_strain=0.4 ### FIXME %
sigmaIso=-5e5 # Pa=N/m2 ### FIXME
particleDensity
#######
################# 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.
s11 = -triax.
s22 = -triax.
s33 = -triax.
i = O.iter,
t = O.time, # virtual (yade) time --- time of simulation
fab = utils.fabricTen
#######
################# DEFINING MATERIALS #################
print ('============ DEFINING MATERIALS ============')
O.materials.
O.materials.
#######
################# DEFINING PACKING #################
print ('============ DEFINING PACKING ============')
walls=aabbWalls
for w in walls:w.
wallIds=
sp=pack.
clumps=False
sp.makeCloud(
O.bodies.
#######
################# DEFINING TRIAXIAL TEST #################
print ('============ DEFINING TRIAXIAL TEST ============')
triax=TriaxialS
maxMultiplier=
finalMaxMultip
thickness = thick,
stressMask = 7,
internalCompac
)
#######
################# DEFINING ENGINES #################
print ('============ DEFINING ENGINES ============')
newton=
O.engines=[
ForceResetter(),
#GravityEng
InsertionSortC
#InsertionS
InteractionLoop(
[Ig2_
[Ip2_
[Law2_
),
GlobalStiffnes
triax,
PyRunner(
PyRunner(
TriaxialStateR
newton
]
#######
################# 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=
ConfStressR
print ('~~~~~~~~~~~~~ Phase_01: Converging to Isotropic Compression, 50kPa ~~~~~~~~~~~~~')
print ('ConfSratio:
print ('mean stress engine:
if unb<stabilityTh
break
export.
e22Check=
print ('Axial Strain',e22Check)
print ('porosity:
print ('void ratio:'
print ('Mean stress engine:
print ('Mean stress (Calculated)
print ('#####
#######
################# REACHING TARGET POROSITY ##################
print ('============ REACHING TARGET POROSITY ============')
import sys
while triax.porosity>
O.run(ORN,True)
compFricDegree = 0.95*compFricDegree
setContactF
print ('\r Friction: ',compFricDegree,' porosity:
sys.
print ('#####
#######
################ PRINT SOME CHECK VARIABLES #################
RRmax=max(
RRmin=min(
print('
print('
print ('Number of elements:', len(O.bodies))
Vt = triax.boxVolume
print ('Vt:', Vt)
Vs = triax.particles
print('Vs:',Vs)
Vv=Vt-Vs
print('Vv:',Vv)
n = triax.porosity
print ('n:',triax.
e=n/(1-n)
print ('e:',e)
print ('Step that starts the deviatoric loading ', O.iter/ORN)
#######
################# DEVIATORIC LOADING #################
print ('============ APPLYING DEVIATORIC LOADING ============')
triax.internalC
setContactFrict
# Allow walls to move. >>> False means fixed , True means free to move
triax.wall_
triax.wall_
triax.wall_
triax.wall_
triax.wall_
triax.wall_
# Make Stress Control on X,Y,Z directions by setting True, False means Strain Control
triax.stressCon
triax.stressCon
triax.stressCon
# set the raite of strain to the desirable axis
#triax.
triax.strainRat
#triax.
# set the deviatoric conditions
triax.stressMask = 5
triax.goal1 = sigmaIso
triax.goal2 = strainRate
triax.goal3 = sigmaIso
newton.damping=0.1
while 1:
O.run(ORN,True)
unb=
axialS=
eps2=
eps2cal=
print ('~~~~~~~~~~~~~ Phase_02: Converging to Deviatoric Compression, Strain Rate ~~~~~~~~~~~~~')
print ('sigma2: ',axialS, ' q: ', axialS-sigmaIso,' step= ', O.iter/ORN,' Time:',O.time, ' iter/sec',O.speed)
print ('Strain2 Calculated: ',eps2cal, ' Strain2 engine: ',eps2,'axial deformation (%)', (eps2-e22Check)
if abs(eps2-
break
print ('#####
#######
#######
RRmax=max(
RRmin=min(
RRavr=sum(
print([
print('
print('
print('
print('
print ('Number of elements:', len(O.bodies))
Vt = triax.boxVolume
print ('Vt:', Vt)
Vs = triax.particles
print('Vs:',Vs)
Vv=Vt-Vs
print('Vv:',Vv)
n = triax.porosity
print ('n:',triax.
e=n/(1-n)
print ('e:',e)
#######
################# END #################
plt.show()
O.pause()
Question information
- Language:
- English Edit question
- Status:
- Solved
- For:
- Yade Edit question
- Assignee:
- No assignee Edit question
- Solved by:
- Jan Stránský
- Solved:
- 2021-02-16
- Last query:
- 2021-02-16
- Last reply:
- 2021-02-16
|
#1 |
Hello,
> r_min=0.15e-3 # m
> r_max=0.45e-3 # m
> r_avr=(
> r_fuz=0.2*r_avr # m
> sp.makeCloud(
r_fuz (rRelFuzz parameter) is **relative** to r_avr [1] (not m, but dimensionless), so put 0.2 (or similar number) instead of 0.2*r_avr (which makes it for "small" r_avr negligible)
cheers
Jan
[1] https:/
ehsan benabbas (ehsanben) said : | #2 |
Thanks Jan Stránský, that solved my question.