API tools faq
paste
Advertisement
nicogodet

Untitled

nicogodet
Jun 28th, 2018
39
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Python 16.75 KB | None | 0 0
raw download clone embed print report
  1. ########################################## encoding: utf-8 #######################################
  2. from yade.gridpfacet import *
  3. from yade import pack,plot
  4. import math
  5. import random as rand
  6. import numpy as np
  7. from scipy.stats import gaussian_kde
  8.  
  9.  
  10. ## Definition of several variables (geometry, density, etc etc)
  11.  
  12. ########################################## MATERIAL #######################################
  13.  
  14. ## I remove the value of the variables
  15. O.materials.append( #Viscous flow material
  16.         FrictMat(
  17.             density=sphDensity,
  18.             frictionAngle=math.radians(sphFrictAng),
  19.             poisson=sphPoisson,
  20.             young=sphYoung,
  21.             label='Frict'
  22.         ))
  23.  
  24. O.materials.append( #Bed channel material
  25.                 FrictMat(
  26.                         density=2500,
  27.                         frictionAngle=math.radians(chanFrictAng),
  28.                         poisson=0.2,
  29.                         young=1e7,
  30.                         label='FrictBase'
  31.                 ))
  32.  
  33. O.materials.append( #Wall channel material
  34.                 FrictMat(
  35.                         density=2500,
  36.                         frictionAngle=math.radians(0.),
  37.                         poisson=0.2,
  38.                         young=1e7,
  39.                         label='FrictWall'
  40.                 ))
  41.  
  42. O.materials.append( #Not used here
  43.                 FrictMat(
  44.                         density=2500,
  45.                         frictionAngle=math.radians(20.),
  46.                         poisson=0.2,
  47.                         young=1e7,
  48.                         label='FrictRigid'
  49.                 ))
  50.  
  51.  
  52.     ###################### FENCE #####################
  53.  
  54. ## I remove the value of the variables
  55.  
  56. O.materials.append( WireMat( young=young_wire,poisson=0.3,frictionAngle=radians(frict_wire),density=7850,diameter=2*Rfilet,isDoubleTwist=False,label='gridNodeNetPackSimple',strainStressValues=SSVal_Simple) ) # Material to create the GridNode
  57.  
  58. O.materials.append( WireMat( young=young_wire,poisson=0.3,frictionAngle=radians(frict_wire),density=7850,diameter=2*Rfilet,isDoubleTwist=False,label='gridNodeNetPackDouble',strainStressValues=SSVal_Simple) ) # Material to create the GridNode
  59.  
  60. O.materials.append(FrictMat( young=young_Conn,poisson=0.3,density=7850,frictionAngle=radians(frict_Con), label='CohMat'))#normalCohesion=200e6,isCohesive=True,shearCohesion=200e2,momentRotationLaw=False, label='CohMat') ) # Material to create the iteraction (Gridconection)
  61.  
  62. ######################################### ENGINES #######################################
  63.  
  64. # Definition of simulation's variable
  65. Factor = 1.5    #Factor ratio for the distance detection of distant interactions
  66. theta_method = 0.55 #1 : Euler Backward, 0.5 : trapezoidal rule, 0 : not used, 0.55 : suggested optimum
  67. resolution = 1      #Change normal component resolution method, 0: Iterative exact resolution (theta method, linear contact), 1: Newton-Rafson dimentionless resolution (theta method, linear contact), 2: Newton-Rafson with nonlinear surface deflection (Hertzian-like contact)
  68. NewtonRafsonTol = 1e-6  #Precision of the resolution of Newton-Rafson's method
  69. vel_front = Vector3(0.,0.,0.)
  70.  
  71. O.engines=[
  72.     ForceResetter(),
  73.     InsertionSortCollider([
  74.         Bo1_Sphere_Aabb(aabbEnlargeFactor=Factor, label="aabb"),
  75.         Bo1_GridConnection_Aabb(),
  76.         Bo1_Wall_Aabb(),
  77.         Bo1_Box_Aabb()],
  78.         verletDist=-0.1,
  79.         allowBiggerThanPeriod=False
  80.     ), 
  81.     InteractionLoop([
  82.         Ig2_Sphere_Sphere_ScGeom(),
  83.         Ig2_GridNode_GridNode_GridNodeGeom6D(),
  84.         Ig2_GridConnection_GridConnection_GridCoGridCoGeom(),  
  85.         Ig2_Sphere_GridConnection_ScGridCoGeom(),      
  86.         Ig2_Sphere_Sphere_ScGeom6D(interactionDetectionFactor=Factor,label="Ig2"),
  87.         Ig2_Box_Sphere_ScGeom6D(),
  88.     ], 
  89.     [
  90.         Ip2_WireMat_WireMat_WirePhys(label='wire_wire'),
  91.         Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True,setCohesionOnNewContacts=False),
  92.        
  93.        
  94.         Ip2_FrictMat_FrictMat_LubricationPhys(eta=m_viscosity,eps=m_roughness,label="Phys_Lub"),
  95.         Ip2_FrictMat_FrictMat_FrictPhys(),
  96.         ## Here I don't know why but changing the order solved one problem....
  97.        
  98.     ],
  99.     [
  100.         Law2_ScGeom6D_CohFrictPhys_CohesionMoment(),
  101.         Law2_ScGeom_FrictPhys_CundallStrack(),
  102.         Law2_ScGridCoGeom_FrictPhys_CundallStrack(),
  103.         Law2_ScGeom_WirePhys_WirePM(),
  104.         Law2_GridCoGridCoGeom_FrictPhys_CundallStrack(),
  105.         Law2_ScGeom_ImplicitLubricationPhys(
  106.                         activateNormalLubrication=True,
  107.                         activateTangencialLubrication=True,
  108.                         activateTwistLubrication=True,
  109.                         activateRollLubrication=True,
  110.                         theta=theta_method,
  111.                         resolution=resolution,
  112.                         warnedOnce=True,
  113.                         maxSubSteps=20,
  114.                         NewtonRafsonTol=NewtonRafsonTol,
  115.                         debug=False)
  116.  
  117.     ]
  118.     ),
  119.     DragEngine(Rho=1200,Cd=0.47,V=vel_front, label='Drag', dead=True), #Modified engine to include computation of drag force using relative velocity based on V instead of absolute velocity of the particle
  120.     NewtonIntegrator(gravity=gravityVector,damping=0.3,label='newton'),
  121.     #Several PyRunner
  122.     VTKRecorder(fileName=folder+'/VTK/vtk-',virtPeriod=0.1),
  123. ]
  124.  
  125. Drag.ids=spIds
  126. Drag.Rho=fluidDensity
  127. Drag.Cd=0.47
  128.  
  129.  
  130. ######################################### CONSTRUCTION DU BASSIN #######################################
  131.  
  132. boxColor = (104/255., 111/255., 140/255.)
  133. chanBack = O.bodies.append(aabbWalls([(0,length_channel/3.,0),(0,length_channel/3.,height_channel)],thickness=0.0,material='FrictWall',oversizeFactor=1.0,color=boxColor))
  134. chanS1 = O.bodies.append(aabbWalls([(0,length_channel/3.,0),(0,0,height_channel)],thickness=0.0,material='FrictWall',oversizeFactor=1.0,color=boxColor))
  135. chanB = O.bodies.append(aabbWalls([(0,length_channel/3.,0),(width_channel,0,0)],thickness=0.0,material='FrictBase',oversizeFactor=1.0,color=boxColor))
  136. chanS2 = O.bodies.append(aabbWalls([(width_channel,0,0),(width_channel,length_channel/3.,height_channel)],thickness=0.0,material='FrictWall',oversizeFactor=1.0,color=boxColor))
  137. chanBottom = O.bodies.append(aabbWalls([(0,0,0),(width_channel,0,height_channel)],thickness=0.0,material='FrictWall',oversizeFactor=1.0,color=boxColor))
  138.  
  139.  
  140. DepBox = O.bodies.append(aabbWalls([(0.,length_channel/3.,-0.1),(width_channel,length_channel/3.-0.1*length_channel,2*height_channel)],thickness=0.0,material='FrictWall',oversizeFactor=1.0,color=boxColor))
  141.  
  142. ######################################### AJOUT ECOULEMENT #######################################
  143.  
  144. partCloud = pack.SpherePack()
  145. partCloud.makeCloud(minCorner=(0.,length_channel/3.-length_channel/10./2.,0.),maxCorner=(width_channel,length_channel/3.-0.095*length_channel,height_channel*2.), rMean=diameterPart/2.0)
  146. spIds = partCloud.toSimulation(material='Frict',color=[205./255,175./255,149./255])
  147. allSp = spIds
  148.  
  149. ######################################### CONSTRUCTION DU MAILLAGE #######################################
  150.  
  151.     ##################### TAILLE REELLE #####################
  152.  
  153.  
  154. for i in range(Nx+1):               # Loops x
  155.     x = xstart + Lx*i
  156.         for j in range(kz):         # loops z
  157.         z = zstart + Lz*j
  158. ################################### ACCROCHE COTE GAUCHE ######################################
  159.     # Note: the radius is 1.01*Rfilet and 0.99*Rfilet because in that way when you create a grid connection betwen two GridNode with diferent characteristics it  takes the inf of the bigest radius
  160.         if i == 0:
  161.                     NetPack += [O.bodies.append( gridNode([x+Lx/4,y,z+Lz/2],1.01*Rfilet, color=[0,0,1],material='gridNodeNetPackDouble',wire=False,highlight=False,fixed=False))]
  162.                     ItLeft.append(NetPack[-1])
  163.                     Agraphes.append(NetPack[-1])
  164.                     if j != 0:
  165.                         NetPack += [ O.bodies.append( gridNode([x,y,z],0.99*Rfilet, color=[1,0,0],material='gridNodeNetPackSimple',wire=False,highlight=False,fixed=False))]
  166.                 ItLeft.append(NetPack[-1])
  167.  
  168.     ##################### ACCROCHE COTE DROIT #####################
  169.         if i == Nx:
  170.                     NetPack += [ O.bodies.append( gridNode([x-Lx/4,y,z],1.01*Rfilet, color=[0,0,1],material='gridNodeNetPackDouble',wire=False,highlight=False,fixed=False))]
  171.                     Agraphes.append(NetPack[-1])
  172.                     ItRight.append(NetPack[-1])
  173.                     if j != kz-1:
  174.                         NetPack += [ O.bodies.append( gridNode([x,y,z+Lz/2],0.99*Rfilet, color=[1,0,0],material='gridNodeNetPackSimple',wire=False,highlight=False,fixed=False))]
  175.                         ItRight.append(NetPack[-1])
  176.  
  177.                     if j == kz-1:
  178.                         if Nz%2 != 0:
  179.                             NetPack +=[ O.bodies.append( gridNode([x,y,z+Lz/2],0.99*Rfilet, color=[1,0,0],material='CohMat',wire=False,highlight=False,fixed=False))]
  180.                             ItRight.append(NetPack[-1])
  181.                             NetPack +=[ O.bodies.append( gridNode([x-Lx/4,y,z+Lz],1.01*Rfilet, color=[0,0,1],material='CohMat',wire=False,highlight=False,fixed=False))]
  182.                             Agraphes.append(NetPack[-1])
  183.                             ItRight.append(NetPack[-1])
  184.  
  185.     ##################### RESTE MAILLAGE #####################
  186.                 if i!= 0 and i!= Nx:
  187.                     NetPack += [ O.bodies.append( gridNode([x-Lx/4,y,z],1.01*Rfilet, color=[0,0,1],material='gridNodeNetPackDouble',wire=False,highlight=False,fixed=False))]
  188.                     Agraphes.append(NetPack[-1])
  189.                     ItCenter.append(NetPack[-1])
  190.                
  191.                     NetPack +=[ O.bodies.append( gridNode([x,y,z+Lz/4],0.99*Rfilet, color=[0,1,0],material='gridNodeNetPackSimple',wire=False,highlight=False,fixed=False))]
  192.                     ItCenter.append(NetPack[-1])
  193.                     NetPack += [ O.bodies.append( gridNode([x+Lx/4,y,z+Lz/2],1.01*Rfilet, color=[0,0,1],material='gridNodeNetPackDouble',wire=False,highlight=False,fixed=False))]
  194.                     Agraphes.append(NetPack[-1])
  195.                     ItCenter.append(NetPack[-1])
  196.                
  197.                     if j != kz-1:
  198.                         NetPack += [ O.bodies.append( gridNode([x,y,z+3*Lz/4],0.99*Rfilet, color=[0,1,0],material='gridNodeNetPackSimple',wire=False,highlight=False,fixed=False))]
  199.                         ItCenter.append(NetPack[-1])
  200.                     if j == kz-1:
  201.                         if  Nz%2!=0:
  202.                              NetPack += [ O.bodies.append( gridNode([x,y,z+3*Lz/4],0.99*Rfilet, color=[0,1,0],material='gridNodeNetPackSimple',wire=False,highlight=False,fixed=False))]
  203.                              ItCenter.append(NetPack[-1])
  204.                              NetPack += [ O.bodies.append( gridNode([x-Lx/4,y,z+Lz],1.01*Rfilet, color=[0,1,0],material='gridNodeNetPackDouble',wire=False,highlight=False,fixed=False))]
  205.                              Agraphes.append(NetPack[-1])
  206.                              ItCenter.append(NetPack[-1])
  207.            
  208.  
  209.     ##################### DEFINITION INTERACTIONS #####################
  210.  
  211. for i,j in zip(ItLeft[:-1], ItLeft[1:]):#it gauche
  212.     O.bodies.append( gridConnection(i,j,Rfilet,material='CohMat',color=[1,2,1],mask=3))
  213. for i,j in zip(ItCenter[:-1], ItCenter[1:]):#it centre
  214.     if O.bodies[i].state.pos[2] != zmax:
  215.         O.bodies.append( gridConnection(i,j,Rfilet,material='CohMat',color=[1,2,1],mask=3))
  216. for i,j in zip(ItRight[:-1], ItRight[1:]):#it Droite
  217.         O.bodies.append( gridConnection(i,j,Rfilet,material='CohMat',color=[1,2,1],mask=3))
  218.  
  219. for i,j in zip(ItCenter[:-2], ItCenter[2:]) :#it centre verticale
  220.         if O.bodies[i].state.pos[0] == O.bodies[j].state.pos[0]:
  221.                 O.bodies.append( gridConnection(i,j,Rfilet,material='CohMat',color=[1,2,1],mask=3))
  222. for i, j in zip(Agraphes[:-Nx], Agraphes[Nx:]):
  223.         O.bodies.append( gridConnection(i,j,Rfilet,material='CohMat',color=[1,2,1],mask=3))
  224.  
  225. ## The entiere mesh is built like above, first GridNode using CohFrictMat then the GridConnections using FrictMat
  226. ## I will not show the entire routine, I'm not sure yet I can show its but it only creates the geometry using GridNode and GridConnections as above.
  227.  
  228.  
  229. def Freins(Pos1, Pos2, Fdec,radius):
  230.    
  231.     ##I remove the value of the variables
  232.  
  233.     freinmat = O.materials.append(WireMat(young=youngf, poisson=poissonf,frictionAngle=radians(10),density=densityf,isDoubleTwist=False,diameter=2*radius,strainStressValues=strainStressValuesf))
  234.    
  235.     frein = []
  236.     frein.extend([
  237.     O.bodies.append(utils.sphere(Pos1,radius,wire=False,fixed=False,material=freinmat,color=colorf,mask=3)),
  238.     O.bodies.append(utils.sphere(Pos2,radius,wire=False,fixed=False,material=freinmat,color=colorf,mask=3))])
  239.     createInteraction(frein[0],frein[1])
  240.  
  241.     return (frein)
  242.  
  243.  
  244. ################################### MANILLE ##################################
  245.  
  246. def Manille_Vert(Rman, ouv, posX, posY, posZ, base):
  247.  
  248.     ## I remove the value of the variables
  249.  
  250.     bouclemat = O.materials.append(CohFrictMat(young=youngb,poisson=poissonb,density=densityb,frictionAngle=radians(frottement),normalCohesion=pi*strenruptb,shearCohesion=shearCohb,momentRotationLaw=True))
  251.     bouclecomat = O.materials.append(FrictMat(young=youngb,poisson=poissonb,density=densityb,frictionAngle=radians(frottement)))
  252.    
  253.     boucle.extend([
  254.         O.bodies.append(gridNode([posX,posY,posZ+dec],Rman,wire=wire,fixed=False,material=bouclemat,color=colorb)),
  255.         O.bodies.append(gridNode([posX,posY+dec,posZ],Rman,wire=wire,fixed=False,material=bouclemat,color=colorb)),
  256.         O.bodies.append(gridNode([posX,posY,posZ-dec],Rman,wire=wire,fixed=False,material=bouclemat,color=colorb)),
  257.         O.bodies.append(gridNode([posX,posY-dec,posZ],Rman,wire=wire,fixed=False,material=bouclemat,color=colorb)) ])
  258.  
  259.     conboucle=[]
  260.     for i,j in zip(boucle[:-1],boucle[1:]):
  261.         conboucle.append( O.bodies.append( gridConnection(i,j,Rman,wire=wire, material = bouclecomat,color=colorb,mask=3) ) )
  262.     conboucle.append( O.bodies.append( gridConnection(boucle[0],boucle[-1],Rman,wire=wire, material = bouclecomat,color=colorb,mask=3) ) )
  263.  
  264.         ############ CLUMP ############
  265.     clumping = []
  266.     clumping.extend([boucle[0],boucle[1],boucle[2],boucle[3],base])
  267.     clumps = O.bodies.clump(clumping)
  268.  
  269.     return boucle,conboucle,clumps
  270.    
  271. def Manille_frein(Rman, ouv, posX, posY, posZ,Dir):
  272.  
  273.     ## I remove the value of the variables
  274.  
  275.     bouclemat = O.materials.append(CohFrictMat(young=youngb,poisson=poissonb,density=densityb,frictionAngle=radians(frottement),normalCohesion=pi*strenruptb,shearCohesion=shearCohb,momentRotationLaw=True))
  276.     bouclecomat = O.materials.append(FrictMat(young=youngb,poisson=poissonb,density=densityb,frictionAngle=radians(frottement)))
  277.    
  278.     if Dir=='Left':
  279.         boucle.extend([
  280.             O.bodies.append(gridNode([posX-dec*cos(2*Lz/Lx),posY,posZ+dec*sin(2*Lz/Lx)],Rman,wire=wire,fixed=False,material=bouclemat,color=colorb)),
  281.             O.bodies.append(gridNode([posX,posY+dec,posZ],Rman,wire=wire,fixed=False,material=bouclemat,color=colorb)),
  282.             O.bodies.append(gridNode([posX+dec*cos(2*Lz/Lx),posY,posZ-dec*sin(2*Lz/Lx)],Rman,wire=wire,fixed=False,material=bouclemat,color=colorb)),
  283.             O.bodies.append(gridNode([posX,posY-dec,posZ],Rman,wire=wire,fixed=False,material=bouclemat,color=colorb)) ])
  284.     else:
  285.         boucle.extend([
  286.             O.bodies.append(gridNode([posX+dec*cos(2*Lz/Lx),posY,posZ+dec*sin(2*Lz/Lx)],Rman,wire=wire,fixed=False,material=bouclemat,color=colorb)),
  287.             O.bodies.append(gridNode([posX,posY+dec,posZ],Rman,wire=wire,fixed=False,material=bouclemat,color=colorb)),
  288.             O.bodies.append(gridNode([posX-dec*cos(2*Lz/Lx),posY,posZ-dec*sin(2*Lz/Lx)],Rman,wire=wire,fixed=False,material=bouclemat,color=colorb)),
  289.             O.bodies.append(gridNode([posX,posY-dec,posZ],Rman,wire=wire,fixed=False,material=bouclemat,color=colorb)) ])
  290.  
  291.     frein = []
  292.     if Dir=='Left':
  293.         frein.extend([
  294.         O.bodies.append(utils.sphere([O.bodies[boucle[0]].state.pos[0],O.bodies[boucle[0]].state.pos[1],O.bodies[boucle[0]].state.pos[2]*1.001],Rman,wire=False,fixed=False,material=freinmat,color=colorf,mask=3)),
  295.         O.bodies.append(utils.sphere([i+j for i,j in zip(O.bodies[boucle[0]].state.pos,[-dec*cos(2*Lz/Lx),0,dec*sin(2*Lz/Lx)])],Rman,wire=False,fixed=False,material=freinmat,color=colorf,mask=3))])
  296.         O.bodies[-1].shape.color.color=[0.5,0.5,0.5]
  297.     else:
  298.         frein.extend([
  299.         O.bodies.append(utils.sphere([O.bodies[boucle[0]].state.pos[0],O.bodies[boucle[0]].state.pos[1],O.bodies[boucle[0]].state.pos[2]*1.001],Rman,wire=False,fixed=False,material=freinmat,color=colorf,mask=3)),
  300.         O.bodies.append(utils.sphere([i+j for i,j in zip(O.bodies[boucle[0]].state.pos,[dec*cos(2*Lz/Lx),0,dec*sin(2*Lz/Lx)])],Rman,wire=False,fixed=False,material=freinmat,color=colorf,mask=3))])
  301.     createInteraction(frein[0],frein[1])
  302.  
  303.     conboucle=[]
  304.     for i,j in zip(boucle[:-1],boucle[1:]):
  305.         conboucle.append( O.bodies.append( gridConnection(i,j,Rman,wire=wire, material = bouclecomat,color=colorb,mask=3) ) )
  306.     conboucle.append( O.bodies.append( gridConnection(boucle[0],boucle[-1],Rman,wire=wire, material = bouclecomat,color=colorb,mask=3) ) )
  307.  
  308.         ############ CLUMP ############
  309.     clumping = []
  310.     clumping.extend([boucle[0],boucle[1],boucle[2],boucle[3],frein[0]])
  311.     clumps = O.bodies.clump(clumping)
  312.     TTT=frein[1]
  313.     return boucle,conboucle,clumps,TTT
  314. boucle,conboucle,shackles=[],[],[]
  315.  
  316.  
  317.  
  318.  
  319. ########################## Functions ##############################
  320.  
  321. ## Lot of functions here that I removed because not relevant to the problem encountered
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!