cartesianSpringDamperUserFunction.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Example with CartesianSpringDamper, using symbolic user function for definition of spring-damper force
  5#
  6# Model:    Nonlinear oscillator with mass point and force user function
  7#
  8# Author:   Johannes Gerstmayr
  9# Date:     2023-12-07
 10#
 11# 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.
 12#
 13# *clean example*
 14#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 15
 16## import exudyn package and utilities, create system
 17import exudyn as exu
 18from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
 19import exudyn.graphics as graphics #only import if it does not conflict
 20
 21SC = exu.SystemContainer()
 22mbs = SC.AddSystem()
 23
 24## set abbreviation for symbolic library
 25esym = exu.symbolic
 26
 27
 28## define parameters for linear spring-damper
 29L=0.5
 30mass = 1.6
 31k = 4000
 32omega0 = 50 # sqrt(4000/1.6)
 33dRel = 0.05
 34d = dRel * 2 * 80 #80=sqrt(1.6*4000)
 35u0=-0.08
 36v0=1
 37f = 80
 38
 39## create ground object
 40objectGround = mbs.CreateGround(referencePosition = [0,0,0])
 41
 42## create mass point with initial displacement and initial velocity
 43massPoint = mbs.CreateMassPoint(referencePosition=[L,0,0],
 44                    initialDisplacement=[u0,0,0],
 45                    initialVelocity=[v0,0,0],
 46                    physicsMass=mass)
 47
 48## create spring damper between ground and mass point
 49csd = mbs.CreateCartesianSpringDamper(bodyNumbers=[objectGround, massPoint],
 50                                stiffness = [k,k,k],
 51                                damping   = [d,0,0],
 52                                offset    = [L,0,0])
 53
 54## add force on mass point
 55load = mbs.CreateForce(bodyNumber=massPoint, loadVector= [f,0,0])
 56
 57## add sensor to measure mass point position
 58sMass = mbs.AddSensor(SensorBody(bodyNumber=massPoint, storeInternal=True,
 59                         outputVariableType=exu.OutputVariableType.Position))
 60
 61## create user function for Cartesian spring damper (can be used as Python or as symbolic user function)
 62#force is just for demonstration and may not represent real behavior
 63def springForceUserFunction(mbs, t, itemNumber, u, v, k, d, offset):
 64    return [0.5*u[0]**2 * k[0]+esym.sign(v[0])*10,
 65            k[1]*u[1],
 66            k[2]*u[2]]
 67
 68CSDuserFunction = springForceUserFunction
 69doSymbolic = False
 70if doSymbolic:
 71    ## create symbolic user function (which can be used in the same way as the Python function)
 72    CSDuserFunction = CreateSymbolicUserFunction(mbs, springForceUserFunction,
 73                                                 'springForceUserFunction', csd)
 74    #check function:
 75    print('user function:\n',CSDuserFunction)
 76
 77mbs.SetObjectParameter(csd, 'springForceUserFunction', CSDuserFunction)
 78
 79## assemble and create simulation settings
 80mbs.Assemble()
 81
 82simulationSettings = exu.SimulationSettings()
 83
 84tEnd = 0.2*10
 85steps = 200000
 86simulationSettings.timeIntegration.numberOfSteps = steps
 87simulationSettings.timeIntegration.endTime = tEnd
 88simulationSettings.timeIntegration.verboseMode = 1
 89#simulationSettings.solutionSettings.solutionWritePeriod = tEnd/steps
 90simulationSettings.solutionSettings.writeSolutionToFile = False
 91simulationSettings.solutionSettings.sensorsWritePeriod = 0.001
 92
 93simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1 #SHOULD work with 0.9 as well
 94
 95## start renderer and solver; use explicit solver to account for switching in spring-damper
 96exu.StartRenderer()
 97mbs.SolveDynamic(simulationSettings,
 98                 solverType=exu.DynamicSolverType.ExplicitMidpoint)
 99# SC.WaitForRenderEngineStopFlag()
100exu.StopRenderer() #safely close rendering window!
101
102## evaluate solution
103n1 = mbs.GetObject(massPoint)['nodeNumber']
104u = mbs.GetNodeOutput(n1, exu.OutputVariableType.Position)
105
106print('u=',u)
107
108mbs.PlotSensor(sMass)