Tips for life, tips for tech

Well this is kinda hard problem. Here is the draft code for the process in Yade. Remember for this version of script, the deformation in the Rendulic plan is not true (it must be lower, since the value recorded is started from the confining state, not the true deviatoric loading)

# -*- coding: utf-8 -*-<br />from yade import pack<br />#from utils import *<br /><br /><br />nRead=utils.readParamsFromTable(<br />    num_spheres=10000,# number of spheres<br />    compFricDegree = 15, # contact friction during the confining phase the final porosity is slightly bigger than 0.409, it means we must lower this slightly, maybe 18.5<br />    unknownOk=True,<br />    isoForce=100000,<br />    conStress=100000<br />)<br />from yade.params import table<br /><br />num_spheres=table.num_spheres<br />## corners of the initial packing<br />mn,mx=Vector3(-0.24,-0.24,-0.24),Vector3(0.24,0.24,0.24)<br />thick = 0.01<br />compFricDegree = table.compFricDegree<br />finalFricDegree=35<br />targetPorosity=0.382<br />rate=0.1<br />damp=0.06<br />stabilityThreshold=0.001<br />compForce=100000<br />stickLoad=140000<br />key='_probing_envelope_e618_'<br /><br />## create material #0, which will be used as default<br />O.materials.append(FrictMat(young=356e6,poisson=0.42,frictionAngle=radians(compFricDegree),density=3000,label='spheres'))<br />O.materials.append(FrictMat(young=18e6,poisson=0.5,frictionAngle=0,density=0,label='walls'))<br /><br />## create walls around the packing<br />walls=utils.aabbWalls([mn,mx],thickness=thick,material='walls')<br />wallIds=O.bodies.append(walls)<br /><br />sp=pack.SpherePack()<br />#psdSizes=[0.002,0.003,0.004,0.005,0.006,0.007,0.008,0.0095] # (sizes or radii of the grains vary from 2mm to 9.5mm) <br />#psdCumm=[0.01,0.09,0.25,0.50,0.69,0.90,0.95,1] # the correspondent amount (percentage) of each diameter<br />psdSizes,psdCumm=[0.001978,0.00218,0.00249,0.002744,0.002894,0.002999,0.003162,0.003305,0.003455,0.00358,0.003709,0.003911,0.004016,0.004273,0.004427,0.004669,0.004968,0.005193,0.005476,0.005826,0.006036,0.00631,0.006595,0.006833,0.007079,0.007335,0.007532,0.007943,0.008526,0.0092],[0.001,0.002,0.005,0.026385,0.047493,0.068602,0.092348,0.12372,0.150396,0.174142,0.200528,0.242744,0.261214,0.311346,0.356201,0.401055,0.469657,0.522427,0.583113,0.654354,0.71504,0.770449,0.823219,0.878628,0.926121,0.968338,0.98153,0.992084,0.997361,1]<br />#---------------------------------------------<br />sp.makeCloud(mn,mx,-1,0,num_spheres,False, 0.95,psdSizes,psdCumm,False,seed=1)<br />#sp.makeCloud(mn,mx,-1,0,-1,False, 0.95,psdSizes,psdCumm,False,seed=1)<br />sp.toSimulation(material='spheres')<br /><br />volume = (mx[0]-mn[0])*(mx[1]-mn[1])*(mx[2]-mn[2])<br />mean_rad = pow(0.09*volume/num_spheres,0.3333)<br /><br />#clumps=False<br />#if clumps:<br />#    c1=pack.SpherePack([((-0.2*mean_rad,0,0),0.5*mean_rad),((0.2*mean_rad,0,0),0.5*mean_rad)])<br />#    sp.makeClumpCloud((-0.24,0.24,-0.24),(0.24,0.24,0.24),[c1],periodic=False)<br />#    O.bodies.append([utils.sphere(center,rad,material='spheres') for center,rad in sp])<br />#    standalone,clumps=sp.getClumps()<br />#    for clump in clumps:<br />#        O.bodies.clump(clump)<br />#        for i in clump[1:]: O.bodies[i].shape.color=O.bodies[clump[0]].shape.color<br />#    #sp.toSimulation()<br />#else:<br />#    O.bodies.append([utils.sphere(center,rad,material='spheres') for center,rad in sp])<br /><br />O.dt=.5*utils.PWaveTimeStep() # initial timestep, to not explode right away<br />O.usesTimeStepper=True<br /><br />triax=TriaxialStressController(<br />    ## ThreeDTriaxialEngine will be used to control stress and strain. It controls particles size and plates positions.<br />    ## this control of boundary conditions was used for instance in http://dx.doi.org/10.1016/j.ijengsci.2008.07.002<br />    maxMultiplier=1.001, # spheres growing factor (fast growth)<br />    finalMaxMultiplier=1.00001, # spheres growing factor (slow growth)<br />    thickness = thick,<br />    radiusControlInterval=20,<br />    ## switch stress/strain control using a bitmask. What is a bitmask, huh?!<br />    ## Say x=1 if stess is controlled on x, else x=0. Same for for y and z, which are 1 or 0.<br />    ## Then an integer uniquely defining the combination of all these tests is: mask = x*1 + y*2 + z*4<br />    ## to put it differently, the mask is the integer whose binary representation is xyz, i.e.<br />    ## "100" (1) means "x", "110" (3) means "x and y", "111" (7) means "x and y and z", etc.<br />    stressMask = 7,<br />    #the value of confining stress for the intitial (growth) phase<br />    goal1=compForce,<br />    goal2=compForce,<br />    goal3=compForce,<br />    max_vel=5,<br />    internalCompaction=False, # If true the confining pressure is generated by growing particles<br />    #Key=key # passed to the engine so that the output file will have the correct name<br />)<br /><br />newton=NewtonIntegrator(damping=damp)<br /><br />O.engines=[<br />    ForceResetter(),<br />    InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()],verletDist=-mean_rad*0.06),<br />    InteractionLoop(<br />        [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],<br />        [Ip2_FrictMat_FrictMat_FrictPhys()],<br />        [Law2_ScGeom_FrictPhys_CundallStrack()]<br />    ),<br />    GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8, defaultDt=4*utils.PWaveTimeStep()),<br />    triax,<br />    TriaxialStateRecorder(iterPeriod=100,file='WallStresses'+key),<br />    newton<br />]<br /><br />O.save('initial'+key+'.xml')<br />#Display spheres with 2 colors for seeing rotations better<br />#Gl1_Sphere.stripes=0<br />#if nRead==0: yade.qt.Controller(), yade.qt.View()<br />print 'Number of elements: ', len(O.bodies)<br />print 'Box Volume: ',  triax.boxVolume<br /><br />#Display spheres with 2 colors for seeing rotations better<br />Gl1_Sphere.stripes=0<br />yade.qt.Controller(), yade.qt.View()<br /><br /><br /># Applying confining force<br />triax.goal1=table.isoForce<br />triax.goal2=table.isoForce<br />triax.goal3=table.isoForce<br /><br />while 1:<br />  O.run(1000, True)<br />  #the global unbalanced force on dynamic bodies, thus excluding boundaries, which are not at equilibrium<br />  unb=unbalancedForce()<br />  #average stress<br />  #note: triax.stress(k) returns a stress vector, so we need to keep only the normal component<br />  meanS=(triax.stress(triax.wall_right_id)[0]+triax.stress(triax.wall_top_id)[1]+triax.stress(triax.wall_front_id)[2])/3<br />  print 'unbalanced force:',unb,' mean stress: ',meanS<br />  if unb&lt;stabilityThreshold and abs(meanS-table.isoForce)/table.isoForce&lt;0.001:<br />    break<br /><br />O.save('confinedState'+key+'.xml')<br />print "###      Isotropic state saved      ###"<br />print 'current porosity=',triax.porosity<br />print 'current void ratio=',triax.porosity/(1-triax.porosity)<br /><br />###### porosity reaching #####<br />#import sys #this is only for the flush() below<br />#while triax.porosity&gt;targetPorosity:<br />#    # we decrease friction value and apply it to all the bodies and contacts<br />#    compFricDegree = 0.95*compFricDegree<br />#    setContactFriction(radians(compFricDegree))<br />#    print "\r Friction: ",compFricDegree," porosity:",triax.porosity,<br />#    sys.stdout.flush()<br />#    # while we run steps, triax will tend to grow particles as the packing<br />#    # keeps shrinking as a consequence of decreasing friction. Consequently<br />#    # porosity will decrease<br />#    O.run(500,1)<br /><br />#O.save('compactedState'+key+'.yade.gz')<br />#print "###    Compacted state saved      ###"<br /><br />print 'current porosity=',triax.porosity<br />print 'current void ratio=',triax.porosity/(1-triax.porosity)<br />print 'step that starts the deviatoric loading ', O.iter<br /><br />from yade import plot<br />from pprint import pprint<br />plot.reset()<br />plot.addData(CheckpointStep=O.iter,Porosity=triax.porosity,e22=triax.strain[1])<br />plot.saveDataTxt('checkpoint.txt',vars=('CheckpointStep','Porosity','e22'))<br />##############################<br />###   DEVIATORIC LOADING   ###<br />##############################<br /><br />#We move to deviatoric loading, let us turn internal compaction off to keep particles sizes constant<br />triax.internalCompaction=False<br /><br /># Change contact friction (remember that decreasing it would generate instantaneous instabilities)<br />setContactFriction(radians(35))<br /><br />#set stress control on x and z, we will impose strain rate on y<br />triax.stressMask = 5<br />#now goal2 is the target strain rate<br />triax.goal2=-rate<br />##we assign constant stress to the other directions<br />triax.goal1=100000<br />triax.goal3=100000<br /><br />##we can change damping here. What is the effect in your opinion?<br />newton.damping=0.1<br /><br />#Save temporary state in live memory. This state will be reloaded from the interface with the "reload" button.<br />O.saveTmp()<br /><br /><br />#First perform a transverse isotropic compression (or a stress controlled drained triaxial compression) to reach an initial state from where stress probes will be applied<br />#... need to active stress control in 3 directions<br />triax.stressMask = 7<br />#... choose the value of axial stress where we want to stop the compression<br />triax.goal2=stickLoad<br />#... fix a maximum strain rate to go progressivly to the desired stress state in direction 2<br />triax.stressMask=0<br />triax.strainRate2=0.01<br />#... fix a high value of maximum strain rate in radial direction to be sure to keep in any conditions a constant confining pressure<br />triax.strainRate1=triax.strainRate3=1000.0<br /><br /><br />while 1:<br />  triax.stressMask=7<br />  O.run(100, True)<br />  #the global unbalanced force on dynamic bodies, thus excluding boundaries, which are not at equilibrium<br />  unb=unbalancedForce()<br />  #note: triax.stress(k) returns a stress vector, so we need to keep only the normal component<br />  axialS=triax.stress(triax.wall_top_id)[1]<br />  print 'unbalanced force:',unb,' sigma2: ',axialS<br />  if unb&lt;stabilityThreshold and abs(axialS-triax.goal2)/triax.goal2&lt;0.001: #triax.sigma2<br />    break<br /><br />O.save('anisotropicState'+key+'.xml')<br /><br /><br />#Perform stress probes from the anisotropic state<br /><br />dSnorm = 100.0 #norm of the stress increment<br />nbProbes = 36 #number of stress directions tested<br />rampIte = 20 #nb iterations to increase the stress state until the final desired stress value<br />#LUC: je fixe des nombres d'iterations c'est moins elegant que de chercher explicitement un etat d'equilibre mais ca permet de poursuivre le calcul meme si un etat de contrainte n'est pas correctement atteint pour un stress probe et qu'il est difficile de se stabiliser a cet etat de contrainte (i.e. attendre longtemps...)<br />stabIte = 5000 #nb iterations to stabilize sample after reaching the final stress value<br /><br /><br /># an array for saving stress increments and strain responses; arrays are in "numpy" extension<br />import numpy<br />probings=numpy.zeros((3,nbProbes))<br /><br />def increment(dsr=0,dsa=1):<br />    for ite in range(rampIte):# progressivaly increase of stress state<br />        O.run(20, True)<br />        #incrementation of stress state<br />        triax.goal2 = initSa+dsa/rampIte*ite<br />        triax.goal1 = triax.goal3 = initSr+dsr/rampIte*ite<br />        print triax.goal1, triax.goal2<br /><br />    # fix the stress value for stabilization at the final state<br />    triax.goal2 = initSa+dsa<br />    triax.goal1 = triax.goal3 = initSr+dsr<br /><br />    while 1:<br />        O.run(100, True)<br />        unb=unbalancedForce()<br />        print 'unbalanced force:',unb,' strain: ',triax.strain<br />        if unb&lt;stabilityThreshold: break<br /><br /><br /><br /># loop over all the stress directions<br />for i in range(nbProbes):<br /><br />    # computation of the stress direction of the current stress probe<br />    alphaS = 2*pi/nbProbes*(i-1)<br />    print 'stress probe nb:',i,' stress direction (deg): ',degrees(alphaS)<br /><br />    # computation of the stress increment in the axial direction<br />    dSa = dSnorm*sin(alphaS) <br /><br />    # computation of the stress increment in the radial direction<br />    dSr = dSnorm*cos(alphaS)/sqrt(2.0)<br /><br />    #Load the initial anisotropic state before running a new stress probe<br />    O.load('anisotropicState'+key+'.xml')<br />    #We redefine the "triax" label, else it would point to inactive engine from previous simulation that is still in memory<br />    triax=O.engines[4]<br /><br />    initSa=triax.goal2  #save of the initial axial stress<br />    initSr=triax.goal1  #save of the initial radial stress<br /><br />    # define the final stress state to be reached<br />    finalSa = initSa+dSa<br />    finalSr = initSr+dSr<br /><br />    #... need to active stress control in 3 directions if not yet done<br />#    triax.stressControl_1=triax.stressControl_2=triax.stressControl_3=True<br />    triax.stressMask=7<br /><br />    # fix a high value of maximum strain rate, the progressive loading will be done by progressively increasing the desired stress state at each iteration<br />#    triax.strainRate1=triax.strainRate2=triax.strainRate3=1000.0<br />    triax.goal1=triax.goal2=triax.goal3=1000.0<br />    increment(dSr,dSa)<br />    probings[:,i]=triax.strain<br />    <br /><br /><br />#open a file for writing probing results<br />a=open('probings'+key,'w')<br />for i in range(nbProbes): a.write(str(probings[0][i])+' '+str(probings[1][i])+' '+str(probings[2][i])+'\n')<br />a.close()<br /><br />#plot<br />from pylab import *<br />plot(probings[0,:]*sqrt(2),probings[1,:],'bo--')<br />ylabel(r'$\epsilon_{22}$')<br />xlabel(r'$\epsilon_{11} \times \sqrt{2}$')<br />title('response envelop')<br />grid()<br /><br />show()<br />

 

Reply in English

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Tag Cloud

%d bloggers like this: