1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
|
# -*- encoding=utf-8 -*-
# 2022 © Vasileios Angelidakis <vasileios.angelidakis@ncl.ac.uk>
# Model multi-drum column under harmonic excitation, using the Potential Blocks
from yade import pack
from potential_utils import *
import math, os, errno
from numpy import array
# -----------------------------------------------------------------------------
m = O.materials.append(FrictMat(young=-1, poisson=-1, frictionAngle=atan(0.4), density=2000, label='frictionless'))
Kn = 1e8 #Pa/m
Ks = Kn * 2 / 3 #Pa/m
# -----------------------------------------------------------------------------
# Engines
O.engines = [
ForceResetter(),
InsertionSortCollider([PotentialBlock2AABB()], verletDist=0.01, avoidSelfInteractionMask=2),
InteractionLoop(
[Ig2_PB_PB_ScGeom(twoDimension=False, unitWidth2D=1.0, calContactArea=True)],
[Ip2_FrictMat_FrictMat_KnKsPBPhys(kn_i=Kn, ks_i=Ks, Knormal=Kn, Kshear=Ks, useFaceProperties=False, viscousDamping=0.3)],
[Law2_SCG_KnKsPBPhys_KnKsPBLaw(label='law', neverErase=False, allowViscousAttraction=False)]
),
NewtonIntegrator(damping=0.0, exactAsphericalRot=True, gravity=[0, 0, -9.81], label='newton'),
]
# -----------------------------------------------------------------------------
# Create drum column with variable cross-sections
height = 10.08
maxR = 1.65
minR = 1.28
#drumHeight=0.25
drumNo = 10
drumHeight = height / drumNo
for i in range(0, drumNo):
zBtm = (1 - (i + 0) * drumHeight / height)
zTop = (1 - (i + 1) * drumHeight / height)
btmRadius = minR + zBtm * (maxR - minR)
topRadius = minR + zTop * (maxR - minR)
b = prism(
material=m,
radius1=btmRadius,
radius2=topRadius,
thickness=drumHeight,
numFaces=10,
r=0.01,
R=0.0,
center=None,
color=[-1, -1, -1],
mask=1,
isBoundary=False,
fixed=False
)
b.state.pos = b.shape.position + [0, 0, i * (drumHeight)] + [0, 0, drumHeight / 2]
O.bodies.append(b)
maxR = 1.28
minR = 1.65
for i in range(drumNo, drumNo + 2):
zBtm = (1 - (i + 1) * drumHeight / height)
zTop = (1 - (i + 0) * drumHeight / height)
btmRadius = maxR + zBtm * (maxR - minR)
topRadius = maxR + zTop * (maxR - minR)
b = prism(
material=m,
radius1=topRadius,
radius2=btmRadius,
thickness=drumHeight,
numFaces=10,
r=0.01,
R=0.0,
center=None,
color=[-1, -1, -1],
mask=1,
isBoundary=False,
fixed=False
)
b.state.pos = b.shape.position + [0, 0, i * (drumHeight)] + [0, 0, drumHeight / 2]
O.bodies.append(b)
# Top plate
length = 3 * minR
thickness = maxR / 1.5
pos = O.bodies[-1].state.pos + [0, 0, drumHeight / 2 + thickness / 2.]
pb = cuboid(material=m, edges=[length, length, thickness], r=0.01, R=0.0, center=pos, mask=1, isBoundary=False, fixed=False, color=[1, 0, 0])
O.bodies.append(pb)
# Bottom plate
length = 10 * minR
thickness = maxR
pos = [0, 0, -thickness / 2]
pb = cuboid(material=m, edges=[length, length, thickness], r=0.01, R=0.0, center=pos, mask=1, isBoundary=False, fixed=True, color=[0, 1, 0])
pb_ID = O.bodies.append(pb)
O.engines = O.engines + [HarmonicMotionEngine(A=[0.1, 0, 0], f=[3, 0, 0], fi=[pi / 2, 0, 0], ids=[pb_ID], label='hME')] # Harmonic Motion Engine
# -----------------------------------------------------------------------------
# Visualisation
from yade import qt
qt.Controller()
v = qt.View()
v.sceneRadius = 20
v.ortho = True
v.eyePosition = Vector3(0, -23, 6)
v.upVector = Vector3(0, 0, 1)
v.viewDir = Vector3(0, 1, 0)
rndr = yade.qt.Renderer()
rndr.wire = True
# -----------------------------------------------------------------------------
# Time step
O.dt = 0.001 * sqrt(O.bodies[0].state.mass / Kn)
O.saveTmp()
O.run()
|
