fourBarMechanismTest.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Test model for four bar mechanism
  5#
  6# Author:   Johannes Gerstmayr
  7# Date:     2019-08-15
  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 *
 15
 16useGraphics = True #without test
 17#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 18#you can erase the following lines and all exudynTestGlobals related operations if this is not intended to be used as TestModel:
 19try: #only if called from test suite
 20    from modelUnitTests import exudynTestGlobals #for globally storing test results
 21    useGraphics = exudynTestGlobals.useGraphics
 22except:
 23    class ExudynTestGlobals:
 24        pass
 25    exudynTestGlobals = ExudynTestGlobals()
 26#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 27
 28SC = exu.SystemContainer()
 29mbs = SC.AddSystem()
 30
 31nBodies = 4
 32mass = 3 #kg
 33g = 9.81 # m/s^2
 34omega0 = 2 #initial angular velocity of first link
 35
 36#add nodes:
 37n0=mbs.AddNode(PointGround(referenceCoordinates=[0, 0, 0]))
 38n1=mbs.AddNode(Point(referenceCoordinates=[0, 2, 0], initialVelocities= [2*omega0,0,0]))
 39n2=mbs.AddNode(Point(referenceCoordinates = [4, 5, 0], initialVelocities= [5*omega0,0,0]))
 40n3=mbs.AddNode(Point(referenceCoordinates = [4, 0, 0]))
 41
 42#nodetypes Point2DSlope1Slope2 Point2DSlope1 RigidBody RigidBody2D
 43
 44L=4
 45Lo=0.
 46#add mass points and ground object:
 47rect = [-6,-6,6,6] #xmin,ymin,xmax,ymax
 48background = {'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
 49#oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0], visualization=VObjectGround(graphicsData= [background])))
 50#graphics1 = {'type':'Line', 'color':[0.1,0.1,0.8,1], 'data':[-(L+Lo),-(L+Lo),0, (L+Lo),-(L+Lo),0, (L+Lo),Lo,0, -(L+Lo),Lo,0, -(L+Lo),-(L+Lo), 0]} #background
 51o0=mbs.AddObject(MassPoint(physicsMass=0,nodeNumber=n0,visualization=VObjectMassPoint(graphicsData= [background])))
 52o1=mbs.AddObject(MassPoint(physicsMass=mass,nodeNumber=n1))
 53o2=mbs.AddObject(MassPoint(physicsMass=mass,nodeNumber=n2))
 54o3=mbs.AddObject(MassPoint(physicsMass=mass,nodeNumber=n3))
 55
 56#use NodePosition markers:
 57for i in range(nBodies):
 58    mbs.AddMarker(MarkerNodePosition(nodeNumber = i))
 59
 60#Coordinate Marker:
 61mc0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = n0, coordinate=0)) #Ground node ==> no action
 62mc1 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = n3, coordinate=1))
 63mc2 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = n3, coordinate=0))
 64mbs.AddObject(CoordinateConstraint(markerNumbers=[mc0,mc1]))
 65mbs.AddObject(CoordinateConstraint(markerNumbers=[mc0,mc2]))
 66
 67useConstraint = True
 68if useConstraint:
 69    #add distance:
 70    mbs.AddObject(ObjectConnectorDistance(markerNumbers=[0,1], distance=2, visualization=VObjectConnectorDistance(drawSize=0.01)))
 71    mbs.AddObject(ObjectConnectorDistance(markerNumbers=[1,2], distance=5, visualization=VObjectConnectorDistance(drawSize=0.01)))
 72    mbs.AddObject(ObjectConnectorDistance(markerNumbers=[2,3], distance=5, visualization=VObjectConnectorDistance(drawSize=0.01)))
 73else:
 74    #add spring-dampers:
 75    k = 40000
 76    d = 20
 77    mbs.AddObject({'objectType': 'ConnectorSpringDamper', 'stiffness': k, 'damping': d, 'force': 0, 'referenceLength':2, 'markerNumbers': [0,1], 'VdrawSize': 0.01})
 78    mbs.AddObject({'objectType': 'ConnectorSpringDamper', 'stiffness': k, 'damping': d, 'force': 0, 'referenceLength':5, 'markerNumbers': [1,2], 'VdrawSize': 0.01})
 79    mbs.AddObject({'objectType': 'ConnectorSpringDamper', 'stiffness': k, 'damping': d, 'force': 0, 'referenceLength':5, 'markerNumbers': [2,3], 'VdrawSize': 0.01})
 80
 81#add loads:
 82mbs.AddLoad(Force(markerNumber = 1, loadVector = [0, -mass*g, 0]))
 83mbs.AddLoad(Force(markerNumber = 2, loadVector = [0, -mass*g, 0]))
 84
 85
 86mbs.Assemble()
 87#exu.Print(mbs)
 88
 89simulationSettings = exu.SimulationSettings() #takes currently set values or default values
 90
 91f = 2000
 92simulationSettings.timeIntegration.numberOfSteps = 1*f
 93simulationSettings.timeIntegration.endTime = 0.001*f
 94simulationSettings.solutionSettings.writeSolutionToFile = True
 95simulationSettings.solutionSettings.solutionWritePeriod = simulationSettings.timeIntegration.endTime/500
 96simulationSettings.displayComputationTime = False
 97simulationSettings.displayStatistics = False
 98simulationSettings.timeIntegration.verboseMode = 1
 99
100simulationSettings.timeIntegration.newton.relativeTolerance = 1e-8*100 #10000
101simulationSettings.timeIntegration.newton.absoluteTolerance = 1e-10
102
103simulationSettings.timeIntegration.newton.useModifiedNewton = False
104simulationSettings.timeIntegration.newton.maxModifiedNewtonIterations = 5
105simulationSettings.timeIntegration.newton.numericalDifferentiation.minimumCoordinateSize = 1
106simulationSettings.timeIntegration.newton.numericalDifferentiation.relativeEpsilon = 6.055454452393343e-06*10 #eps^(1/3)
107simulationSettings.timeIntegration.generalizedAlpha.useIndex2Constraints = False
108simulationSettings.timeIntegration.generalizedAlpha.useNewmark = False
109simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.5
110
111#SC.visualizationSettings.nodes.showNumbers = True
112SC.visualizationSettings.bodies.showNumbers = True
113#SC.visualizationSettings.connectors.showNumbers = True
114SC.visualizationSettings.nodes.defaultSize = 0.05
115
116simulationSettings.solutionSettings.solutionInformation = "Planar four-bar-mechanism with initial angular velocity and gravity"
117
118#useGraphics = True #uncomment this line to visualize the example!
119if useGraphics:
120    exu.StartRenderer()
121    #mbs.WaitForUserToContinue()
122
123mbs.SolveDynamic(simulationSettings)
124
125if useGraphics:
126    SC.WaitForRenderEngineStopFlag()
127    exu.StopRenderer() #safely close rendering window!
128
129#compute error for test suite:
130sol = mbs.systemData.GetODE2Coordinates();
131u = sol[1]; #y-displacement of first node of four bar mechanism
132exu.Print('solution of fourbar mechanism =',u)
133
134exudynTestGlobals.testError = u - (-2.354666317492353) #2020-01-09: -2.354666317492353; 2019-12-15: (-2.3546596670554125); 2019-11-22:(-2.354659593986869);  previous: (-2.354659593986899)
135exudynTestGlobals.testResult = u