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.utilities import * #includes itemInterface and rigidBodyUtilities
16import exudyn.graphics as graphics #only import if it does not conflict
17
18import matplotlib.pyplot as plt
19from exudyn.interactive import InteractiveDialog
20
21import numpy as np
22from math import sin, cos, pi, sqrt
23
24import time #for sleep()
25SC = exu.SystemContainer()
26mbs = SC.AddSystem()
27
28#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
29#
30N = 12; #number of masses
31spring = 800; #stiffness [800]
32mass = 1; #mass
33damper = 2; #old:0.1; damping parameter
34force = 1; #force amplitude
35
36d0 = damper*spring/(2*sqrt(mass*spring)) #dimensionless damping for single mass
37
38
39caseHarmonic = 1
40#damper=2
41#mode1:force=0.52
42#mode2:force=2
43#mode3:force=4
44#mode12: force=100, damping=0.05, period=0.005
45
46caseStep = 2
47#damper=1
48#force=20
49
50
51mbs.variables['loadCase'] = caseHarmonic
52mbs.variables['resetMotion'] = 0 #run
53mbs.variables['forceAmplitude'] = force
54eigenMode = 1
55
56h = 0.002 #step size
57deltaT = 0.01 #time period to be simulated between every update
58
59
60
61# if (mode < 2) h=h*2; F=0.4*F; end
62# if (mode < 5) h=h*2; F=0.4*F; end
63# if (mode < 6) h=h*2.5; end
64# if (mode == 6) F=F*2; end
65# if (mode == 3) h=h*0.5; end
66
67# if (mode > 16) F=2*F; end
68# if (mode > 10) F=2*F; end
69# if (N < 11) h=h/2; l_mass = 2*l_mass; r_mass=2*r_mass; F=0.5*F; end
70# if (N < 6) h=h/2; l_mass = 2*l_mass; r_mass=2*r_mass; end
71
72# om=sqrt(diag(ew))
73
74
75
76
77
78#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
79#create model for linear and nonlinear oscillator:
80# L=0.5
81# load0 = 80
82
83# omega0=np.sqrt(spring/mass)
84# f0 = 0.*omega0/(2*np.pi)
85# f1 = 1.*omega0/(2*np.pi)
86
87# tEnd = 200 #end time of simulation
88# steps = 20000 #number of steps
89
90omegaInit = 3.55
91omegaMax = 40 #for plots
92mbs.variables['mode'] = 0 #0=linear, 1=cubic nonlinear
93mbs.variables['omega'] = omegaInit #excitation frequency changed by user
94#mbs.variables['omega'] = omegaInit #excitation, changed in simFunction
95mbs.variables['phi'] = 0 #excitation phase, used to get smooth excitations
96mbs.variables['stiffness'] = spring
97mbs.variables['damping'] = damper
98mbs.variables['dampingPrev'] = damper
99
100
101# #user function for spring force
102# def springForce(mbs, t, itemIndex, u, v, k, d, offset, mu, muPropZone):
103# k=mbs.variables['stiffness']
104# d=mbs.variables['damping']
105# if mbs.variables['mode'] == 0:
106# return k*u + v*d
107# else:
108# #return 0.1*k*u+k*u**3+v*d
109# return k*u+1000*k*u**3+v*d #breaks down at 13.40Hz
110
111# mode = 0 #0...forward, 1...backward
112
113#user function for load
114def userLoad(mbs, t, load):
115 f = mbs.variables['forceAmplitude']
116 fact = 1
117 if mbs.variables['loadCase']==caseHarmonic:
118 fact = sin(mbs.GetSensorValues(mbs.variables['sensorPhi']))
119 return f*fact
120
121def userLoad3D(mbs,t, load):
122 f = mbs.variables['forceAmplitude']
123 fact = 10
124 # if mbs.variables['loadCase']==caseHarmonic:
125 # fact = sin(mbs.GetSensorValues(mbs.variables['sensorPhi']))
126 # mbs.SetLoadParameter(0,'loadVector',[fact,0,0])
127 return [f*fact,0,0]
128
129#dummy user function for frequency
130def userFrequency(mbs, t, load):
131 return mbs.variables['omega']
132
133#user function used in GenericODE2 to integrate current omega
134def UFintegrateOmega(mbs, t, itemIndex, q, q_t):
135 return [mbs.variables['omega']] #current frequency*2*pi is integrated into phi, return vector!
136
137#ground node
138nGround=mbs.AddNode(NodePointGround(referenceCoordinates = [0,0,0]))
139
140#drawing parameters:
141l_mass = 0.2 #spring length
142r_mass = 0.030*2 #radius of mass
143r_spring = r_mass*1.2
144L0 = l_mass*1
145L = N * l_mass + 4*l_mass
146z=-r_mass-0.1
147hy=0.25*L
148hy1=2*hy - 4*r_mass
149hy0=-4*r_mass
150maxAmp0 = 0.1
151maxAmpN = 0.1*N
152
153background = [graphics.Quad([[-L0,hy0,z],[ L,hy0,z],[ L,hy1,z],[-L0,hy1,z]],
154 color=graphics.color.lightgrey)]
155offCircleY = 1*hy
156# for i in range(N):
157# t=r_mass*0.5
158# ox = l_mass*(i+1)
159# oy = offCircleY
160# line0 = {'type':'Line', 'data':[ox-t,oy+0,0, ox+t,oy+0,0], 'color':graphics.color.grey}
161# line1 = {'type':'Line', 'data':[ox+0,oy-t,0, ox+0,oy+t,0], 'color':graphics.color.grey}
162# background += [line0, line1]
163
164oGround=mbs.AddObject(ObjectGround(visualization=VObjectGround(graphicsData=background)))
165#marker for ground (=fixed):
166groundMarker=mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= nGround, coordinate = 0))
167prevMarker = groundMarker
168nMass = []
169mbs.variables['springDamperList'] = []
170
171for i in range(N):
172 #node for 3D mass point:
173 col = graphics.color.steelblue
174 if i==0:
175 col = graphics.color.green
176 elif i==N-1:
177 col = graphics.color.lightred
178
179 gSphere = graphics.Sphere(point=[0,0,0], radius=r_mass, color=col, nTiles=16)
180 node = mbs.AddNode(Node1D(referenceCoordinates = [l_mass*(1+len(nMass))],
181 initialCoordinates=[0.],
182 initialVelocities=[0.]))
183 massPoint = mbs.AddObject(Mass1D(nodeNumber = node, physicsMass=mass,
184 referencePosition=[0,0,0],
185 visualization=VMass1D(graphicsData=[gSphere])))
186
187 # gCircle = {'type':'Circle','position':[0,0,0],'radius':0.5*r_mass, 'color':col}
188 # massPoint2 = mbs.AddObject(Mass1D(nodeNumber = node, physicsMass=0,
189 # referencePosition=[l_mass*(len(nMass)+1),offCircleY-l_mass*(len(nMass)+1),0],
190 # referenceRotation=[[0,1,0],[1,0,0],[0,0,1]],
191 # visualization=VMass1D(graphicsData=[gCircle])))
192
193
194 # node = mbs.AddNode(Point(referenceCoordinates = [l_mass*(1+len(nMass)),0,0]))
195
196 # massPoint = mbs.AddObject(MassPoint(physicsMass = mass, nodeNumber = node,
197 # visualization=VMassPoint(graphicsData=[gSphere])))
198
199 nMass += [node]
200 #marker for springDamper for first (x-)coordinate:
201 nodeMarker =mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= node, coordinate = 0))
202
203 #Spring-Damper between two marker coordinates
204 sd = mbs.AddObject(CoordinateSpringDamper(markerNumbers = [prevMarker, nodeMarker],
205 stiffness = spring, damping = damper,
206 #springForceUserFunction = springForce,
207 visualization=VCoordinateSpringDamper(drawSize=r_spring)))
208 mbs.variables['springDamperList'] += [sd]
209 prevMarker = nodeMarker
210
211#add load to last mass:
212if False: #scalar load
213 mbs.AddLoad(LoadCoordinate(markerNumber = nodeMarker,
214 load = 0, loadUserFunction=userLoad)) #load set in user function
215else:
216 mMassN = mbs.AddMarker(MarkerBodyPosition(bodyNumber= massPoint, localPosition=[0,0,0]))
217 mbs.AddLoad(Force(markerNumber=mMassN, loadVector=[1,0,0],
218 loadVectorUserFunction=userLoad3D))
219
220# #dummy load applied to ground marker, just to record/integrate frequency
221lFreq = mbs.AddLoad(LoadCoordinate(markerNumber = groundMarker,
222 load = 0, loadUserFunction=userFrequency))
223
224sensPos0 = mbs.AddSensor(SensorNode(nodeNumber=nMass[0], fileName='solution/nMassPos0.txt',
225 outputVariableType=exu.OutputVariableType.Coordinates))
226sensPosN = mbs.AddSensor(SensorNode(nodeNumber=nMass[-1], fileName='solution/nMassPosN.txt',
227 outputVariableType=exu.OutputVariableType.Coordinates))
228sensFreq = mbs.AddSensor(SensorLoad(loadNumber=lFreq, fileName='solution/nMassFreq.txt',
229 visualization=VSensorLoad(show=False)))
230
231#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
232#compute eigenvalues
233from exudyn.solver import ComputeODE2Eigenvalues
234mbs.Assemble()
235[values, vectors] = ComputeODE2Eigenvalues(mbs)
236print('omegas (rad/s)=', np.sqrt(values))
237
238#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
239#integrate omega: node used to integrate omega into phi for excitation function
240nODE2=mbs.AddNode(NodeGenericODE2(referenceCoordinates=[0], initialCoordinates=[0],initialCoordinates_t=[0],
241 numberOfODE2Coordinates=1))
242
243oODE2=mbs.AddObject(ObjectGenericODE2(nodeNumbers=[nODE2],massMatrix=np.eye(1),
244 forceUserFunction=UFintegrateOmega,
245 visualization=VObjectGenericODE2(show=False)))
246#improved version, using integration of omega:
247mbs.variables['sensorPhi'] = mbs.AddSensor(SensorNode(nodeNumber=nODE2, fileName='solution/nonlinearPhi.txt',
248 outputVariableType = exu.OutputVariableType.Coordinates_t,
249 visualization=VSensorNode(show=False)))
250#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
251#finalize model and settings
252mbs.Assemble()
253
254
255SC.visualizationSettings.general.textSize = 12
256SC.visualizationSettings.openGL.lineWidth = 2
257SC.visualizationSettings.openGL.multiSampling = 4
258SC.visualizationSettings.general.graphicsUpdateInterval = 0.005
259#SC.visualizationSettings.window.renderWindowSize=[1024,900]
260SC.visualizationSettings.window.renderWindowSize=[1600,1000]
261SC.visualizationSettings.general.showSolverInformation = False
262SC.visualizationSettings.general.drawCoordinateSystem = False
263
264SC.visualizationSettings.loads.fixedLoadSize=0
265SC.visualizationSettings.loads.loadSizeFactor=0.5
266SC.visualizationSettings.loads.drawSimplified=False
267SC.visualizationSettings.loads.defaultSize=1
268SC.visualizationSettings.loads.defaultRadius=0.01
269
270SC.visualizationSettings.general.autoFitScene = True #otherwise, renderState not accepted for zoom
271
272#++++++++++++++++++++++++++++++++++++++++
273#setup simulation settings and run interactive dialog:
274simulationSettings = exu.SimulationSettings()
275simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1
276simulationSettings.solutionSettings.writeSolutionToFile = True
277simulationSettings.solutionSettings.solutionWritePeriod = 0.05 #data not used
278simulationSettings.solutionSettings.sensorsWritePeriod = 0.1 #data not used
279simulationSettings.solutionSettings.solutionInformation = 'n-mass-oscillatior'
280simulationSettings.timeIntegration.verboseMode = 1 #turn off, because of lots of output
281
282simulationSettings.timeIntegration.numberOfSteps = int(1000)
283simulationSettings.timeIntegration.endTime = 5
284simulationSettings.timeIntegration.newton.useModifiedNewton = True
285
286simulationSettings.displayComputationTime = True
287# simulationSettings.numberOfThreads = 1
288
289#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
290#set up interactive window
291
292
293mbs.SolveDynamic(simulationSettings=simulationSettings)
294
295mbs.SolutionViewer()