SpringDamperMasspointSystem.py
You can view and download this file on Github: SpringDamperMasspointSystem.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is a EXUDYN example
3#
4# Details: The 3D movement of a point mass system is simulated.
5# The point masses are connected by spring-damper elements.
6# The system ist modelled in accordance with the Bachelor thesis of Thomas
7# Pfurtscheller (SS/WS19)
8#
9# Author: Holzinger Stefan
10# Date: 2019-10-14
11#
12# Copyright:This file is part of Exudyn. Exudyn is free software. You can redistribute it and/or modify it under the terms of the Exudyn license. See 'LICENSE.txt' for more details.
13#
14#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
15
16import exudyn as exu
17from exudyn.itemInterface import *
18import exudyn
19from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
20import exudyn.graphics as graphics
21
22import numpy as np
23
24print('EXUDYN version='+exu.config.Version())
25
26#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
27# Modelling of MBS starts here
28
29# define number of nodes in mesh
30numberOfNodesX = 6
31numberOfNodesY = 4
32numberOfNodes = numberOfNodesX*numberOfNodesY
33
34# define undeformed spring length
35springDamperElementLength = 1
36springDamperDiagnoalElementLength = np.sqrt( springDamperElementLength*springDamperElementLength + springDamperElementLength*springDamperElementLength )
37
38# define mass point mass and spring-damper properties
39massMassPoint = 1/2 # each node has half mass, since elements have mass
40elementStiffness = 1000
41elementDamping = 5
42activateElement = 1
43
44# Flags
45clambRightSide = True # clamb right side of mesh. If false, only left side of system is clambed
46useSpringDamperElementsX = True
47useSpringDamperElementsY = True
48useSpringDamperElementsDiagonal = True
49useSpringDamperElementsOffDiagonal = True
50createSolutionFile = False
51
52# this part of the code is used for convergence analysis
53# if hend < h0, step size will be reduced until this condition is false
54h0 = 1.0e-5
55hend = 1.0e-5
56stepSizeReductionFactor = 2 # factor which is used to reduce time step size (increase number of steps)
57
58stepSizeList = list()
59stepSizeList.append(h0)
60ctr1 = 0
61while stepSizeList[ctr1] > hend:
62 stepSizeList.append( stepSizeList[ctr1]/stepSizeReductionFactor )
63 ctr1 = ctr1 + 1
64
65
66# simulation settings
67tEnd = 5.0 # s # simulate system until tend
68
69# actual system definition starts here
70for i in range( len(stepSizeList) ):
71
72 SC = exu.SystemContainer()
73 mbs = SC.AddSystem()
74
75 simulationSettings = exu.SimulationSettings() #takes currently set values or default values
76 simulationSettings.solution.file.name = 'solution/BASpringDamperSystem h=' + str(stepSizeList[i]) + '.txt'
77
78 print(int( tEnd/stepSizeList[i] ) )
79
80 writeStep = 0.01
81 simulationSettings.timeIntegration.numberOfSteps = int( tEnd/stepSizeList[i] )
82 simulationSettings.timeIntegration.endTime = tEnd
83 simulationSettings.solution.file.write = createSolutionFile
84 simulationSettings.solution.file.writePeriod = writeStep
85 simulationSettings.solution.precision = 16
86 simulationSettings.show.computationTime = True
87 simulationSettings.timeIntegration.verboseMode = 1
88 #simulationSettings.linearSolver.solverType = exu.LinearSolverType.EigenSparse
89 #simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1
90 #simulationSettings.show.statistics = True
91
92
93 SC.visualizationSettings.bodies.showNumbers = True
94 SC.visualizationSettings.nodes.defaultSize = 0.1
95 SC.visualizationSettings.markers.defaultSize = 0.05
96 SC.visualizationSettings.connectors.defaultSize = 0.1
97
98
99
100####################### define nodes #########################################
101 nodeCounter = 0 # zero based
102 for i in range(numberOfNodesY):
103 for j in range(numberOfNodesX):
104
105 nodeName = "node " + str(nodeCounter)
106
107 # increase node counter
108 nodeCounter = nodeCounter + 1
109
110 if j == 0:
111 nodeDict = {"nodeType": "PointGround",
112 "referenceCoordinates": [0.0, i*springDamperElementLength, 0.0],
113 "name": nodeName}
114
115 elif nodeCounter == (i+1)*numberOfNodesX and clambRightSide == True:
116 nodeDict = {"nodeType": "PointGround",
117 "referenceCoordinates": [j*springDamperElementLength, i*springDamperElementLength, 0.0],
118 "name": nodeName}
119
120 else:
121 nodeDict = {"nodeType": "Point",
122 "referenceCoordinates": [j*springDamperElementLength, i*springDamperElementLength, 0.0],
123 "initialCoordinates": [0.0, 0.0, 0.0],
124 "name": nodeName}
125
126 # add node to mbs
127 mbs.AddNode(nodeDict)
128
129
130####################### create objects #######################################
131 for i in range(numberOfNodes):
132
133 nodeDict = mbs.GetNode(i)
134 nodeName = nodeDict["name"]
135 nodeNumber = mbs.GetNodeNumber(nodeName)
136
137 massPointName = "mass point - " + nodeName
138
139 objectDict = {"objectType": "MassPoint",
140 "mass": massMassPoint,
141 "nodeNumber": nodeNumber,
142 "name": massPointName}
143
144 mbs.AddObject(objectDict)
145
146
147####################### add markers ##########################################
148 bodyMassMarkerName = list()
149 nodePositionMarkerName = list()
150 for i in range(numberOfNodes):
151
152 nodeDict = mbs.GetNode(i)
153 nodeName = nodeDict["name"]
154 nodeNumber = mbs.GetNodeNumber(nodeName)
155
156 bodyMassMarkerName.append("body mass marker - " + nodeName)
157 nodePositionMarkerName.append("node position marker - " + nodeName)
158
159 # body mass - used for garvitational load
160 bodyMassMarkerDict = {"markerType": "BodyMass",
161 "bodyNumber": exu.ObjectIndex(nodeNumber), #nodeNumber=bodyNumber
162 "name": bodyMassMarkerName[i]}
163
164 nodePositionMarkerDict = {"markerType": "NodePosition",
165 "nodeNumber": nodeNumber,
166 "name": nodePositionMarkerName[i]}
167
168 mbs.AddMarker(bodyMassMarkerDict)
169 mbs.AddMarker(nodePositionMarkerDict)
170
171
172####################### add loads ############################################
173 for i in range( len(bodyMassMarkerName) ):
174
175 markerNumberOfBodyMassMarker = mbs.GetMarkerNumber(bodyMassMarkerName[i])
176 loadName = "gravity " + str(bodyMassMarkerName[i])
177
178 loadDict = {"loadType": "MassProportional",
179 "markerNumber": markerNumberOfBodyMassMarker,
180 "loadVector": [0.0, 0.0, -9.81],
181 "name": loadName}
182
183 mbs.AddLoad(loadDict)
184
185
186####################### add connectors #######################################
187 ctr1 = 0
188 elementCtr = 0;
189
190 # spring-damper elements in x-direction
191 if useSpringDamperElementsX == True:
192
193 for i in range(numberOfNodesY):
194
195 for j in range(numberOfNodesX-1):
196
197 markerNumberOfPositionMarkerL = mbs.GetMarkerNumber(nodePositionMarkerName[ctr1])
198 markerNumberOfPositionMarkerR = mbs.GetMarkerNumber(nodePositionMarkerName[ctr1+1])
199
200 ctr1 = ctr1 + 1
201
202 springDamperElementName = "spring damper" + str(elementCtr)
203
204 objectDict = {"objectType": "ConnectorSpringDamper",
205 "markerNumbers": [markerNumberOfPositionMarkerL, markerNumberOfPositionMarkerR],
206 "stiffness": elementStiffness,
207 "damping": elementDamping,
208 "force": 0,
209 "referenceLength": springDamperElementLength,
210 "activeConnector": activateElement,
211 "name": springDamperElementName,
212 "VdrawSize": 0.0}
213
214 mbs.AddObject(objectDict)
215
216 elementCtr = elementCtr + 1
217
218 ctr1 = ctr1 + 1
219
220 # spring-damper elements in y-direction
221 if useSpringDamperElementsY == True:
222
223 ctr1 = 0
224 for i in range(numberOfNodesY-1):
225
226 for j in range(numberOfNodesX):
227
228 markerNumberOfPositionMarkerL = mbs.GetMarkerNumber(nodePositionMarkerName[ctr1])
229 markerNumberOfPositionMarkerR = mbs.GetMarkerNumber(nodePositionMarkerName[ctr1+numberOfNodesX])
230
231 ctr1 = ctr1 + 1
232
233 springDamperElementName = "spring damper" + str(elementCtr)
234
235 objectDict = {"objectType": "ConnectorSpringDamper",
236 "markerNumbers": [markerNumberOfPositionMarkerL, markerNumberOfPositionMarkerR],
237 "stiffness": elementStiffness,
238 "damping": elementDamping,
239 "force": 0,
240 "referenceLength": springDamperElementLength,
241 "activeConnector": activateElement,
242 "name": springDamperElementName,
243 "VdrawSize": 0.0}
244
245 mbs.AddObject(objectDict)
246
247 elementCtr = elementCtr + 1
248
249
250 # spring-damper elements in off-diagnoal-direction
251 if useSpringDamperElementsOffDiagonal == True:
252
253 ctr1 = 0
254 for i in range(numberOfNodesY-1): #range(numberOfNodes):
255
256 for j in range(numberOfNodesX-1):
257
258 markerNumberOfPositionMarkerL = mbs.GetMarkerNumber(nodePositionMarkerName[ctr1])
259 markerNumberOfPositionMarkerR = mbs.GetMarkerNumber(nodePositionMarkerName[ctr1+numberOfNodesX+1])
260
261 ctr1 = ctr1 + 1
262
263 springDamperElementName = "spring damper" + str(elementCtr)
264
265 objectDict = {"objectType": "ConnectorSpringDamper",
266 "markerNumbers": [markerNumberOfPositionMarkerL, markerNumberOfPositionMarkerR],
267 "stiffness": elementStiffness,
268 "damping": elementDamping,
269 "force": 0,
270 "referenceLength": springDamperDiagnoalElementLength,
271 "activeConnector": activateElement,
272 "name": springDamperElementName,
273 "VdrawSize": 0.0}
274
275 mbs.AddObject(objectDict)
276
277 elementCtr = elementCtr + 1
278
279 ctr1 = ctr1 + 1
280
281
282 # spring-damper elements in diagnoal-direction
283 if useSpringDamperElementsDiagonal == True:
284
285 ctr1 = 0
286 for i in range(numberOfNodesY-1):
287
288 for j in range(numberOfNodesX-1):
289
290 markerNumberOfPositionMarkerL = mbs.GetMarkerNumber(nodePositionMarkerName[ctr1+1])
291 markerNumberOfPositionMarkerR = mbs.GetMarkerNumber(nodePositionMarkerName[ctr1+numberOfNodesX])
292
293 ctr1 = ctr1 + 1
294
295 springDamperElementName = "spring damper" + str(elementCtr)
296
297 objectDict = {"objectType": "ConnectorSpringDamper",
298 "markerNumbers": [markerNumberOfPositionMarkerL, markerNumberOfPositionMarkerR],
299 "stiffness": elementStiffness,
300 "damping": elementDamping,
301 "force": 0,
302 "referenceLength": springDamperDiagnoalElementLength,
303 "activeConnector": activateElement,
304 "name": springDamperElementName,
305 "VdrawSize": 0.0}
306
307 mbs.AddObject(objectDict)
308
309 elementCtr = elementCtr + 1
310
311 ctr1 = ctr1 + 1
312
313
314
315####################### assemble and solve mbs ###############################
316 mbs.Assemble()
317 #print(mbs)
318
319
320 SC.renderer.Start()
321 mbs.SolveDynamic(simulationSettings, solverType = exudyn.DynamicSolverType.ODE23)
322 SC.renderer.DoIdleTasks()
323 SC.renderer.Stop() #safely close rendering window!