New question #708020 on Yade: https://answers.launchpad.net/yade/+question/708020
Hello, I am having some difficulties with reassigning the rolling friction using the PeriTriaxController after compaction. In the example code below, I use an initial rolling friction 'eta_roll_init' of 0 (sliding fric also set to 0) to achieve a dense packing, then reassigning a new rolling friction 'eta_roll_final' at end of compaction. The issue seems to be that using different 'eta_roll_final' values yields identical specimen responses at the end of shearing. I am not sure exactly what the issue is and am wondering if someone with experience with this model can comment if I've doing the material and existing contact re-assignment correctly? For reference, I used these previous answers as a guide: https://answers.launchpad.net/yade/+question/701184 https://answers.launchpad.net/yade/+question/251218 https://answers.launchpad.net/yade/+question/251218 My understanding is that for existing contacts at end of compaction, I need to only update/re-calculate the maxRollPl value which is calculated as: min(rolling_friction * ball1_radius, rolling_friction * ball2_radius) Example Code: #================================================================================ from __future__ import print_function from yade import pack, qt, plot, export # Define model parameters density = 2650 # Density (kg/m3) young = 1e8 # Contact modulus (Pa) poisson = 0.25 # Shear to normal stiffness ratio, Ks/Kn fric_final = 20.0 # Sliding friction (deg) alpha_roll = 1.0 # Rolling resistance ratio, Kr/(R1*R2*Ks) eta_roll_final = 1.0 # Rolling friction (rad) # Define PSD and RVE parameters d_min, d_max = 20*0.13/1e3, 20*0.3/1e3 # ..log-linear Lx, Ly, Lz = 15*d_max, 15*d_max, 15*d_max numbins = 25 # Define specimen preparation parameters sigma_iso = -1e5 fric_init = 0.0 eta_roll_init = 0.0 packing = 0 max_strain_rate_iso = 0.20 max_unbalanced_iso = 1e-5 rel_stress_tol_iso = 1e-5 # Define controller parameters strain_target = 0.20 max_strain_rate_horiz = 1.0 max_strain_rate_vertical = 0.20 max_unbalanced_force = 10 # Define materials O.materials.append(CohFrictMat( density=density, # Density [kg/m3] young=young, # Particle modulus [Pa] poisson=poisson, # Ks/Kn ratio frictionAngle=radians(fric_init), # Local friction [rad] isCohesive=False, # Turn off adhesion momentRotationLaw=True, # Turn on rotational stiffness alphaKr=alpha_roll, # Dimensionless rolling stiffness alphaKtw=0, # Dimensionless twist stiffness etaRoll=eta_roll_init, # Rotational friction [rad] etaTwist=0, # Turn off twisting label="granr" # Material label )) # Define periodic boundaries O.periodic = True # Prepare specimen packing sp = pack.SpherePack() # ..log-linear distribution start, end = numpy.log10(d_min), numpy.log10(d_max) psdSizes = numpy.logspace(start, end, num=numbins, endpoint=True) psdCumm = numpy.linspace(0, 1., len(psdSizes)) sp.makeCloud((0, 0, 0), (Lx, Ly, Lz), psdSizes=psdSizes, psdCumm=psdCumm, periodic=True, seed=7) # Insert the packing sp.toSimulation(material="granr") # Define engines # more info: https://yade-dem.org/doc/introduction.html O.engines = [ # Reset forces ForceResetter(), # Approximate collision detection, create interactions InsertionSortCollider( # Sphere axis-aligned bounding box [Bo1_Sphere_Aabb()], label='collider'), # Handle interactions InteractionLoop( # Geom functors to create interaction geometry [Ig2_Sphere_Sphere_ScGeom()], # Phys functors to create interaction from materials [Ip2_FrictMat_FrictMat_FrictPhys(), Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()], # Law fcunctors to resolve interaction by computing forces per constituitive law [Law2_ScGeom_FrictPhys_CundallStrack(), Law2_ScGeom6D_CohFrictPhys_CohesionMoment()], label='interactionLoop'), GlobalStiffnessTimeStepper( # Use adaptive stiffness-based timestepper timestepSafetyCoefficient=0.8, timeStepUpdateInterval=100, # Set to True for domain decomp parallelMode=False, label='timeStepper' ), # Update positions uising Newton's equations NewtonIntegrator( # Non-viscous newton damping damping=.2, # Set gravity gravity=(0,0,0), # Create label for integrator engine label='newton' ), PeriTriaxController( # Create label/var name for triax engine label='triax', # Specify target values and whether they are strains or stresses goal=(sigma_iso, sigma_iso, sigma_iso), stressMask=7, # Type of servo-control dynCell=True, maxStrainRate=(max_strain_rate_iso, max_strain_rate_iso, max_strain_rate_iso), # Wait until the unbalanced force goes below this value maxUnbalanced=max_unbalanced_iso, relStressTol=rel_stress_tol_iso, # Call this function when goal is reached and the packing is stable doneHook='compactionFinished()' ), PyRunner(command='addPlotData()', iterPeriod=100), ] # Add custom data to plot def addPlotData(): plot.addData( unbalanced=unbalancedForce(), i=O.iter, sxx=triax.stress[0], syy=triax.stress[1], szz=triax.stress[2], exx=triax.strain[0], eyy=triax.strain[1], ezz=triax.strain[2], n0=yade._utils.porosity(), #n1=yade._utils.voxelPorosity(200,*yade._utils.aabbExtrema()), n1=yade._utils.voxelPorosity(200,yade._utils.aabbExtrema()[0]*1.05,yade._utils.aabbExtrema()[0]*0.95), # mechanical coordination number Zm Zm=utils.avgNumInteractions(skipFree=True), # coordination number Cn Cn=utils.avgNumInteractions(skipFree=False), # save all available energy data Etot=O.energy.total(), **O.energy ) # Enable energy tracking in the code O.trackEnergy = True # Define what to plot plot.plots = { 'i': ('unbalanced',), 'i ': ('sxx', 'syy', 'szz'), ' i': ('exx', 'eyy', 'ezz'), # Energy plot ' i ': (O.energy.keys, None, 'Etot'), ' i': ('Zm', None, 'n0', 'n1'), 'ezz': ('sxx','syy','szz') } # Show the plot plot.plot() def compactionFinished(): """ Defines hook for setting up shearing stage after compaction is finished. """ # Reassign friction values for shearing stage # ..for future contacts change material O.materials[0].frictionAngle = radians(fric_final) # radians O.materials[0].etaRoll = eta_roll_final # ..for existing contacts, set contact friction directly for i in O.interactions: i.phys.tangensOfFrictionAngle = tan(radians(fric_final)) i.phys.maxRollPl = min(eta_roll_final*O.bodies[i.id1].shape.radius, eta_roll_final*O.bodies[i.id2].shape.radius) # Set the current cell configuration to be the reference one O.cell.trsf = Matrix3.Identity # Change control type: keep constant confinement in x,y, 20% compression in z triax.goal = (sigma_iso, sigma_iso, -strain_target) triax.stressMask = 3 # Allow faster deformation along x,y to better maintain stresses triax.maxStrainRate = (max_strain_rate_horiz, max_strain_rate_horiz, max_strain_rate_vertical) # Next time, call triaxFinished instead of compactionFinished triax.doneHook = 'triaxFinished()' # Do not wait for stabilization before calling triaxFinished triax.maxUnbalanced = 10 def triaxFinished(): print('Test Finished') O.pause() O.run() -- You received this question notification because your team yade-users is an answer contact for Yade. _______________________________________________ Mailing list: https://launchpad.net/~yade-users Post to : yade-users@lists.launchpad.net Unsubscribe : https://launchpad.net/~yade-users More help : https://help.launchpad.net/ListHelp