springDamperTutorialNew.py

You can view and download this file on Github: springDamperTutorialNew.py

#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# This is an EXUDYN example
#
# Details:  This is the file for the EXUDYN first tutorial example showing a simple masspoint a SpringDamper
#
# Author:   Johannes Gerstmayr
# Date:     2023-05-15
#
# 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.utilities import *

import numpy as np #for postprocessing

SC = exu.SystemContainer()
mbs = SC.AddSystem()

print('EXUDYN version='+exu.GetVersionString())

L=0.5
mass = 1.6          #mass in kg
spring = 4000       #stiffness of spring-damper in N/m
damper = 8          #damping constant in N/(m/s)

u0=-0.08            #initial displacement
v0=1                #initial velocity
f =80               #force on mass
x0=f/spring         #static displacement

print('resonance frequency = '+str(np.sqrt(spring/mass)))
print('static displacement = '+str(x0))

oMass = mbs.CreateMassPoint(referencePosition=[L,0,0],
                            initialDisplacement = [u0,0,0],
                            initialVelocity= [v0,0,0],
                            physicsMass=mass) #force created via gravity

oGround = mbs.AddObject(ObjectGround())

#create spring damper with reference length computed from reference positions (=L)
oSD = mbs.CreateSpringDamper(bodyOrNodeList=[oMass, oGround],
                             stiffness = spring, damping = damper)

#add load via marker:
bodyMarker = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oMass))
mbs.AddLoad(LoadForceVector(markerNumber = bodyMarker, loadVector = [f,0,0]))


#add sensor:
sForce = mbs.AddSensor(SensorObject(objectNumber=oSD, storeInternal=True,
                                    outputVariableType=exu.OutputVariableType.ForceLocal))
sDisp = mbs.AddSensor(SensorBody(bodyNumber=oMass, storeInternal=True,
                                  outputVariableType=exu.OutputVariableType.Displacement))

print(mbs)
mbs.Assemble()

tEnd = 1     #end time of simulation
h = 0.001    #step size; leads to 1000 steps

simulationSettings = exu.SimulationSettings()
simulationSettings.solutionSettings.solutionWritePeriod = 5e-3  #output interval general
simulationSettings.solutionSettings.sensorsWritePeriod = 5e-3  #output interval of sensors
simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h) #must be integer
simulationSettings.timeIntegration.endTime = tEnd

#add some drawing parameters for this example
SC.visualizationSettings.nodes.drawNodesAsPoint=False
SC.visualizationSettings.nodes.defaultSize=0.1

exu.StartRenderer()              #start graphics visualization
mbs.WaitForUserToContinue()    #wait for pressing SPACE bar to continue

#start solver:
mbs.SolveDynamic(simulationSettings,
                 solverType=exu.DynamicSolverType.TrapezoidalIndex2)

SC.WaitForRenderEngineStopFlag()#wait for pressing 'Q' to quit
exu.StopRenderer()               #safely close rendering window!

#evaluate final (=current) output values
u = mbs.GetNodeOutput(0, exu.OutputVariableType.Position) #Node 0 is first node
print('displacement=',u)

#+++++++++++++++++++++++++++++++++++++++++++++++++++++
#compute exact solution:

omega0 = np.sqrt(spring/mass)     #eigen frequency of undamped system
dRel = damper/(2*np.sqrt(spring*mass)) #dimensionless damping
omega = omega0*np.sqrt(1-dRel**2) #eigen frequency of damped system
C1 = u0-x0 #static solution needs to be considered!
C2 = (v0+omega0*dRel*C1) / omega #C1, C2 are coeffs for solution
steps = int(tEnd/h)

refSol = np.zeros((steps+1,2))
for i in range(0,steps+1):
    t = tEnd*i/steps
    refSol[i,0] = t
    refSol[i,1] = np.exp(-omega0*dRel*t)*(C1*np.cos(omega*t) + C2*np.sin(omega*t))+x0

#use PlotSensor functionality to plot data:
mbs.PlotSensor(sensorNumbers=[refSol], components=[0], labels='displacement (m); exact solution',
               colorCodeOffset=2, closeAll=True) #color code offset to have same colors as in original example
mbs.PlotSensor(sensorNumbers=[sDisp], components=[0], labels='displacement (m); numerical solution',
               colorCodeOffset=0, newFigure=False)

mbs.PlotSensor(sensorNumbers=[sForce], labels='force (kN)',
               colorCodeOffset=1, factors=[1e-3], newFigure=False)