nMassOscillatorEigenmodes.py

You can view and download this file on Github: nMassOscillatorEigenmodes.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
 26useConstraint = False
 27#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 28#
 29N = 12;                  #number of masses
 30spring = 800;           #stiffness [800]
 31mass = 1;               #mass
 32damper = 2;             #old:0.1; damping parameter
 33force = 1;              #force amplitude
 34
 35d0 = damper*spring/(2*sqrt(mass*spring))  #dimensionless damping for single mass
 36
 37stepSize = 0.002            #step size
 38endTime = 10 #time period to be simulated between every update
 39
 40#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 41
 42omegaInit = 3.55
 43omegaMax = 40 #for plots
 44
 45#ground node
 46nGround=mbs.AddNode(NodePointGround(referenceCoordinates = [0,0,0]))
 47
 48#drawing parameters:
 49l_mass = 0.2          #spring length
 50r_mass = 0.030*2       #radius of mass
 51r_spring = r_mass*1.2
 52L0 = l_mass*1
 53L = N * l_mass + 4*l_mass
 54z=-r_mass-0.1
 55hy=0.25*L
 56hy1= hy
 57hy0=-hy
 58maxAmp0 = 0.1
 59maxAmpN = 0.1*N
 60
 61background = [graphics.Quad([[-L0,hy0,z],[ L,hy0,z],[ L,hy1,z],[-L0,hy1,z]],
 62                              color=graphics.color.lightgrey)]
 63
 64oGround=mbs.AddObject(ObjectGround(visualization=VObjectGround(graphicsData=background)))
 65
 66groundMarker=mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= nGround, coordinate = 0))
 67prevMarker = groundMarker
 68if useConstraint:
 69    prevMarker = None
 70
 71nMass = []
 72markerList = []
 73
 74for i in range(N+useConstraint):
 75    #node for 3D mass point:
 76    col = graphics.color.steelblue
 77    if i==int(useConstraint):
 78        col = graphics.color.green
 79    elif i==N-int(not useConstraint):
 80        col = graphics.color.lightred
 81    elif i==0 and useConstraint:
 82        col = graphics.color.lightgrey
 83
 84    gSphere = graphics.Sphere(point=[0,0,0], radius=r_mass, color=col, nTiles=32)
 85
 86    node = mbs.AddNode(Node1D(referenceCoordinates = [l_mass*(len(nMass)+1-useConstraint)],
 87                              initialCoordinates=[0.],
 88                              initialVelocities=[0.]
 89                              ))
 90    nMass += [node]
 91    massPoint = mbs.AddObject(Mass1D(nodeNumber = node, physicsMass=mass,
 92                                     referencePosition=[0,0,0],
 93                                     visualization=VMass1D(graphicsData=[gSphere])
 94                                     ))
 95
 96    #marker for springDamper for first (x-)coordinate:
 97    nodeMarker =mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= node, coordinate = 0))
 98    markerList += [nodeMarker]
 99
100    #Spring-Damper between two marker coordinates
101    if prevMarker!=None:
102        sd = mbs.AddObject(CoordinateSpringDamper(markerNumbers = [prevMarker, nodeMarker],
103                                             stiffness = spring, damping = damper,
104                                             visualization=VCoordinateSpringDamper(drawSize=r_spring)))
105
106    prevMarker = nodeMarker
107
108if useConstraint: #with constraints
109    mbs.AddObject(CoordinateConstraint(markerNumbers=[groundMarker,markerList[0]],
110                                       visualization=VCoordinateConstraint(show=False)))
111
112#add load to last mass:
113if True: #scalar load
114    mbs.AddLoad(LoadCoordinate(markerNumber = nodeMarker,
115                               load = 10)) #load set in user function
116
117
118sensPos0 = mbs.AddSensor(SensorNode(nodeNumber=nMass[0], storeInternal=True,
119                                    outputVariableType=exu.OutputVariableType.Coordinates))
120sensPosN = mbs.AddSensor(SensorNode(nodeNumber=nMass[-1], storeInternal=True,
121                                    outputVariableType=exu.OutputVariableType.Coordinates))
122
123#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
124#compute eigenvalues
125mbs.Assemble()
126[values, vectors] = mbs.ComputeODE2Eigenvalues()
127print('omegas (rad/s)=', np.sqrt(values))
128
129#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
130#finalize model and settings
131mbs.Assemble()
132
133
134SC.visualizationSettings.general.textSize = 12
135SC.visualizationSettings.openGL.lineWidth = 2
136SC.visualizationSettings.openGL.multiSampling = 4
137SC.visualizationSettings.general.graphicsUpdateInterval = 0.005
138#SC.visualizationSettings.window.renderWindowSize=[1024,900]
139SC.visualizationSettings.window.renderWindowSize=[1600,1000]
140SC.visualizationSettings.general.showSolverInformation = False
141SC.visualizationSettings.general.drawCoordinateSystem = False
142
143SC.visualizationSettings.loads.fixedLoadSize=0
144SC.visualizationSettings.loads.loadSizeFactor=0.5
145SC.visualizationSettings.loads.drawSimplified=False
146SC.visualizationSettings.loads.defaultSize=1
147SC.visualizationSettings.loads.defaultRadius=0.01
148
149
150#++++++++++++++++++++++++++++++++++++++++
151#setup simulation settings and run interactive dialog:
152simulationSettings = exu.SimulationSettings()
153simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1
154simulationSettings.solutionSettings.writeSolutionToFile = False
155simulationSettings.solutionSettings.solutionWritePeriod = 0.1 #data not used
156simulationSettings.solutionSettings.sensorsWritePeriod = 0.01 #data not used
157simulationSettings.solutionSettings.solutionInformation = 'n-mass-oscillatior'
158simulationSettings.timeIntegration.verboseMode = 0 #turn off, because of lots of output
159
160simulationSettings.timeIntegration.numberOfSteps = int(endTime/stepSize)
161simulationSettings.timeIntegration.endTime = endTime
162simulationSettings.timeIntegration.newton.useModifiedNewton = True
163#simulationSettings.timeIntegration.simulateInRealtime = True
164
165simulationSettings.displayComputationTime = True
166
167#plot FFT
168if False:
169    #%%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
170    #exu.StartRenderer()
171    #mbs.WaitForUserToContinue()
172    mbs.SolveDynamic(simulationSettings=simulationSettings)
173    #exu.StopRenderer() #safely close rendering window!
174
175    from exudyn.signalProcessing import ComputeFFT
176    from exudyn.plot import PlotFFT
177    data = mbs.GetSensorStoredData(sensPosN)
178    [freq, amp, phase] = ComputeFFT(data[:,0], data[:,1])
179    PlotFFT(freq, amp)
180
181#%%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
182if True:
183    from exudyn.interactive import AnimateModes
184    [values, systemEigenVectors] = mbs.ComputeODE2Eigenvalues()
185    AnimateModes(SC, mbs, nodeNumber=None, systemEigenVectors=systemEigenVectors,
186                 runOnStart=True,)