SliderCrank2DTest.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  A planar slider-crank of two rigid bodies and a sliding mass, integrated over 1 s **as an
  5#           index 3 and as an index 2 system**: the x-position of the slider from both, added up, so
  6#           that the test covers the two constraint formulations at once.
  7#           The model compares against a reference value written into it, so its
  8#           result is that difference and its reference solution is 0 (#2632).
  9#
 10# Author:   Johannes Gerstmayr
 11# Date:     2019-11-01, reworked 2026-09-24
 12#
 13# 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.
 14#
 15#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 16
 17import exudyn as exu
 18from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
 19import exudyn.graphics as graphics
 20
 21testIsActive = exu.sys.get('testIsActive', False)
 22exu.sys['testTolerance'] = 4e-13 #the tolerance RunAllModelUnitTests used for these ten
 23
 24SC = exu.SystemContainer()
 25mbs = SC.AddSystem()
 26
 27
 28#++++++++++++++++++++++++++++++++
 29#ground object/node:
 30
 31rect = [-1,-2,3,2] #xmin,ymin,xmax,ymax
 32background = 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.])
 33#{'type':'Line', 'color':[0.1,0.1,0.8,1], 'data':[rect[0],rect[1],0, rect[2],rect[1],0, rect[2],rect[3],0, rect[0],rect[3],0, rect[0],rect[1],0]} #background
 34oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0], visualization=VObjectGround(graphicsData= [background])))
 35nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0])) #ground node for coordinate constraint
 36
 37#++++++++++++++++++++++++++++++++
 38#nodes and bodies
 39#crank is mounted at (0,0,0); crank length = 2*a0, connecting rod length = 2*a1
 40a0 = 0.25     #half x-dim of body
 41b0 = 0.05    #half y-dim of body
 42massRigid0 = 2
 43inertiaRigid0 = massRigid0/12*(2*a0)**2
 44graphics0 = graphics.Lines([[-a0,-b0,0], [a0,-b0,0], [a0,b0,0], [-a0,b0,0], [-a0,-b0,0]])
 45#{'type':'Line', 'color':[0.1,0.1,0.8,1], 'data':[-a0,-b0,0, a0,-b0,0, a0,b0,0, -a0,b0,0, -a0,-b0,0]} #background
 46
 47a1 = 0.5     #half x-dim of body
 48b1 = 0.05    #half y-dim of body
 49massRigid1 = 4
 50inertiaRigid1 = massRigid1/12*(2*a1)**2
 51graphics1 = graphics.Lines([[-a1,-b1,0], [a1,-b1,0], [a1,b1,0], [-a1,b1,0], [-a1,-b1,0]])
 52
 53nRigid0 = mbs.AddNode(Rigid2D(referenceCoordinates=[a0,0,0], initialVelocities=[0,0,0]));
 54oRigid0 = mbs.AddObject(RigidBody2D(mass=massRigid0, inertia=inertiaRigid0,nodeNumber=nRigid0,visualization=VObjectRigidBody2D(graphicsData= [graphics0])))
 55
 56nRigid1 = mbs.AddNode(Rigid2D(referenceCoordinates=[2*a0+a1,0,0], initialVelocities=[0,0,0]));
 57oRigid1 = mbs.AddObject(RigidBody2D(mass=massRigid1, inertia=inertiaRigid1,nodeNumber=nRigid1,visualization=VObjectRigidBody2D(graphicsData= [graphics1])))
 58
 59c=0.05 #dimension of mass
 60sliderMass = 1
 61graphics2 = graphics.Lines([[-c,-c,0], [c,-c,0], [c,c,0], [-c,c,0], [-c,-c,0]])
 62
 63nMass = mbs.AddNode(Point2D(referenceCoordinates=[2*a0+2*a1,0]))
 64oMass = mbs.AddObject(MassPoint2D(mass=sliderMass, nodeNumber=nMass,visualization=VObjectRigidBody2D(graphicsData= [graphics2])))
 65
 66#++++++++++++++++++++++++++++++++
 67#markers for joints:
 68mR0Left = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oRigid0, localPosition=[-a0,0.,0.])) #support point # MUST be a rigidBodyMarker, because a torque is applied
 69mR0Right = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid0, localPosition=[ a0,0.,0.])) #end point; connection to connecting rod
 70
 71mR1Left = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid1, localPosition=[-a1,0.,0.])) #connection to crank
 72mR1Right = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid1, localPosition=[ a1,0.,0.])) #end point; connection to slider
 73
 74mMass = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oMass, localPosition=[ 0.,0.,0.]))
 75mG0 = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oGround, localPosition=[0,0,0.]))
 76
 77#++++++++++++++++++++++++++++++++
 78#joints:
 79mbs.AddObject(RevoluteJoint2D(markerNumbers=[mG0,mR0Left]))
 80mbs.AddObject(RevoluteJoint2D(markerNumbers=[mR0Right,mR1Left]))
 81mbs.AddObject(RevoluteJoint2D(markerNumbers=[mR1Right,mMass]))
 82
 83#++++++++++++++++++++++++++++++++
 84#markers for node constraints:
 85mGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nGround, coordinate=0)) #Ground node ==> no action
 86mNodeSlider = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nMass, coordinate=1)) #y-coordinate is constrained
 87
 88#++++++++++++++++++++++++++++++++
 89#coordinate constraints
 90mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mNodeSlider]))
 91
 92#loads and driving forces:
 93mbs.AddLoad(Torque(markerNumber = mR0Left, loadVector = [0, 0, 10])) #apply torque at crank
 94
 95#++++++++++++++++++++++++++++++++
 96#assemble, adjust settings and start time integration
 97mbs.Assemble()
 98
 99simulationSettings = exu.SimulationSettings() #takes currently set values or default values
100simulationSettings.solution.file.write=False
101
102simulationSettings.timeIntegration.numberOfSteps = 1000
103simulationSettings.timeIntegration.endTime = 1
104simulationSettings.timeIntegration.newton.useModifiedNewton = True
105
106simulationSettings.timeIntegration.newton.relativeTolerance = 1e-10 #10000
107simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.5
108simulationSettings.show.statistics = False
109
110if not testIsActive: 
111    SC.renderer.Start()
112
113#solve generalized alpha / index3:
114mbs.SolveDynamic(simulationSettings)
115
116u = mbs.GetNodeOutput(nMass, exu.OutputVariableType.Position) #tip node
117if True:
118    errorSliderCrankIndex3 = u[0] - 1.3550008762955048 #2021-02-04: 1.3550008762955048
119else:
120    errorSliderCrankIndex3 = u[0] - 1.353298442702153  #2019-12-26: 1.353298442702153; 15.12.2019: 1.3513750614337234; before 15.12.2019: 1.3513750614326427 #2019-11-22; previous: 1.3513750614331235 #x-position of slider
121exu.Print('solution SliderCrankIndex3  =',u[0])
122exu.Print('error errorSliderCrankIndex3=',errorSliderCrankIndex3)
123
124simulationSettings.timeIntegration.generalizedAlpha.useNewmark = True
125simulationSettings.timeIntegration.generalizedAlpha.useIndex2Constraints = True
126
127#solve index 2 / trapezoidal rule:
128mbs.SolveDynamic(simulationSettings)
129
130u = mbs.GetNodeOutput(nMass, exu.OutputVariableType.Position) #tip node
131errorSliderCrankIndex2 = u[0] - 1.3550413308333111 #2019-12-26: 1.3550413308333111; 15.12.2019: 1.352878631961969; before 15.12.2019: 1.3528786319585846 #2019-11-22; previous: 1.3528786319585837 #x-position of slider
132exu.Print('solution SliderCrankIndex2  =',u[0])
133exu.Print('error errorSliderCrankIndex2=',errorSliderCrankIndex2)
134
135if not testIsActive: 
136    SC.renderer.DoIdleTasks()
137    SC.renderer.Stop() #safely close rendering window!
138
139
140testResult = abs(errorSliderCrankIndex3)+abs(errorSliderCrankIndex2)
141exu.Print('solution of SliderCrank2DTest=', testResult)
142exu.sys['testResult'] = testResult