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?

[1] https://gitlab.com/yade-dev/trunk/blob/master/examples/triax-tutorial/script-session1.py

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
from yade import pack

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

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

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

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

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

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

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

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

######################################
################# 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')
yade.timing.stats()

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

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

Thank you for your helps.

Question information

Language:
English Edit question
Status:
Answered
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
2020-02-11
Last reply:
2020-02-11
ehsan benabbas (ehsanben) said : #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 : #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 : #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
from yade import pack, plot,qt
fIDSI=O.bodies.append(utils.geom.facetBox((.15,.15,.135),(.15,.15,.045),wallMask=15,material=steel))
fIDSII=O.bodies.append(utils.geom.facetBox((.15,.15,.045),(.15,.15,.045),wallMask=31,material=steel))
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]

[1] https://answers.launchpad.net/yade/+question/210746

Jan Stránský (honzik) said : #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

Can you help with this problem?

Provide an answer of your own, or ask ehsan benabbas for more information if necessary.

To post a message you must log in.