springDamperUserFunctionTest.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Test with user-defined load function and user-defined spring-damper function (Duffing oscillator)
  5#
  6# Author:   Johannes Gerstmayr
  7# Date:     2019-11-15
  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 *
 15
 16import numpy as np
 17
 18useGraphics = True #without test
 19#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 20#you can erase the following lines and all exudynTestGlobals related operations if this is not intended to be used as TestModel:
 21try: #only if called from test suite
 22    from modelUnitTests import exudynTestGlobals #for globally storing test results
 23    useGraphics = exudynTestGlobals.useGraphics
 24except:
 25    class ExudynTestGlobals:
 26        pass
 27    exudynTestGlobals = ExudynTestGlobals()
 28#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 29
 30SC = exu.SystemContainer()
 31mbs = SC.AddSystem()
 32exu.Print('EXUDYN version='+exu.GetVersionString())
 33
 34L=0.5
 35mass = 1.6          #mass in kg
 36spring = 4000       #stiffness of spring-damper in N/m
 37damper = 4          #damping constant in N/(m/s)
 38load0 = 80
 39
 40omega0=np.sqrt(spring/mass)
 41f0 = 0.*omega0/(2*np.pi)
 42f1 = 1.*omega0/(2*np.pi)
 43
 44exu.Print('resonance frequency = '+str(omega0))
 45tEnd = 50     #end time of simulation
 46steps = 5000  #number of steps
 47
 48#tEnd = 0.05     #end time of simulation
 49#steps = 5  #number of steps
 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    return 0.1*k*u+k*u**3+v*d
 55
 56#linear frequency sweep in time interval [0, t1] and frequency interval [f0,f1];
 57def Sweep(t, t1, f0, f1):
 58    k = (f1-f0)/t1
 59    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!!!
 60
 61#user function for load
 62def userLoad(mbs, t, load):
 63    #return load*np.sin(0.5*omega0*t) #gives resonance
 64    #exu.Print(t)
 65    return load*Sweep(t, tEnd, f0, f1)
 66    #return load*Sweep(t, tEnd, f1, f0) #backward sweep
 67
 68#node for 3D mass point:
 69n1=mbs.AddNode(Point(referenceCoordinates = [L,0,0]))
 70
 71#ground node
 72nGround=mbs.AddNode(NodePointGround(referenceCoordinates = [0,0,0]))
 73
 74#add mass point (this is a 3D object with 3 coordinates):
 75massPoint = mbs.AddObject(MassPoint(physicsMass = mass, nodeNumber = n1))
 76
 77#marker for ground (=fixed):
 78groundMarker=mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= nGround, coordinate = 0))
 79#marker for springDamper for first (x-)coordinate:
 80nodeMarker  =mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= n1, coordinate = 0))
 81
 82#Spring-Damper between two marker coordinates
 83mbs.AddObject(CoordinateSpringDamper(markerNumbers = [groundMarker, nodeMarker],
 84                                     stiffness = spring, damping = damper,
 85                                     springForceUserFunction = springForce))
 86
 87#add load:
 88loadC = mbs.AddLoad(LoadCoordinate(markerNumber = nodeMarker,
 89                           load = load0, loadUserFunction=userLoad))
 90
 91writeSensorFile = False
 92if useGraphics:
 93    writeSensorFile = True
 94
 95sLoad=mbs.AddSensor(SensorLoad(loadNumber=loadC, writeToFile = writeSensorFile,
 96                         storeInternal=True,#fileName="solution/userFunctionLoad.txt"
 97                         ))
 98#mbs.AddSensor(SensorNode(nodeNumber=n1, writeToFile = writeSensorFile, fileName="solution/userFunctionNode.txt"))
 99sCoords=mbs.AddSensor(SensorNode(nodeNumber=n1, writeToFile = writeSensorFile,
100                         outputVariableType=exu.OutputVariableType.Coordinates,
101                         storeInternal=True,#fileName="solution/userFunctionNode.txt"
102                         ))
103
104#exu.Print(mbs)
105mbs.Assemble()
106
107simulationSettings = exu.SimulationSettings()
108simulationSettings.solutionSettings.solutionWritePeriod = 2e-3  #output interval
109simulationSettings.timeIntegration.numberOfSteps = steps
110simulationSettings.timeIntegration.endTime = tEnd
111
112simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1
113
114simulationSettings.displayStatistics = True
115simulationSettings.timeIntegration.verboseMode = 1
116
117#exu.StartRenderer()              #start graphics visualization
118#mbs.WaitForUserToContinue()    #wait for pressing SPACE bar to continue
119
120#start solver:
121mbs.SolveDynamic(simulationSettings)
122
123#SC.WaitForRenderEngineStopFlag()#wait for pressing 'Q' to quit
124#exu.StopRenderer()               #safely close rendering window!
125
126#evaluate final (=current) output values
127u = mbs.GetNodeOutput(n1, exu.OutputVariableType.Position)
128exu.Print('displacement=',u[0])
129
130exudynTestGlobals.testError = u[0] - (0.5062872273010898) #2019-12-18: 0.5062872273010898; #2019-12-15: 0.5062872272996835; 2019-12-13:0.5062872273014417; 2019-12-01: 0.5152217339585201
131exudynTestGlobals.testResult = u[0]
132
133#+++++++++++++++++++++++++++++++++++++++++++++++++++++
134
135if useGraphics:
136
137
138    mbs.PlotSensor(sCoords, components=[0], closeAll=True)