API tools faq
paste
lingran

plastic energy dissipation

lingran
Aug 19th, 2013
130
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 7.55 KB | None | 0 0
raw download clone embed print report
  1. ### this simulation models impact between a boulder and a 2D soil layer.
  2. ### the soil layer is divided into several parts and plastic energy dissipation inside each parts are tracked in order to investigate the impact mechanisms.
  3.  
  4.  
  5.  
  6. ##### part 1 : generate a soil sample by gravity deposit
  7.  
  8. from yade import pack,plot
  9. from yade import utils
  10. import random
  11.  
  12.  
  13. finalFricDegree = 30
  14. O.trackEnergy=True
  15. young=100e6
  16. rmean=0.1
  17.  
  18. rmin=0.05
  19. rmax=0.15
  20.  
  21. mn=Vector3(0,0,-2*rmean)
  22. mx=Vector3(30*2*rmean,30*2*rmean,2*rmean)
  23.  
  24. number=500
  25. fallingHeight=5
  26.  
  27. O.materials.append(CohFrictMat(young=100e6,poisson=0.3,frictionAngle=radians(0),density=2600,normalCohesion=0,shearCohesion=0,alphaKr=0.5,alphaKtw=0.01,momentRotationLaw=False,isCohesive=False,etaRoll=-1,label='spheres'))
  28.  
  29. O.materials.append(FrictMat(young=100e6,poisson=0.3,frictionAngle=0,density=0,label='walls'))
  30.  
  31.  
  32. walls=utils.aabbWalls([mn,mx],thickness=0.5*rmean, oversizeFactor=1.0, material='walls')
  33. wallIds=O.bodies.append(walls)
  34.  
  35.  
  36.  
  37. ssp=pack.SpherePack()
  38.  
  39. for i in range(1,number):
  40. r =rmin + (rmax-rmin)*(-0.15*log(1-random.random()))**1.33333
  41. if r>rmax:
  42. r =rmin + (rmax-rmin)*(-0.15*log(1-random.random()))**1.33333
  43. mbegin=(0,0,-1*r)
  44. mend=(mx[0],mx[1],1*r)
  45. ssp.makeCloud(minCorner=mbegin,maxCorner=mend,rMean=r,rRelFuzz=0,num=1)
  46.  
  47. O.bodies.append([utils.sphere(center,rad,material='spheres') for center,rad in ssp])
  48.  
  49.  
  50. newton=NewtonIntegrator(damping=0.9,kinSplit=True,gravity=(0,0,0))
  51.  
  52. O.engines=[
  53. ForceResetter(),
  54. InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
  55. InteractionLoop(
  56. [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom()],
  57. [Ip2_FrictMat_FrictMat_FrictPhys(),Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
  58. [Law2_ScGeom_FrictPhys_CundallStrack(label='law1'),Law2_ScGeom6D_CohFrictPhys_CohesionMoment(
  59. useIncrementalForm=False,
  60. always_use_moment_law=False,
  61. label='law2')]
  62. ),
  63. GravityEngine(gravity=(0,-9.81,0)),
  64. newton
  65.  
  66. ]
  67.  
  68. for b in O.bodies:
  69. if isinstance(b.shape,Sphere): b.state.blockedDOFs='zXY' # blocked movement along Z
  70.  
  71.  
  72. O.dt=utils.PWaveTimeStep()
  73.  
  74. O.run(150000, True)
  75.  
  76. ## reset boundaries
  77. Mz=max([b.state.pos[1]+b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)])
  78. h=Mz/4
  79.  
  80. mn=(0,0,-2*rmean)
  81. mxafter=Vector3(mx[0],Mz+2*rmean,2*rmean)
  82.  
  83. O.bodies.erase(0)
  84. O.bodies.erase(1)
  85. O.bodies.erase(2)
  86. O.bodies.erase(3)
  87. O.bodies.erase(4)
  88. O.bodies.erase(5)
  89.  
  90. O.bodies.append(utils.aabbWalls([mn,mxafter],thickness=0.1*rmean, oversizeFactor=1.0, material='walls'))
  91.  
  92. O.bodies.erase(3)
  93. O.bodies.erase(4)
  94. O.bodies.erase(5)
  95.  
  96. O.run(10000, True)
  97.  
  98. ## add friction and moment law
  99. O.materials[0].frictionAngle=radians(finalFricDegree)
  100. O.materials[0].momentRotationLaw=True
  101. O.engines[2].lawDispatcher.functors[1].always_use_moment_law = True
  102.  
  103. O.engines[4].damping=0.000
  104.  
  105. O.run(80000, True)
  106. O.saveTmp()
  107.  
  108. ##### part 2: divide the soil layer into different parts
  109.  
  110.  
  111. mxafter=Vector3(mx[0],Mz,2*rmean)
  112.  
  113. centerx=0.5*mxafter[0]
  114. centery=mxafter[1]
  115.  
  116.  
  117.  
  118. O.materials.append(CohFrictMat(young=100e6,poisson=0.3,frictionAngle=radians(finalFricDegree),density=2600,normalCohesion=0,shearCohesion=0,alphaKr=0.5,alphaKtw=0.01,
  119. momentRotationLaw=True,isCohesive=False,etaRoll=-1,label='heart'))
  120.  
  121. for i in range(1,4):
  122. for j in range(1,3):
  123. O.materials.append(CohFrictMat(young=100e6,poisson=0.3,frictionAngle=radians(finalFricDegree),density=2600,normalCohesion=0,shearCohesion=0,alphaKr=0.5,alphaKtw=0.01,
  124. momentRotationLaw=True,isCohesive=False,etaRoll=-1,label='crown_'+str(i)+str(j)))
  125.  
  126.  
  127. import numpy,random
  128. a=numpy.zeros(6).reshape((3, 2))
  129.  
  130.  
  131. for i in O.bodies:
  132. x=i.state.pos[0]
  133. y=i.state.pos[1]
  134. D=sqrt((x-centerx)**2+(y-centery)**2)
  135. xy=atan2((y-centery),(x-centerx))
  136. k=0
  137. if (D/h <= 1):
  138. i.material=O.materials['heart']
  139. i.shape.color=Vector3(1,0,0)
  140. if (1 < D/h < 4):
  141. j=int(D/h)
  142. if -2*pi/3 <= xy <= -pi/3 :
  143. k=1
  144. else:
  145. k=0
  146. a[j-1][k]+=1
  147. i.material=O.materials['crown_'+str(j)+str(k+1)]
  148.  
  149.  
  150. ## mark the soil layers with different colors
  151. color=Vector3(random.random(),random.random(),random.random())
  152.  
  153. for i in O.bodies:
  154. i.shape.color=color
  155.  
  156. for i in O.bodies:
  157.  
  158. if i.material==O.materials['crown_11']:
  159. i.shape.color=Vector3(0.2,0.5,0)
  160.  
  161. if i.material==O.materials['crown_12']:
  162. i.shape.color=Vector3(0.6,0.2,1)
  163.  
  164. if i.material==O.materials['crown_21']:
  165. i.shape.color=Vector3(1,0.5,0.6)
  166.  
  167. if i.material==O.materials['crown_22']:
  168. i.shape.color=Vector3(1,1,0)
  169.  
  170. if i.material==O.materials['crown_31']:
  171. i.shape.color=Vector3(0.2,0.9,0.3)
  172.  
  173. if i.material==O.materials['crown_32']:
  174.  
  175. i.shape.color=Vector3(0.6,1,1)
  176.  
  177.  
  178. if i.material==O.materials['heart']:
  179. i.shape.color=Vector3(1,0,0)
  180.  
  181.  
  182. Gl1_Sphere.stripes=1
  183. O.run(10000,True)
  184. O.save('2dlayer.xml')
  185.  
  186. #### part 3: add the boulder and conduct impact modelling, track energy dissipations inside each soil layer
  187.  
  188. boulderradius=3*rmean
  189.  
  190. sp2=pack.SpherePack()
  191.  
  192. mbegin=(0.5*mxafter[0]*1.0-boulderradius,mxafter[1],-1*boulderradius)
  193. mend=(0.5*mxafter[0]*1.0+boulderradius,mxafter[1]+2
  194. *boulderradius,1*boulderradius)
  195. sp2.makeCloud(minCorner=mbegin,maxCorner=mend,rMean=boulderradius,rRelFuzz=0,num=1)
  196. O.bodies.append([utils.sphere(center,rad,material='spheres') for center,rad in sp2])
  197.  
  198. O.bodies[3].state.blockedDOFs='zXY'
  199.  
  200. a=sin(0*3.14/180)
  201. b=cos(0*3.14/180)
  202. O.bodies[3].state.vel=Vector3(-sqrt(2*9.8*fallingHeight)*a,-sqrt(2*9.8*fallingHeight)*b,0)
  203.  
  204. O.dt=0.1*utils.PWaveTimeStep()
  205.  
  206. for i in range (1,200):
  207. O.run(100,True)
  208. if -O.forces.f(3)[1]< O.bodies[3].state.mass*9.81:
  209. tt=O.time
  210. break
  211.  
  212. ## time when the impact begins
  213. tt=O.time
  214.  
  215. O.saveTmp()
  216.  
  217. aa=dict(O.energy.items())
  218. E=aa['plastDissip']
  219. f=open('eptrack.txt','w')
  220. p_11=0
  221. p_12=0
  222. p_21=0
  223. p_22=0
  224. p_31=0
  225. p_32=0
  226. p=0
  227. psum=0
  228.  
  229. for k in range (1,4000):
  230.  
  231. O.run(1,True)
  232. aa=dict(O.energy.items())
  233. for i in O.interactions:
  234.  
  235. Ft=i.phys.shearForce
  236. du=i.geom.shearInc
  237. kt=i.phys.ks
  238. Fn=i.phys.normalForce
  239. Ftt=Ft-kt*du
  240. maxFs=Fn.norm()*i.phys.tangensOfFrictionAngle
  241. if(Ftt.norm() > maxFs):
  242. ratio = maxFs / Ftt.norm()
  243. trialForce=Ftt
  244. Ftt *= ratio
  245. e=((1/kt)*(trialForce-Ftt)).dot(Ftt)
  246. psum+=e
  247. if (O.bodies[i.id1].material== O.materials['crown_11']) or (O.bodies[i.id2].material== O.materials['crown_11']):
  248. p_11+=e
  249. elif (O.bodies[i.id1].material== O.materials['crown_12']) or (O.bodies[i.id2].material== O.materials['crown_12']):
  250. p_12+=e
  251. elif (O.bodies[i.id1].material== O.materials['crown_21']) or (O.bodies[i.id2].material== O.materials['crown_21']):
  252. p_21+=e
  253. elif (O.bodies[i.id1].material== O.materials['crown_22']) or (O.bodies[i.id2].material== O.materials['crown_22']):
  254. p_22+=e
  255. elif (O.bodies[i.id1].material== O.materials['crown_31']) or (O.bodies[i.id2].material== O.materials['crown_31']):
  256. p_31+=e
  257. elif (O.bodies[i.id1].material== O.materials['crown_32']) or (O.bodies[i.id2].material== O.materials['crown_32']):
  258. p_32+=e
  259. else:
  260. p+=e
  261. f.write('%f %f %f %f %f %f %f %f %f %f\n' %(O.time-tt,p_11,p_12,p_21,p_22,p_31,p_32,p,psum,aa['plastDissip']-E))
  262.  
  263.  
  264. f.close()
Advertisement
Add Comment
Please, Sign In to add comment
Public Pastes
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!