# How to make the boundaries flexible?

Asked by ehsan benabbas on 2020-02-10

Hello all,

I am using Ubuntu 18.04, and Yade 2019-08-08.git-775ae74

I developed a Triaxial test based on [1] in the following form. This code is working with rigid walls/boundaries which move in the "deviatoric" part of the test and this leads to the deformation of the specimen. In this case, I will not be able to capture the shear band and localization. To do so, I need to define the walls/boundaries as the deformable/flexible boundaries. In fact, I should not let the boundaries to move, but I should let them be deformed based on their flexibility. My question is that what do I need to do to make these rigid walls to flexible boundaries?

My code:

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

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

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

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 #################

O.materials.append(FrictMat(young=young,poisson=poisson,frictionAngle=0,density=0,label='walls'))

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

walls=aabbWalls([mn,mx],thickness=thick,material='walls')
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)

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

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

triax=TriaxialStressController(
maxMultiplier=1.+2e4/young,
finalMaxMultiplier=1.+2e3/young,
thickness = thick,
internalCompaction=True,
)

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

def history():
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 #################

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

##########################################################
################# 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)

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

print ('============ REACHING TARGET POROSITY ============')
import sys
while triax.porosity>targetPorosity:
compFricDegree = 0.95*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)

##################################################

triax.internalCompaction=False

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.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

######################################
################# 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()

%%%%%%%%%%%%%%%%%%%%%

## Question information

Language:
English Edit question
Status:
For:
Assignee:
No assignee Edit question
Last query:
2020-02-11
2020-02-11
 ehsan benabbas (ehsanben) said on 2020-02-10: #1

Sorry just to not get any error in this code, please remove the line"record_micro_data(str(O.iter/ORN))" in the deviatoric part.

 Jan Stránský (honzik) said on 2020-02-11: #2

Hi,

I think it is not possible directly in Yade at the moment.. and also in general it might be a tricky task (e.g. the mass of the flexible boundary is low, not very good for explicit dynamic simulations, the formulation how the membrane should behave etc.).

One option is to discretize the boundaries with facets. Based on the resultant forces on facets, you can change their vertices. How? It would be difficult answer :-)

cheers
Jan

 ehsan benabbas (ehsanben) said on 2020-02-11: #3

As the facet elements are flexible, isn't it possible to define the box (walls) with facets? Something like:

# create box with free top, and ceate loose packing inside the box
sp=pack.SpherePack()
sp.makeCloud((0,0,0),(0.3,0.3,0.3250),rMean=rMean,rRelFuzz=rRelFuzz)
sp.toSimulation(material=spheremat)

I have read this is [1]

 Jan Stránský (honzik) said on 2020-02-11: #4

> isn't it possible to define the box (walls) with facets?

sure, this is easy

> As the facet elements are flexible

no, facet in Yade is a rigid triangle.
But you can change its position or set new coordinates of vertices arbitrarily (e.g. according to some computation how a flexible membrane should deform)

cheers
Jan