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
 18testIsActive = exu.sys.get('testIsActive', False)
 19
 20SC = exu.SystemContainer()
 21mbs = SC.AddSystem()
 22exu.Print('EXUDYN version='+exu.config.Version())
 23
 24L=0.5
 25mass = 1.6          #mass in kg
 26spring = 4000       #stiffness of spring-damper in N/m
 27damper = 4          #damping constant in N/(m/s)
 28load0 = 80
 29
 30omega0=np.sqrt(spring/mass)
 31f0 = 0.*omega0/(2*np.pi)
 32f1 = 1.*omega0/(2*np.pi)
 33
 34exu.Print('resonance frequency = '+str(omega0))
 35tEnd = 50     #end time of simulation
 36steps = 5000  #number of steps
 37
 38#tEnd = 0.05     #end time of simulation
 39#steps = 5  #number of steps
 40
 41
 42#user function for spring force
 43def springForce(mbs, t, itemIndex, u, v, k, d, offset): #changed 2023-01-21:, mu, muPropZone):
 44    return 0.1*k*u+k*u**3+v*d
 45
 46#linear frequency sweep in time interval [0, t1] and frequency interval [f0,f1];
 47def Sweep(t, t1, f0, f1):
 48    k = (f1-f0)/t1
 49    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!!!
 50
 51#user function for load
 52def userLoad(mbs, t, load):
 53    #return load*np.sin(0.5*omega0*t) #gives resonance
 54    #exu.Print(t)
 55    return load*Sweep(t, tEnd, f0, f1)
 56    #return load*Sweep(t, tEnd, f1, f0) #backward sweep
 57
 58#node for 3D mass point:
 59n1=mbs.AddNode(Point(referenceCoordinates = [L,0,0]))
 60
 61#ground node
 62nGround=mbs.AddNode(NodePointGround(referenceCoordinates = [0,0,0]))
 63
 64#add mass point (this is a 3D object with 3 coordinates):
 65massPoint = mbs.AddObject(MassPoint(mass = mass, nodeNumber = n1))
 66
 67#marker for ground (=fixed):
 68groundMarker=mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= nGround, coordinate = 0))
 69#marker for springDamper for first (x-)coordinate:
 70nodeMarker  =mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= n1, coordinate = 0))
 71
 72#Spring-Damper between two marker coordinates
 73mbs.AddObject(CoordinateSpringDamper(markerNumbers = [groundMarker, nodeMarker], 
 74                                     stiffness = spring, damping = damper, 
 75                                     springForceUserFunction = springForce)) 
 76
 77#add load:
 78loadC = mbs.AddLoad(LoadCoordinate(markerNumber = nodeMarker, 
 79                           load = load0, loadUserFunction=userLoad))
 80
 81writeSensorFile = False
 82if not testIsActive:
 83    writeSensorFile = True
 84
 85sLoad=mbs.AddSensor(SensorLoad(loadNumber=loadC, writeToFile = writeSensorFile, 
 86                         storeInternal=True,#fileName="solution/userFunctionLoad.txt"
 87                         ))
 88#mbs.AddSensor(SensorNode(nodeNumber=n1, writeToFile = writeSensorFile, fileName="solution/userFunctionNode.txt"))
 89sCoords=mbs.AddSensor(SensorNode(nodeNumber=n1, writeToFile = writeSensorFile, 
 90                         outputVariableType=exu.OutputVariableType.CoordinatesTotal, 
 91                         storeInternal=True,#fileName="solution/userFunctionNode.txt"
 92                         ))
 93    
 94#exu.Print(mbs)
 95mbs.Assemble()
 96
 97simulationSettings = exu.SimulationSettings()
 98simulationSettings.timeIntegration.newton.useModifiedNewton = False #Just for the test; modified Newton is usually faster
 99simulationSettings.solution.file.writePeriod = 2e-3  #output interval
100simulationSettings.timeIntegration.numberOfSteps = steps
101simulationSettings.timeIntegration.endTime = tEnd
102
103simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1
104
105simulationSettings.show.statistics = True
106simulationSettings.timeIntegration.verboseMode = 1
107
108#SC.renderer.Start()              #start graphics visualization
109#SC.renderer.DoIdleTasks()    #wait for pressing SPACE bar to continue
110
111#start solver:q
112mbs.SolveDynamic(simulationSettings)
113
114#SC.renderer.DoIdleTasks()#wait for pressing 'Q' to quit
115#SC.renderer.Stop()               #safely close rendering window!
116
117#evaluate final (=current) output values
118# u = mbs.GetNodeOutput(n1, exu.OutputVariableType.Position)
119# exu.Print('u     =',u)
120uTotal = mbs.GetNodeOutput(n1, exu.OutputVariableType.CoordinatesTotal)
121exu.Print('uTotal=',uTotal[0])
122
123exu.sys['testResult'] = uTotal[0]
124
125#+++++++++++++++++++++++++++++++++++++++++++++++++++++
126
127if not testIsActive:
128    
129    
130    mbs.PlotSensor(sCoords, components=[0], closeAll=True)