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