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