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