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