nMassOscillator.py
You can view and download this file on Github: nMassOscillator.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, offset, 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 = 10
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
154# for 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
269
270#++++++++++++++++++++++++++++++++++++++++
271#setup simulation settings and run interactive dialog:
272simulationSettings = exu.SimulationSettings()
273simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1
274simulationSettings.solutionSettings.writeSolutionToFile = True
275simulationSettings.solutionSettings.solutionWritePeriod = 0.05 #data not used
276simulationSettings.solutionSettings.sensorsWritePeriod = 0.1 #data not used
277simulationSettings.solutionSettings.solutionInformation = 'n-mass-oscillatior'
278simulationSettings.timeIntegration.verboseMode = 1 #turn off, because of lots of output
279
280simulationSettings.timeIntegration.numberOfSteps = int(1000)
281simulationSettings.timeIntegration.endTime = 5
282simulationSettings.timeIntegration.newton.useModifiedNewton = True
283
284simulationSettings.displayComputationTime = True
285# simulationSettings.numberOfThreads = 1
286
287#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
288#set up interactive window
289
290
291mbs.SolveDynamic(simulationSettings=simulationSettings)
292
293mbs.SolutionViewer()