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