simulateInteractively.py

You can view and download this file on Github: simulateInteractively.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
 21
 22import time #for sleep()
 23SC = exu.SystemContainer()
 24mbs = SC.AddSystem()
 25
 26#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 27#create model for linear and nonlinear oscillator:
 28L=0.5
 29mass = 1.6          #mass in kg
 30spring = 4000       #stiffness of spring-damper in N/m
 31damper = 20         #damping constant in N/(m/s)
 32load0 = 80
 33
 34omega0=np.sqrt(spring/mass)
 35f0 = 0.*omega0/(2*np.pi)
 36f1 = 1.*omega0/(2*np.pi)
 37
 38print('resonance frequency = '+str(omega0/(2*pi)))
 39tEnd = 200     #end time of simulation
 40steps = 20000  #number of steps
 41
 42omegaInit = omega0*0.5
 43mbs.variables['mode'] = 0           #0=linear, 1=cubic nonlinear
 44mbs.variables['frequency'] = omegaInit/(2*pi)  #excitation frequency changed by user
 45#mbs.variables['omega'] = omegaInit  #excitation, changed in simFunction
 46mbs.variables['phi'] = 0            #excitation phase, used to get smooth excitations
 47mbs.variables['stiffness'] = spring
 48mbs.variables['damping'] = damper
 49
 50
 51
 52#user function for spring force
 53def springForce(mbs, t, itemIndex, u, v, k, d, offset): #changed 2023-01-21:, mu, muPropZone):
 54    k=mbs.variables['stiffness']
 55    d=mbs.variables['damping']
 56    if mbs.variables['mode'] == 0:
 57        return k*u + v*d
 58    else:
 59        #return 0.1*k*u+k*u**3+v*d
 60        return k*u+1000*k*u**3+v*d #breaks down at 13.40Hz
 61
 62mode = 0 #0...forward, 1...backward
 63
 64# #linear frequency sweep in time interval [0, t1] and frequency interval [f0,f1];
 65# def Sweep(t, t1, f0, f1):
 66#     k = (f1-f0)/t1
 67#     return np.sin(2*np.pi*(f0+k*0.5*t)*t) #take care of factor 0.5 in k*0.5*t, in order to obtain correct frequencies!!!
 68
 69# #user function for load
 70# def userLoad(mbs, t, load):
 71#     #return load*np.sin(0.5*omega0*t) #gives resonance
 72#     #print(t)
 73#     if mode==0:
 74#         return load*Sweep(t, tEnd, f0, f1)
 75#     else:
 76#         return load*Sweep(t, tEnd, f1, f0) #backward sweep
 77
 78#user function for load
 79def userLoad(mbs, t, load):
 80    #return load*sin(t*mbs.variables['frequency']*2*pi+mbs.variables['phi'])
 81    return load*sin(mbs.GetSensorValues(mbs.variables['sensorPhi']))
 82
 83#dummy user function for frequency
 84def userFrequency(mbs, t, load):
 85    return mbs.variables['frequency']
 86
 87#user function used in GenericODE2 to integrate current omega
 88def UFintegrateOmega(mbs, t, itemIndex, q, q_t):
 89    return [mbs.variables['frequency']*(2*pi)] #current frequency*2*pi is integrated into phi, return vector!
 90
 91#node for 3D mass point:
 92nMass=mbs.AddNode(Point(referenceCoordinates = [L,0,0]))
 93
 94#ground node
 95nGround=mbs.AddNode(NodePointGround(referenceCoordinates = [0,0,0]))
 96a=L
 97z=-0.1*L
 98background = graphics.Quad([[-0,-a,z],[ 2*a,-a,z],[ 2*a, a,z],[0, a,z]],
 99                              color=graphics.color.lightgrey, alternatingColor=graphics.color.white)
100oGround=mbs.AddObject(ObjectGround(visualization=VObjectGround(graphicsData=[background])))
101
102#add mass point (this is a 3D object with 3 coordinates):
103gCube = graphics.Brick([0.1*L,0,0], [0.2*L]*3, graphics.color.steelblue)
104massPoint = mbs.AddObject(MassPoint(physicsMass = mass, nodeNumber = nMass,
105                                    visualization=VMassPoint(graphicsData=[gCube])))
106
107#marker for ground (=fixed):
108groundMarker=mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= nGround, coordinate = 0))
109#marker for springDamper for first (x-)coordinate:
110nodeMarker  =mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= nMass, coordinate = 0))
111
112#Spring-Damper between two marker coordinates
113mbs.AddObject(CoordinateSpringDamper(markerNumbers = [groundMarker, nodeMarker],
114                                     stiffness = spring, damping = damper,
115                                     springForceUserFunction = springForce,
116                                     visualization=VCoordinateSpringDamper(drawSize=0.05)))
117
118#add load:
119mbs.AddLoad(LoadCoordinate(markerNumber = nodeMarker,
120                           load = load0, loadUserFunction=userLoad))
121
122#dummy load applied to ground marker, just to record/integrate frequency
123lFreq = mbs.AddLoad(LoadCoordinate(markerNumber = groundMarker,
124                           load = load0, loadUserFunction=userFrequency))
125
126sensPos = mbs.AddSensor(SensorNode(nodeNumber=nMass, fileName='solution/nonlinearPos.txt',
127                                   outputVariableType=exu.OutputVariableType.Displacement))
128sensVel = mbs.AddSensor(SensorNode(nodeNumber=nMass, fileName='solution/nonlinearVel.txt',
129                                   outputVariableType=exu.OutputVariableType.Velocity))
130sensFreq = mbs.AddSensor(SensorLoad(loadNumber=lFreq, fileName='solution/nonlinearFreq.txt',
131                                    visualization=VSensorLoad(show=False)))
132
133#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
134#node used to integrate omega into phi for excitation function
135nODE2=mbs.AddNode(NodeGenericODE2(referenceCoordinates=[0], initialCoordinates=[0],initialCoordinates_t=[0],
136                                  numberOfODE2Coordinates=1))
137
138oODE2=mbs.AddObject(ObjectGenericODE2(nodeNumbers=[nODE2],massMatrix=np.diag([1]),
139                                      forceUserFunction=UFintegrateOmega,
140                                      visualization=VObjectGenericODE2(show=False)))
141#improved version, using integration of omega:
142mbs.variables['sensorPhi'] = mbs.AddSensor(SensorNode(nodeNumber=nODE2, fileName='solution/nonlinearPhi.txt',
143                                    outputVariableType = exu.OutputVariableType.Coordinates_t,
144                                    visualization=VSensorNode(show=False)))
145#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
146
147
148mbs.Assemble()
149
150
151SC.visualizationSettings.general.textSize = 16
152SC.visualizationSettings.openGL.lineWidth = 2
153SC.visualizationSettings.openGL.multiSampling = 4
154SC.visualizationSettings.general.graphicsUpdateInterval = 0.02
155#SC.visualizationSettings.window.renderWindowSize=[1024,900]
156SC.visualizationSettings.window.renderWindowSize=[1200,1080]
157SC.visualizationSettings.general.showSolverInformation = False
158
159
160
161#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
162#this is an exemplary simulation function, which adjusts some values for simulation
163def SimulationUF(mbs, dialog):
164    #next two commands to zoom all ...:
165    if mbs.variables['mode'] == 1:
166        dialog.plots['limitsY'][0] = (-0.055,0.055)
167    else:
168        dialog.plots['limitsY'][0] = (-0.1,0.1)
169
170
171
172SC.visualizationSettings.general.autoFitScene = False #otherwise, renderState not accepted for zoom
173exu.StartRenderer()
174
175SC.SetRenderState({'centerPoint': [0.500249445438385, -0.02912527695298195, 0.0],
176 'maxSceneSize': 0.5,
177 'zoom': 0.428807526826858,
178 'currentWindowSize': [1400, 1200],
179 'modelRotation': [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]})
180time.sleep(0.5) #allow window to adjust view
181
182h = 1e-3      #step size of solver
183deltaT = 0.01 #time period to be simulated between every update
184
185#++++++++++++++++++++++++++++
186#define items for dialog
187dialogItems = [{'type':'label', 'text':'Nonlinear oscillation simulator', 'grid':(0,0,2), 'options':['L']},
188               {'type':'radio', 'textValueList':[('linear',0),('nonlinear (f=k*u+1000*k*u**3+d*v)',1)], 'value':0, 'variable':'mode', 'grid': [(2,0),(2,1)]},
189               {'type':'label', 'text':'excitation frequency (Hz):', 'grid':(5,0)},
190               {'type':'slider', 'range':(3*f1/800, 2.2*f1), 'value':omegaInit/(2*pi), 'steps':600, 'variable':'frequency', 'grid':(5,1)},
191               {'type':'label', 'text':'damping:', 'grid':(6,0)},
192               {'type':'slider', 'range': (0, 40), 'value':damper, 'steps':600, 'variable':'damping', 'grid':(6,1)},
193               {'type':'label', 'text':'stiffness:', 'grid':(7,0)},
194               {'type':'slider', 'range':(0, 10000), 'value':spring, 'steps':600, 'variable':'stiffness', 'grid':(7,1)}]
195
196#++++++++++++++++++++++++++++++++++++++++
197#specify subplots to be shown interactively
198plt.close('all')
199
200if False: #with phase
201    deltaT*=0.5 #higher resolution for phase
202    plots={'fontSize':16,'sizeInches':(12,12),'nPoints':200,
203           'subplots':(2,2), 'sensors':[[(sensPos,0),(sensPos,1),'time','mass position'],
204                                        [(sensFreq,0),(sensFreq,1),'time','excitation frequency'],
205                                        [(sensPos,1),(sensVel,1),'position (phase space)','velocity (phase space)']
206                                        ],
207           'limitsX':[(),(),()], #omit if time auto-range
208           'limitsY':[(-0.1,0.1),(0,2.2*f1*1.01),()]}
209else:
210    plots={'fontSize':16,'sizeInches':(12,12),'nPoints':400,
211           'subplots':(2,1), 'sensors':[[(sensPos,0),(sensPos,1),'time','mass position'],
212                                        [(sensFreq,0),(sensFreq,1),'time','excitation frequency']],
213           'limitsX':[(),()], #omit if time auto-range
214           'limitsY':[(-0.1,0.1),(0,2.2*f1*1.01)]}
215
216#++++++++++++++++++++++++++++++++++++++++
217#setup simulation settings and run interactive dialog:
218simulationSettings = exu.SimulationSettings()
219simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1
220simulationSettings.solutionSettings.writeSolutionToFile = False
221simulationSettings.solutionSettings.solutionWritePeriod = 0.1 #data not used
222simulationSettings.solutionSettings.sensorsWritePeriod = 0.1 #data not used
223simulationSettings.solutionSettings.solutionInformation = 'Nonlinear oscillations: compare linear / nonlinear case'
224simulationSettings.timeIntegration.verboseMode = 0 #turn off, because of lots of output
225
226simulationSettings.timeIntegration.numberOfSteps = int(deltaT/h)
227simulationSettings.timeIntegration.endTime = deltaT
228
229InteractiveDialog(mbs=mbs, simulationSettings=simulationSettings,
230                  simulationFunction=SimulationUF, title='Interactive window',
231                  dialogItems=dialogItems, period=deltaT, realtimeFactor=10,
232                  plots=plots, fontSize=12)
233
234# #stop solver and close render window
235exu.StopRenderer() #safely close rendering window!