nMassOscillatorEigenmodes.py

#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# This is an EXUDYN example
#
# Details:  Nonlinear oscillations interactive simulation
#
# Author:   Johannes Gerstmayr
# Date:     2020-01-16
#
# 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.
#
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

import exudyn as exu
from exudyn.itemInterface import *
from exudyn.graphicsDataUtilities import *
import matplotlib.pyplot as plt
from exudyn.interactive import InteractiveDialog

import numpy as np
from math import sin, cos, pi, sqrt

import time #for sleep()
SC = exu.SystemContainer()
mbs = SC.AddSystem()

useConstraint = False
#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#
N = 12;                  #number of masses
spring = 800;           #stiffness [800]
mass = 1;               #mass
damper = 2;             #old:0.1; damping parameter
force = 1;              #force amplitude

d0 = damper*spring/(2*sqrt(mass*spring))  #dimensionless damping for single mass

stepSize = 0.002            #step size
endTime = 10 #time period to be simulated between every update

#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

omegaInit = 3.55
omegaMax = 40 #for plots

#ground node
nGround=mbs.AddNode(NodePointGround(referenceCoordinates = [0,0,0]))

#drawing parameters:
l_mass = 0.2          #spring length
r_mass = 0.030*2       #radius of mass
r_spring = r_mass*1.2
L0 = l_mass*1
L = N * l_mass + 4*l_mass
z=-r_mass-0.1
hy=0.25*L
hy1= hy
hy0=-hy
maxAmp0 = 0.1
maxAmpN = 0.1*N

background = [graphics.Quad([[-L0,hy0,z],[ L,hy0,z],[ L,hy1,z],[-L0,hy1,z]], 
                              color=graphics.color.lightgrey)]
    
oGround=mbs.AddObject(ObjectGround(visualization=VObjectGround(graphicsData=background)))

groundMarker=mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= nGround, coordinate = 0))
prevMarker = groundMarker
if useConstraint:
    prevMarker = None
    
nMass = []
markerList = []

for i in range(N+useConstraint):
    #node for 3D mass point:
    col = graphics.color.steelblue
    if i==int(useConstraint):
        col = graphics.color.green
    elif i==N-int(not useConstraint):
        col = graphics.color.lightred
    elif i==0 and useConstraint:
        col = graphics.color.lightgrey

    gSphere = graphics.Sphere(point=[0,0,0], radius=r_mass, color=col, nTiles=32)
    
    node = mbs.AddNode(Node1D(referenceCoordinates = [l_mass*(len(nMass)+1-useConstraint)],
                              initialCoordinates=[0.],
                              initialVelocities=[0.]
                              ))
    nMass += [node]
    massPoint = mbs.AddObject(Mass1D(nodeNumber = node, physicsMass=mass,
                                     referencePosition=[0,0,0],
                                     visualization=VMass1D(graphicsData=[gSphere])
                                     ))

    #marker for springDamper for first (x-)coordinate:
    nodeMarker =mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= node, coordinate = 0))
    markerList += [nodeMarker]

    #Spring-Damper between two marker coordinates
    if prevMarker!=None:
        sd = mbs.AddObject(CoordinateSpringDamper(markerNumbers = [prevMarker, nodeMarker], 
                                             stiffness = spring, damping = damper, 
                                             visualization=VCoordinateSpringDamper(drawSize=r_spring))) 

    prevMarker = nodeMarker

if useConstraint: #with constraints
    mbs.AddObject(CoordinateConstraint(markerNumbers=[groundMarker,markerList[0]],
                                       visualization=VCoordinateConstraint(show=False)))    

#add load to last mass:
if True: #scalar load
    mbs.AddLoad(LoadCoordinate(markerNumber = nodeMarker, 
                               load = 10)) #load set in user function


sensPos0 = mbs.AddSensor(SensorNode(nodeNumber=nMass[0], storeInternal=True,
                                    outputVariableType=exu.OutputVariableType.Coordinates))
sensPosN = mbs.AddSensor(SensorNode(nodeNumber=nMass[-1], storeInternal=True,
                                    outputVariableType=exu.OutputVariableType.Coordinates))

#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#compute eigenvalues
mbs.Assemble()
[values, vectors] = mbs.ComputeODE2Eigenvalues()
print('omegas (rad/s)=', np.sqrt(values))

#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#finalize model and settings
mbs.Assemble()


SC.visualizationSettings.general.textSize = 12
SC.visualizationSettings.openGL.lineWidth = 2
SC.visualizationSettings.openGL.multiSampling = 4
SC.visualizationSettings.general.graphicsUpdateInterval = 0.005
#SC.visualizationSettings.window.renderWindowSize=[1024,900]
SC.visualizationSettings.window.renderWindowSize=[1600,1000]
SC.visualizationSettings.general.showSolverInformation = False
SC.visualizationSettings.general.drawCoordinateSystem = False

SC.visualizationSettings.loads.fixedLoadSize=0
SC.visualizationSettings.loads.loadSizeFactor=0.5
SC.visualizationSettings.loads.drawSimplified=False
SC.visualizationSettings.loads.defaultSize=1
SC.visualizationSettings.loads.defaultRadius=0.01


#++++++++++++++++++++++++++++++++++++++++
#setup simulation settings and run interactive dialog:
simulationSettings = exu.SimulationSettings()
simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1
simulationSettings.solutionSettings.writeSolutionToFile = False
simulationSettings.solutionSettings.solutionWritePeriod = 0.1 #data not used
simulationSettings.solutionSettings.sensorsWritePeriod = 0.01 #data not used
simulationSettings.solutionSettings.solutionInformation = 'n-mass-oscillatior'
simulationSettings.timeIntegration.verboseMode = 0 #turn off, because of lots of output

simulationSettings.timeIntegration.numberOfSteps = int(endTime/stepSize)
simulationSettings.timeIntegration.endTime = endTime
simulationSettings.timeIntegration.newton.useModifiedNewton = True
#simulationSettings.timeIntegration.simulateInRealtime = True

simulationSettings.displayComputationTime = True

#plot FFT
if False:
    #%%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    #exu.StartRenderer()
    #mbs.WaitForUserToContinue()
    mbs.SolveDynamic(simulationSettings=simulationSettings)
    #exu.StopRenderer() #safely close rendering window!

    from exudyn.signalProcessing import ComputeFFT
    from exudyn.plot import PlotFFT
    data = mbs.GetSensorStoredData(sensPosN)
    [freq, amp, phase] = ComputeFFT(data[:,0], data[:,1])
    PlotFFT(freq, amp)

#%%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
if True:
    from exudyn.interactive import AnimateModes
    [values, systemEigenVectors] = mbs.ComputeODE2Eigenvalues()
    AnimateModes(SC, mbs, nodeNumber=None, systemEigenVectors=systemEigenVectors, 
                 runOnStart=True,)