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.itemInterface import *
15from exudyn.graphicsDataUtilities import *
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 = [GraphicsDataQuad([[-L0,hy0,z],[ L,hy0,z],[ L,hy1,z],[-L0,hy1,z]],
152 color=color4lightgrey)]
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':color4grey}
159 line1 = {'type':'Line', 'data':[ox+0,oy-t,0, ox+0,oy+t,0], 'color':color4grey}
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 = color4steelblue
172 if i==0:
173 col = color4green
174 elif i==N-1:
175 col = color4lightred
176
177 gSphere = GraphicsDataSphere(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, physicsMass=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, physicsMass=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(physicsMass = 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.general.textSize = 12
254SC.visualizationSettings.openGL.lineWidth = 2
255SC.visualizationSettings.openGL.multiSampling = 4
256SC.visualizationSettings.general.graphicsUpdateInterval = 0.005
257#SC.visualizationSettings.window.renderWindowSize=[1024,900]
258SC.visualizationSettings.window.renderWindowSize=[1600,1000]
259SC.visualizationSettings.general.showSolverInformation = False
260SC.visualizationSettings.general.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
269exu.StartRenderer()
270
271#++++++++++++++++++++++++++++++++++++++++
272#setup simulation settings and run interactive dialog:
273simulationSettings = exu.SimulationSettings()
274simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1
275simulationSettings.solutionSettings.writeSolutionToFile = False
276simulationSettings.solutionSettings.solutionWritePeriod = 0.1 #data not used
277simulationSettings.solutionSettings.sensorsWritePeriod = 0.1 #data not used
278simulationSettings.solutionSettings.solutionInformation = '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.displayComputationTime = 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
353dialog = InteractiveDialog(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
361exu.StopRenderer() #safely close rendering window!