SliderCrank.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Slider crank test model
  5#
  6# Author:   Johannes Gerstmayr
  7# Date:     2019-11-01
  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
 14#a local run must not open a plot window when windows are suppressed (#2477)
 15if exu.special.userInterface.suppressPlots: import matplotlib; matplotlib.use('Agg')
 16from exudyn.itemInterface import *
 17from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
 18import exudyn.graphics as graphics
 19
 20import numpy as np
 21import matplotlib.pyplot as plt
 22import matplotlib.ticker as ticker
 23
 24SC = exu.SystemContainer()
 25mbs = SC.AddSystem()
 26
 27#++++++++++++++++++++++++++++++++
 28#slider-crank
 29# test nonlinear model; index2 and index3-formulation for ConnectorCoordinate and RevoluteJoint2D
 30# crank is mounted at (0,0,0); crank length = 2*a0, connecting rod length = 2*a1
 31
 32#++++++++++++++++++++++++++++++++
 33#ground object/node:
 34
 35background = graphics.Lines([[-1,-2,0], [3,-2,0], [3,2,0], [-1,2,0], [-1,-2,0]], color=[0.9,0.9,0.9,1.])
 36
 37oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0], visualization=VObjectGround(graphicsData= [background])))
 38nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0])) #ground node for coordinate constraint
 39
 40#++++++++++++++++++++++++++++++++
 41#nodes and bodies
 42g = 9.81    # gravity
 43
 44a0 = 0.25     #half x-dim of body
 45b0 = 0.05    #half y-dim of body
 46massRigid0 = 2
 47inertiaRigid0 = massRigid0/12*(2*a0)**2
 48graphics0 = graphics.Lines([[-a0,-b0,0], [a0,-b0,0], [a0,b0,0], [-a0,b0,0], [-a0,-b0,0]])
 49
 50a1 = 0.5     #half x-dim of body
 51b1 = 0.05    #half y-dim of body
 52massRigid1 = 4
 53inertiaRigid1 = massRigid1/12*(2*a1)**2
 54graphics1 = graphics.Lines([[-a1,-b1,0], [a1,-b1,0], [a1,b1,0], [-a1,b1,0], [-a1,-b1,0]])
 55
 56nRigid0 = mbs.AddNode(Rigid2D(referenceCoordinates=[a0,0,0], 
 57                              initialVelocities=[0,0,0]));
 58oRigid0 = mbs.AddObject(RigidBody2D(mass=massRigid0, 
 59                                    inertia=inertiaRigid0,
 60                                    nodeNumber=nRigid0,
 61                                    visualization=VObjectRigidBody2D(graphicsData= [graphics0])))
 62
 63nRigid1 = mbs.AddNode(Rigid2D(referenceCoordinates=[2*a0+a1,0,0], initialVelocities=[0,0,0]));
 64oRigid1 = mbs.AddObject(RigidBody2D(mass=massRigid1, inertia=inertiaRigid1,nodeNumber=nRigid1,visualization=VObjectRigidBody2D(graphicsData= [graphics1])))
 65
 66c=0.05 #dimension of mass
 67sliderMass = 1
 68graphics2 = graphics.Lines([[-c,-c,0], [c,-c,0], [c,c,0], [-c,c,0], [-c,-c,0]])
 69
 70nMass = mbs.AddNode(Point2D(referenceCoordinates=[2*a0+2*a1,0]))
 71oMass = mbs.AddObject(MassPoint2D(mass=sliderMass, nodeNumber=nMass,visualization=VObjectRigidBody2D(graphicsData= [graphics2])))
 72
 73#++++++++++++++++++++++++++++++++
 74#markers for joints:
 75mR0Left = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oRigid0, localPosition=[-a0,0.,0.])) #support point # MUST be a rigidBodyMarker, because a torque is applied
 76mR0Right = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid0, localPosition=[ a0,0.,0.])) #end point; connection to connecting rod
 77
 78mR1Left = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid1, localPosition=[-a1,0.,0.])) #connection to crank
 79mR1Right = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid1, localPosition=[ a1,0.,0.])) #end point; connection to slider
 80
 81mMass = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oMass, localPosition=[ 0.,0.,0.]))
 82mG0 = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oGround, localPosition=[0,0,0.]))
 83
 84#++++++++++++++++++++++++++++++++
 85#joints:
 86mbs.AddObject(RevoluteJoint2D(markerNumbers=[mG0,mR0Left]))
 87mbs.AddObject(RevoluteJoint2D(markerNumbers=[mR0Right,mR1Left]))
 88mbs.AddObject(RevoluteJoint2D(markerNumbers=[mR1Right,mMass]))
 89
 90#++++++++++++++++++++++++++++++++
 91#markers for node constraints:
 92mGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nGround, coordinate=0)) #Ground node ==> no action
 93mNodeSlider = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nMass, coordinate=1)) #y-coordinate is constrained
 94
 95#++++++++++++++++++++++++++++++++
 96#coordinate constraints
 97mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mNodeSlider]))
 98
 99#loads and driving forces:
100mbs.AddLoad(Torque(markerNumber = mR0Left, loadVector = [0, 0, 10])) #apply torque at crank
101
102#++++++++++++++++++++++++++++++++
103#assemble, adjust settings and start time integration
104mbs.Assemble()
105
106#now as system is assembled, nodes know their global coordinate index (for reading the coordinate out of the solution file):
107#deprecated: globalIndex = mbs.CallNodeFunction(nMass, 'GetGlobalODE2CoordinateIndex')
108globalIndex = mbs.GetNodeODE2Index(nMass) 
109print('global ODE2 coordinate index of mass:', globalIndex)
110#alternatively: use mbs.systemData.GetObjectLTGODE2(oMass)[0] to obtain e.g. first coordinate index of sliding mass object
111
112simulationSettings = exu.SimulationSettings() #takes currently set values or default values
113
114simulationSettings.timeIntegration.numberOfSteps = 2*100000 #1000 steps for test suite/error
115simulationSettings.timeIntegration.endTime = 2              #1s for test suite / error
116#simulationSettings.timeIntegration.newton.relativeTolerance = 1e-10 #10000
117simulationSettings.timeIntegration.verboseMode = 1 #10000
118
119simulationSettings.solution.file.writePeriod = 1e-3
120
121simulationSettings.timeIntegration.newton.useModifiedNewton = True
122
123# simulationSettings.timeIntegration.generalizedAlpha.useNewmark = True
124# simulationSettings.timeIntegration.generalizedAlpha.useIndex2Constraints = True
125simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.5
126
127
128SC.renderer.Start()
129SC.renderer.DoIdleTasks()
130
131#++++++++++++++++++++++++++++++++++++++++++
132#solve generalized alpha / index3:
133mbs.SolveDynamic(simulationSettings)
134
135SC.renderer.DoIdleTasks()
136SC.renderer.Stop() #safely close rendering window!
137
138
139u = mbs.GetNodeOutput(nMass, exu.OutputVariableType.Position) #tip node
140errorSliderCrankIndex3 = u[0] - 1.3513750614331235 #x-position of slider
141print('error errorSliderCrankIndex3=',errorSliderCrankIndex3)
142
143dataIndex3 = np.loadtxt(OutputFilePath('solution/coordinatesSolution.txt'), comments='#', delimiter=',')
144
145#++++++++++++++++++++++++++++++++++++++++++
146##solve index 2 / trapezoidal rule:
147#simulationSettings.timeIntegration.generalizedAlpha.useNewmark = True
148#simulationSettings.timeIntegration.generalizedAlpha.useIndex2Constraints = True
149#
150#mbs.SolveDynamic(simulationSettings)
151#
152#u = mbs.GetNodeOutput(nMass, exu.OutputVariableType.Position) #tip node
153#errorSliderCrankIndex2 = u[0] - 1.3528786319585837 #x-position of slider
154#print('error errorSliderCrankIndex2=',errorSliderCrankIndex2)
155#
156#dataIndex2 = np.loadtxt(OutputFilePath('solution/coordinatesSolution.txt'), comments='#', delimiter=',')
157#plt.plot(dataIndex2[:,0], dataIndex2[:,1+globalIndex], 'r-') #plot x-coordinate of slider
158
159plt.plot(dataIndex3[:,0], dataIndex3[:,1+globalIndex], 'b-') #plot x-coordinate of slider
160
161ax=plt.gca() # get current axes
162ax.grid(True, 'major', 'both')
163ax.xaxis.set_major_locator(ticker.MaxNLocator(10)) #use maximum of 8 ticks on y-axis
164ax.yaxis.set_major_locator(ticker.MaxNLocator(10)) #use maximum of 8 ticks on y-axis
165plt.tight_layout()
166plt.show() 
167
168##animate solution
169#fileName = 'solution/coordinatesSolution.txt'
170#solution = LoadSolutionFile(OutputFilePath('solution/coordinatesSolution.txt'))
171#AnimateSolution(mbs, solution, 10, 0.05)