pendulumFriction.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Mathematical pendulum with friction;
  5#           Remark: uses two friction models: CoordinateSpringDamper and CartesianSpringDamper
  6#
  7# Author:   Johannes Gerstmayr
  8# Date:     2019-12-26
  9#
 10# 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.
 11#
 12#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 13
 14import exudyn as exu
 15from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
 16import exudyn.graphics as graphics #only import if it does not conflict
 17from exudyn.FEM import *
 18
 19import numpy as np
 20
 21useGraphics = True #without test
 22#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 23#you can erase the following lines and all exudynTestGlobals related operations if this is not intended to be used as TestModel:
 24try: #only if called from test suite
 25    from modelUnitTests import exudynTestGlobals #for globally storing test results
 26    useGraphics = exudynTestGlobals.useGraphics
 27except:
 28    class ExudynTestGlobals:
 29        pass
 30    exudynTestGlobals = ExudynTestGlobals()
 31#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 32
 33SC = exu.SystemContainer()
 34mbs = SC.AddSystem()
 35
 36
 37L = 0.8 #length of arm
 38mass = 2.5
 39g = 9.81
 40
 41r = 0.05 #just for graphics
 42d = r/2
 43
 44#add ground object and mass point:
 45graphicsBackground = GraphicsDataRectangle(-1.2*L,-1.2*L, 1.2*L, 0.2*L, [1,1,1,1]) #for appropriate zoom
 46oGround = mbs.AddObject(ObjectGround(referencePosition = [0,0,0],
 47                           visualization = VObjectGround(graphicsData = [graphicsBackground])))
 48
 49
 50graphicsSphere = graphics.Sphere(point=[L/2,0,0], radius=r, color=[1.,0.2,0.2,1], nTiles = 16)
 51graphicsSphere2 = graphics.Sphere(point=[0,0,0], radius=r, color=graphics.color.steelblue, nTiles = 16)
 52graphicsLink = graphics.BrickXYZ(-L/2,-d/2,-d/2, L/2,d/2, d/2, [0.5,0.5,0.5,0.5])
 53
 54inertia = InertiaCuboid(density=mass/(L*d*d), sideLengths=[L,d,d])
 55exu.Print("mass=",inertia.mass)
 56
 57nR0 = mbs.AddNode(Rigid2D(referenceCoordinates=[L/2,0,0])) #body goes from [0,0,0] to [L,0,0]
 58oR0 = mbs.AddObject(RigidBody2D(nodeNumber=nR0, physicsMass = inertia.mass, physicsInertia=inertia.inertiaTensor[2][2],
 59                                  visualization = VObjectRigidBody2D(graphicsData = [graphicsLink,graphicsSphere])))
 60
 61#markers:
 62mGround0 = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oGround, localPosition = [0,0,0]))
 63mR0 = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oR0, localPosition=[-L/2,0,0]))
 64mTip0 = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oR0, localPosition=[L/2,0,0]))
 65mNodeR0 = mbs.AddMarker(MarkerNodePosition(nodeNumber=nR0))
 66mR0COM = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oR0, localPosition=[0,0,0]))
 67
 68
 69oRJ0 = mbs.AddObject(RevoluteJoint2D(markerNumbers=[mGround0,mR0]))
 70#
 71mbs.AddLoad(Force(markerNumber = mNodeR0, loadVector = [0, -mass*g, 0]))
 72
 73zeroZoneFriction = 1e-3 #zero-zone for velocity in friction
 74fFriction = 1          #friction force (norm); acts against velocity
 75#user function for friction against velocity vector, including zeroZone
 76def UserFunctionSpringDamper(mbs, t, itemIndex, u, v, k, d, offset):
 77    vNorm = NormL2(v)
 78    f=[v[0],v[1],v[2]]
 79    if abs(vNorm) < offset[0]:
 80        f = ScalarMult(offset[1]/offset[0], f)
 81    else:
 82        f = ScalarMult(offset[1]/vNorm, f)
 83    return f
 84
 85mbs.AddObject(CartesianSpringDamper(markerNumbers=[mGround0, mTip0],
 86                                    offset=[zeroZoneFriction, fFriction, 0],
 87                                    springForceUserFunction=UserFunctionSpringDamper))
 88
 89if useGraphics:
 90    sRot1 = mbs.AddSensor(SensorBody(bodyNumber = oR0, fileName='solution/pendulumFrictionRotation0.txt',
 91                             outputVariableType=exu.OutputVariableType.Rotation))
 92
 93    sRot2 = mbs.AddSensor(SensorMarker(markerNumber = mR0COM, fileName='solution/pendulumFrictionRotation0marker.txt',
 94                               writeToFile = useGraphics,
 95                               outputVariableType=exu.OutputVariableType.Rotation))
 96
 97sPos = mbs.AddSensor(SensorMarker(markerNumber = mR0COM, writeToFile = False,
 98                           outputVariableType=exu.OutputVariableType.Position))
 99
100mbs.Assemble()
101
102simulationSettings = exu.SimulationSettings()
103
104f = 4000
105simulationSettings.timeIntegration.numberOfSteps = int(1*f)
106simulationSettings.timeIntegration.endTime = 0.0001*f
107simulationSettings.solutionSettings.solutionWritePeriod = simulationSettings.timeIntegration.endTime/5000
108simulationSettings.solutionSettings.sensorsWritePeriod = simulationSettings.timeIntegration.endTime/2000
109#simulationSettings.displayComputationTime = True
110simulationSettings.timeIntegration.verboseMode = 1
111
112#simulationSettings.timeIntegration.newton.useModifiedNewton = False
113simulationSettings.timeIntegration.generalizedAlpha.useNewmark = True
114simulationSettings.timeIntegration.generalizedAlpha.useIndex2Constraints = simulationSettings.timeIntegration.generalizedAlpha.useNewmark
115simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.60 #0.62 is approx. the limit
116
117simulationSettings.timeIntegration.generalizedAlpha.computeInitialAccelerations = True
118simulationSettings.solutionSettings.coordinatesSolutionFileName= "solution/coordinatesSolution.txt"
119simulationSettings.solutionSettings.writeSolutionToFile=False
120
121#simulationSettings.displayStatistics = True
122
123SC.visualizationSettings.nodes.defaultSize = 0.05
124
125if useGraphics:
126    exu.StartRenderer()
127    mbs.WaitForUserToContinue()
128
129mbs.SolveDynamic(simulationSettings)
130
131p0=mbs.GetObjectOutputBody(oR0, exu.OutputVariableType.Position, localPosition=[0,0,0])
132exu.Print("p0=", p0)
133
134p0 = mbs.GetSensorValues(sPos) #obtain values from marker
135exu.Print("p0=", p0, '(marker)')
136u=NormL2(p0)
137exu.Print('solution of pendulumFriction=',u)
138
139exudynTestGlobals.testError = u - (0.3999999877698205) #2020-04-22: 0.3999999877698205
140exudynTestGlobals.testResult = u
141
142
143if useGraphics:
144    SC.WaitForRenderEngineStopFlag()
145    exu.StopRenderer() #safely close rendering window!
146
147
148
149    mbs.PlotSensor([sRot1, sRot2], components=[0,2], closeAll=True, markerStyles=['x','+'])
150
151    # import matplotlib.pyplot as plt
152    # import matplotlib.ticker as ticker
153
154    # data = np.loadtxt('solution/pendulumFrictionRotation0.txt', comments='#', delimiter=',')
155    # plt.plot(data[:,0], data[:,1], 'b-', label='rotation 0') #ccordinate 1 = rotation, scalar for ObjectRigidBody2D
156    # data = np.loadtxt('solution/pendulumFrictionRotation0marker.txt', comments='#', delimiter=',')
157    # plt.plot(data[:,0], data[:,3], 'r-', label='rotation 0') #ccordinate 3 = rotation, Z-coordinate because marker always 3D
158
159    # ax=plt.gca() # get current axes
160    # ax.grid(True, 'major', 'both')
161    # ax.xaxis.set_major_locator(ticker.MaxNLocator(10))
162    # ax.yaxis.set_major_locator(ticker.MaxNLocator(10))
163    # plt.xlabel("time (s)")
164    # plt.ylabel("angle (rad)")
165    # plt.tight_layout() #better arrangement of plot
166    # plt.legend()
167    # plt.show()