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
 16testIsActive = exu.sys.get('testIsActive', False)
 17
 18SC = exu.SystemContainer()
 19mbs = SC.AddSystem()
 20
 21nBodies = 4
 22mass = 3 #kg
 23g = 9.81 # m/s^2
 24omega0 = 2 #initial angular velocity of first link
 25
 26#add nodes:
 27n0=mbs.AddNode(PointGround(referenceCoordinates=[0, 0, 0]))
 28n1=mbs.AddNode(Point(referenceCoordinates=[0, 2, 0], initialVelocities= [2*omega0,0,0]))
 29n2=mbs.AddNode(Point(referenceCoordinates = [4, 5, 0], initialVelocities= [5*omega0,0,0]))
 30n3=mbs.AddNode(Point(referenceCoordinates = [4, 0, 0]))
 31
 32#nodetypes Point2DSlope1Slope2 Point2DSlope1 RigidBody RigidBody2D
 33
 34L=4
 35Lo=0.
 36#add mass points and ground object:
 37rect = [-6,-6,6,6] #xmin,ymin,xmax,ymax
 38background = {'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
 39#oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0], visualization=VObjectGround(graphicsData= [background])))
 40#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
 41o0=mbs.AddObject(MassPoint(mass=0,nodeNumber=n0,visualization=VObjectMassPoint(graphicsData= [background])))
 42o1=mbs.AddObject(MassPoint(mass=mass,nodeNumber=n1))
 43o2=mbs.AddObject(MassPoint(mass=mass,nodeNumber=n2))
 44o3=mbs.AddObject(MassPoint(mass=mass,nodeNumber=n3))
 45
 46#use NodePosition markers:
 47for i in range(nBodies): 
 48    mbs.AddMarker(MarkerNodePosition(nodeNumber = i))
 49
 50#Coordinate Marker:
 51mc0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = n0, coordinate=0)) #Ground node ==> no action
 52mc1 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = n3, coordinate=1))
 53mc2 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = n3, coordinate=0))
 54mbs.AddObject(CoordinateConstraint(markerNumbers=[mc0,mc1]))
 55mbs.AddObject(CoordinateConstraint(markerNumbers=[mc0,mc2]))
 56
 57useConstraint = True
 58if useConstraint:
 59    #add distance:
 60    mbs.AddObject(ObjectConnectorDistance(markerNumbers=[0,1], distance=2, visualization=VObjectConnectorDistance(drawSize=0.01)))
 61    mbs.AddObject(ObjectConnectorDistance(markerNumbers=[1,2], distance=5, visualization=VObjectConnectorDistance(drawSize=0.01)))
 62    mbs.AddObject(ObjectConnectorDistance(markerNumbers=[2,3], distance=5, visualization=VObjectConnectorDistance(drawSize=0.01)))
 63else:
 64    #add spring-dampers:
 65    k = 40000
 66    d = 20
 67    mbs.AddObject({'objectType': 'ConnectorSpringDamper', 'stiffness': k, 'damping': d, 'force': 0, 'referenceLength':2, 'markerNumbers': [0,1], 'VdrawSize': 0.01})
 68    mbs.AddObject({'objectType': 'ConnectorSpringDamper', 'stiffness': k, 'damping': d, 'force': 0, 'referenceLength':5, 'markerNumbers': [1,2], 'VdrawSize': 0.01})
 69    mbs.AddObject({'objectType': 'ConnectorSpringDamper', 'stiffness': k, 'damping': d, 'force': 0, 'referenceLength':5, 'markerNumbers': [2,3], 'VdrawSize': 0.01})
 70
 71#add loads:
 72mbs.AddLoad(Force(markerNumber = 1, loadVector = [0, -mass*g, 0]))
 73mbs.AddLoad(Force(markerNumber = 2, loadVector = [0, -mass*g, 0]))
 74
 75
 76mbs.Assemble()
 77#exu.Print(mbs)
 78
 79simulationSettings = exu.SimulationSettings() #takes currently set values or default values
 80
 81f = 2000
 82simulationSettings.timeIntegration.numberOfSteps = 1*f
 83simulationSettings.timeIntegration.endTime = 0.001*f
 84simulationSettings.solution.file.write = (not testIsActive) #only the SolutionViewer reads it (#2492)
 85simulationSettings.solution.file.writePeriod = simulationSettings.timeIntegration.endTime/500
 86simulationSettings.show.computationTime = False
 87simulationSettings.show.statistics = False
 88simulationSettings.timeIntegration.verboseMode = 1
 89
 90simulationSettings.timeIntegration.newton.relativeTolerance = 1e-8*100 #10000
 91simulationSettings.timeIntegration.newton.absoluteTolerance = 1e-10
 92
 93simulationSettings.timeIntegration.newton.useModifiedNewton = False
 94simulationSettings.timeIntegration.newton.maxModifiedNewtonIterations = 5
 95simulationSettings.timeIntegration.newton.numericalDifferentiation.minimumCoordinateSize = 1
 96simulationSettings.timeIntegration.newton.numericalDifferentiation.relativeEpsilon = 6.055454452393343e-06*10 #eps^(1/3)
 97simulationSettings.timeIntegration.generalizedAlpha.useIndex2Constraints = False
 98simulationSettings.timeIntegration.generalizedAlpha.useNewmark = False
 99simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.5
100
101#SC.visualizationSettings.nodes.showNumbers = True
102SC.visualizationSettings.bodies.showNumbers = True
103#SC.visualizationSettings.connectors.showNumbers = True
104SC.visualizationSettings.nodes.defaultSize = 0.05
105
106simulationSettings.solution.file.information = "Planar four-bar-mechanism with initial angular velocity and gravity"
107
108#testIsActive = False #uncomment this line to visualize the example!
109if not testIsActive: 
110    SC.renderer.Start()
111    #SC.renderer.DoIdleTasks()
112
113mbs.SolveDynamic(simulationSettings)
114
115if not testIsActive: 
116    SC.renderer.DoIdleTasks()
117    SC.renderer.Stop() #safely close rendering window!
118
119#compute error for test suite:
120sol = mbs.systemData.GetODE2Coordinates(); 
121u = sol[1]; #y-displacement of first node of four bar mechanism
122exu.Print('solution of fourbar mechanism =',u)
123
124exu.sys['testResult'] = u