nMassOscillatorInteractive.py
You can view and download this file on Github: nMassOscillatorInteractive.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: Nonlinear oscillations interactive simulation
5#
6# Author: Johannes Gerstmayr
7# Date: 2020-01-16
8#
9# 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.
10#
11#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
12
13import exudyn as exu
14from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
15import exudyn.graphics as graphics
16import matplotlib.pyplot as plt
17from exudyn.interactive import InteractiveDialog
18
19import numpy as np
20from math import sin, cos, pi, sqrt
21
22import time #for sleep()
23SC = exu.SystemContainer()
24mbs = SC.AddSystem()
25
26#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
27#
28N = 12; #number of masses
29spring = 800; #stiffness [800]
30mass = 1; #mass
31damper = 2; #old:0.1; damping parameter
32force = 1; #force amplitude
33
34d0 = damper*spring/(2*sqrt(mass*spring)) #dimensionless damping for single mass
35
36
37caseHarmonic = 1
38#damper=2
39#mode1:force=0.52
40#mode2:force=2
41#mode3:force=4
42#mode12: force=100, damping=0.05, period=0.005
43
44caseStep = 2
45#damper=1
46#force=20
47
48
49mbs.variables['loadCase'] = caseHarmonic
50mbs.variables['resetMotion'] = 0 #run
51mbs.variables['forceAmplitude'] = force
52eigenMode = 1
53
54h = 0.002 #step size
55deltaT = 0.01 #time period to be simulated between every update
56
57
58
59# if (mode < 2) h=h*2; F=0.4*F; end
60# if (mode < 5) h=h*2; F=0.4*F; end
61# if (mode < 6) h=h*2.5; end
62# if (mode == 6) F=F*2; end
63# if (mode == 3) h=h*0.5; end
64
65# if (mode > 16) F=2*F; end
66# if (mode > 10) F=2*F; end
67# if (N < 11) h=h/2; l_mass = 2*l_mass; r_mass=2*r_mass; F=0.5*F; end
68# if (N < 6) h=h/2; l_mass = 2*l_mass; r_mass=2*r_mass; end
69
70# om=sqrt(diag(ew))
71
72
73
74
75
76#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
77#create model for linear and nonlinear oscillator:
78# L=0.5
79# load0 = 80
80
81# omega0=np.sqrt(spring/mass)
82# f0 = 0.*omega0/(2*np.pi)
83# f1 = 1.*omega0/(2*np.pi)
84
85# tEnd = 200 #end time of simulation
86# steps = 20000 #number of steps
87
88omegaInit = 3.55
89omegaMax = 40 #for plots
90mbs.variables['mode'] = 0 #0=linear, 1=cubic nonlinear
91mbs.variables['omega'] = omegaInit #excitation frequency changed by user
92#mbs.variables['omega'] = omegaInit #excitation, changed in simFunction
93mbs.variables['phi'] = 0 #excitation phase, used to get smooth excitations
94mbs.variables['stiffness'] = spring
95mbs.variables['damping'] = damper
96mbs.variables['dampingPrev'] = damper
97
98
99# #user function for spring force
100# def springForce(mbs, t, itemIndex, u, v, k, d): #changed 2023-01-21:, mu, muPropZone):
101# k=mbs.variables['stiffness']
102# d=mbs.variables['damping']
103# if mbs.variables['mode'] == 0:
104# return k*u + v*d
105# else:
106# #return 0.1*k*u+k*u**3+v*d
107# return k*u+1000*k*u**3+v*d #breaks down at 13.40Hz
108
109# mode = 0 #0...forward, 1...backward
110
111#user function for load
112def userLoad(mbs, t, load):
113 f = mbs.variables['forceAmplitude']
114 fact = 1
115 if mbs.variables['loadCase']==caseHarmonic:
116 fact = sin(mbs.GetSensorValues(mbs.variables['sensorPhi']))
117 return f*fact
118
119def userLoad3D(mbs,t, load):
120 f = mbs.variables['forceAmplitude']
121 fact = 1
122 if mbs.variables['loadCase']==caseHarmonic:
123 fact = sin(mbs.GetSensorValues(mbs.variables['sensorPhi']))
124 mbs.SetLoadParameter(0,'loadVector',[fact,0,0])
125 return [f*fact,0,0]
126
127#dummy user function for frequency
128def userFrequency(mbs, t, load):
129 return mbs.variables['omega']
130
131#user function used in GenericODE2 to integrate current omega
132def UFintegrateOmega(mbs, t, itemIndex, q, q_t):
133 return [mbs.variables['omega']] #current frequency*2*pi is integrated into phi, return vector!
134
135#ground node
136nGround=mbs.AddNode(NodePointGround(referenceCoordinates = [0,0,0]))
137
138#drawing parameters:
139l_mass = 0.2 #spring length
140r_mass = 0.030*2 #radius of mass
141r_spring = r_mass*1.2
142L0 = l_mass*1
143L = N * l_mass + 4*l_mass
144z=-r_mass-0.1
145hy=0.25*L
146hy1=2*hy - 4*r_mass
147hy0=-4*r_mass
148maxAmp0 = 0.1
149maxAmpN = 0.1*N
150
151background = [graphics.Quad([[-L0,hy0,z],[ L,hy0,z],[ L,hy1,z],[-L0,hy1,z]],
152 color=graphics.color.lightgrey)]
153offCircleY = 1*hy
154for i in range(N):
155 t=r_mass*0.5
156 ox = l_mass*(i+1)
157 oy = offCircleY
158 line0 = {'type':'Line', 'data':[ox-t,oy+0,0, ox+t,oy+0,0], 'color':graphics.color.grey}
159 line1 = {'type':'Line', 'data':[ox+0,oy-t,0, ox+0,oy+t,0], 'color':graphics.color.grey}
160 background += [line0, line1]
161
162oGround=mbs.AddObject(ObjectGround(visualization=VObjectGround(graphicsData=background)))
163#marker for ground (=fixed):
164groundMarker=mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= nGround, coordinate = 0))
165prevMarker = groundMarker
166nMass = []
167mbs.variables['springDamperList'] = []
168
169for i in range(N):
170 #node for 3D mass point:
171 col = graphics.color.steelblue
172 if i==0:
173 col = graphics.color.green
174 elif i==N-1:
175 col = graphics.color.lightred
176
177 gSphere = graphics.Sphere(point=[0,0,0], radius=r_mass, color=col, nTiles=16)
178 node = mbs.AddNode(Node1D(referenceCoordinates = [l_mass*(1+len(nMass))],
179 initialCoordinates=[0.],
180 initialVelocities=[0.]))
181 massPoint = mbs.AddObject(Mass1D(nodeNumber = node, mass=mass,
182 referencePosition=[0,0,0],
183 visualization=VMass1D(graphicsData=[gSphere])))
184
185 gCircle = {'type':'Circle','position':[0,0,0],'radius':0.5*r_mass, 'color':col}
186 massPoint2 = mbs.AddObject(Mass1D(nodeNumber = node, mass=0,
187 referencePosition=[l_mass*(len(nMass)+1),offCircleY-l_mass*(len(nMass)+1),0],
188 referenceRotation=[[0,1,0],[1,0,0],[0,0,1]],
189 visualization=VMass1D(graphicsData=[gCircle])))
190
191
192 # node = mbs.AddNode(Point(referenceCoordinates = [l_mass*(1+len(nMass)),0,0]))
193
194 # massPoint = mbs.AddObject(MassPoint(mass = mass, nodeNumber = node,
195 # visualization=VMassPoint(graphicsData=[gSphere])))
196
197 nMass += [node]
198 #marker for springDamper for first (x-)coordinate:
199 nodeMarker =mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= node, coordinate = 0))
200
201 #Spring-Damper between two marker coordinates
202 sd = mbs.AddObject(CoordinateSpringDamper(markerNumbers = [prevMarker, nodeMarker],
203 stiffness = spring, damping = damper,
204 #springForceUserFunction = springForce,
205 visualization=VCoordinateSpringDamper(drawSize=r_spring)))
206 mbs.variables['springDamperList'] += [sd]
207 prevMarker = nodeMarker
208
209#add load to last mass:
210if False: #scalar load
211 mbs.AddLoad(LoadCoordinate(markerNumber = nodeMarker,
212 load = 0, loadUserFunction=userLoad)) #load set in user function
213else:
214 mMassN = mbs.AddMarker(MarkerBodyPosition(bodyNumber= massPoint, localPosition=[0,0,0]))
215 mbs.AddLoad(Force(markerNumber=mMassN, loadVector=[1,0,0],
216 loadVectorUserFunction=userLoad3D))
217
218# #dummy load applied to ground marker, just to record/integrate frequency
219lFreq = mbs.AddLoad(LoadCoordinate(markerNumber = groundMarker,
220 load = 0, loadUserFunction=userFrequency))
221
222sensPos0 = mbs.AddSensor(SensorNode(nodeNumber=nMass[0], fileName='solution/nMassPos0.txt',
223 outputVariableType=exu.OutputVariableType.Coordinates))
224sensPosN = mbs.AddSensor(SensorNode(nodeNumber=nMass[-1], fileName='solution/nMassPosN.txt',
225 outputVariableType=exu.OutputVariableType.Coordinates))
226sensFreq = mbs.AddSensor(SensorLoad(loadNumber=lFreq, fileName='solution/nMassFreq.txt',
227 visualization=VSensorLoad(show=False)))
228
229#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
230#compute eigenvalues
231from exudyn.solver import ComputeODE2Eigenvalues
232mbs.Assemble()
233[values, vectors] = ComputeODE2Eigenvalues(mbs)
234print('omegas (rad/s)=', np.sqrt(values))
235
236#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
237#integrate omega: node used to integrate omega into phi for excitation function
238nODE2=mbs.AddNode(NodeGenericODE2(referenceCoordinates=[0], initialCoordinates=[0],initialCoordinates_t=[0],
239 numberOfODE2Coordinates=1))
240
241oODE2=mbs.AddObject(ObjectGenericODE2(nodeNumbers=[nODE2],massMatrix=np.eye(1),
242 forceUserFunction=UFintegrateOmega,
243 visualization=VObjectGenericODE2(show=False)))
244#improved version, using integration of omega:
245mbs.variables['sensorPhi'] = mbs.AddSensor(SensorNode(nodeNumber=nODE2, fileName='solution/nonlinearPhi.txt',
246 outputVariableType = exu.OutputVariableType.Coordinates_t,
247 visualization=VSensorNode(show=False)))
248#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
249#finalize model and settings
250mbs.Assemble()
251
252
253SC.visualizationSettings.view0.window.globalFontSize = 12
254SC.visualizationSettings.openGL.lineWidth = 2
255SC.visualizationSettings.openGL.multiSampling = 4
256SC.visualizationSettings.general.graphicsUpdateInterval = 0.005
257#SC.visualizationSettings.view0.window.renderWindowSize=[1024,900]
258SC.visualizationSettings.view0.window.renderWindowSize=[1600,1000]
259SC.visualizationSettings.general.showSolverInformation = False
260SC.visualizationSettings.view0.scene.drawCoordinateSystem = False
261
262SC.visualizationSettings.loads.fixedLoadSize=0
263SC.visualizationSettings.loads.loadSizeFactor=0.5
264SC.visualizationSettings.loads.drawSimplified=False
265SC.visualizationSettings.loads.defaultSize=1
266SC.visualizationSettings.loads.defaultRadius=0.01
267
268SC.visualizationSettings.general.autoFitScene = True #otherwise, renderState not accepted for zoom
269SC.renderer.Start()
270
271#++++++++++++++++++++++++++++++++++++++++
272#setup simulation settings and run interactive dialog:
273simulationSettings = exu.SimulationSettings()
274simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1
275simulationSettings.solution.file.write = False
276simulationSettings.solution.file.writePeriod = 0.1 #data not used
277simulationSettings.solution.sensors.writePeriod = 0.1 #data not used
278simulationSettings.solution.file.information = 'n-mass-oscillatior'
279simulationSettings.timeIntegration.verboseMode = 0 #turn off, because of lots of output
280
281simulationSettings.timeIntegration.numberOfSteps = int(deltaT/h)
282simulationSettings.timeIntegration.endTime = deltaT
283simulationSettings.timeIntegration.newton.useModifiedNewton = True
284
285simulationSettings.show.computationTime = True
286simulationSettings.parallel.numberOfThreads = 2
287
288#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
289#set up interactive window
290
291#interactive function
292def SimulationUF(mbs, dialog):
293 if mbs.variables['resetMotion']:
294 u = mbs.systemData.GetODE2Coordinates()
295 u *= 0
296 mbs.systemData.SetODE2Coordinates(u, configuration=exu.ConfigurationType.Initial)
297 v = mbs.systemData.GetODE2Coordinates_t()
298 v *= 0
299 mbs.systemData.SetODE2Coordinates_t(v, configuration=exu.ConfigurationType.Initial)
300 #clear plot data:
301 nPoints = dialog.plots['nPoints']
302 t = mbs.systemData.GetTime()
303 p = mbs.variables['period']
304 for j in range(len(dialog.plots['sensors'])):
305 #dialog.plots['data'][j][:,1] *= 0
306 # tEnd = dialog.plots['data'][j][:,0]
307 dialog.plots['data'][j] = np.zeros((nPoints,2))
308 dialog.plots['data'][j][:,0] = p*np.arange(0,nPoints) + ( t-2*p*(nPoints))
309 dialog.plots['currentIndex'] = nPoints-1
310
311 if mbs.variables['dampingPrev'] != mbs.variables['damping']:
312 #print('here')
313 for sd in mbs.variables['springDamperList']:
314 mbs.SetObjectParameter(sd,'damping',mbs.variables['damping'])
315 mbs.variables['dampingPrev'] = mbs.variables['damping']
316
317 dialog.period = mbs.variables['period']
318
319#++++++++++++++++++++++++++++
320#define items for dialog
321dialogItems = [{'type':'label', 'text':'n-mass oscillatior simulator', 'grid':(0,0,2), 'options':['L']},
322 {'type':'radio', 'textValueList':[('harmonic',caseHarmonic),('step',caseStep)], 'value': mbs.variables['loadCase'], 'variable':'loadCase', 'grid': [(2,0),(2,1)]},
323 {'type':'label', 'text':'excitation frequency (rad/s):', 'grid':(5,0)},
324 {'type':'slider', 'range':(0, 60), 'value':omegaInit, 'steps':601, 'variable':'omega', 'grid':(5,1)},
325 {'type':'radio', 'textValueList':[('run',0),('reset motion',1)], 'value': mbs.variables['resetMotion'], 'variable':'resetMotion', 'grid': [(6,0),(6,1)]},
326 {'type':'label', 'text':'force amplitude:', 'grid':(7,0)},
327 {'type':'slider', 'range': (0, 100), 'value':mbs.variables['forceAmplitude'], 'steps':1001, 'variable':'forceAmplitude', 'grid':(7,1)},
328 {'type':'label', 'text':'damping:', 'grid':(8,0)},
329 {'type':'slider', 'range': (0, 40), 'value':damper, 'steps':801, 'variable':'damping', 'grid':(8,1)},
330 {'type':'label', 'text':'period:', 'grid':(9,0)},
331 {'type':'slider', 'range': (0.001, 0.1), 'value':deltaT, 'steps':100, 'variable':'period', 'grid':(9,1)},
332 # {'type':'label', 'text':'stiffness:', 'grid':(7,0)},
333 # {'type':'slider', 'range':(0, 1000), 'value':spring, 'steps':500, 'variable':'stiffness', 'grid':(7,1)}
334 ]
335
336#++++++++++++++++++++++++++++++++++++++++
337#specify subplots to be shown interactively
338plt.close('all')
339
340plots={'fontSize':12,'sizeInches':(16,12),'nPoints':400,
341 'subplots':(2,1),
342 'sensors':[[(sensPos0,0),(sensPos0,1),'time','mass0 position'],
343 [(sensPosN,0),(sensPosN,1),'time','massN position'],
344 #[(sensFreq,0),(sensFreq,1),'time','excitation freq.(rad/s)']
345 ],
346 'subplots':False,
347 'lineStyles':['g-','r-'],
348 'limitsX':[(),()], #omit if time auto-range
349 'limitsY':[(),()]}#,(0,omegaMax)]}
350 # 'limitsY':[(-maxAmp0,maxAmp0),(-maxAmpN,maxAmpN)]}#,(0,omegaMax)]}
351
352
353InteractiveDialog(mbs=mbs, simulationSettings=simulationSettings,
354 simulationFunction=SimulationUF, title='Interactive window',
355 dialogItems=dialogItems, period=deltaT,
356 #realtimeFactor=1,
357 plots=plots,
358 fontSize=12)
359
360# #stop solver and close render window
361SC.renderer.Stop() #safely close rendering window!