switchingConstraintsPendulum.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Rigid pendulum with tip mass
  5#           The functionality mbs.SetPreStepUserFunction is used to deactivate a constraint after a given time
  6#           Shows how to load solution to animate an existing solution
  7#
  8# Author:   Johannes Gerstmayr
  9# Date:     2019-11-22
 10#
 11# 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.
 12#
 13#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 14
 15import exudyn as exu
 16from exudyn.itemInterface import *
 17from exudyn.utilities import *
 18
 19SC = exu.SystemContainer()
 20mbs = SC.AddSystem()
 21
 22#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 23#rigid pendulum with initial velocities
 24#constraints and joints are deactivated during simulation
 25
 26rect = [-2.5,-1.5,0.5,1.5] #xmin,ymin,xmax,ymax
 27background = {'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
 28oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0], visualization=VObjectGround(graphicsData= [background])))
 29#nGround=mbs.AddNode(NodePointGround(referencePosition= [0,0,0]))
 30a = 0.5     #x-dim of pendulum
 31b = 0.05    #y-dim of pendulum
 32massRigid = 12
 33mass = 2 #of additional mass
 34inertiaRigid = massRigid/12*(2*a)**2
 35
 36omega0 = 4
 37
 38graphics2 = {'type':'Line', 'color':[0.1,0.1,0.8,1], 'data':[-a,-b,0, a,-b,0, a,b,0, -a,b,0, -a,-b,0]} #background
 39nRigid = mbs.AddNode(Rigid2D(referenceCoordinates=[-0.5,0,0], initialVelocities=[0,omega0*a,omega0]));
 40oRigid = mbs.AddObject(RigidBody2D(physicsMass=massRigid, physicsInertia=inertiaRigid,nodeNumber=nRigid,visualization=VObjectRigidBody2D(graphicsData= [graphics2])))
 41
 42mR1 = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid, localPosition=[-0.5,0.,0.])) #support point
 43mG0 = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oGround, localPosition=[-0.5-a,0.,0.]))
 44oRJoint = mbs.AddObject(RevoluteJoint2D(markerNumbers=[mG0,mR1],activeConnector=True))
 45
 46#mass point is attached with coordinate constraints:
 47mCoordR0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nRigid, coordinate=0))
 48mCoordR1 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nRigid, coordinate=1))
 49
 50#additional mass point attached to COM of rigid body:
 51nMass = mbs.AddNode(Point2D(referenceCoordinates=[-0.5,0], initialVelocities=[0,omega0*a]));
 52oMass = mbs.AddObject(MassPoint2D(physicsMass=mass, nodeNumber=nMass) )
 53
 54mCoordM0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nMass, coordinate=0))
 55mCoordM1 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nMass, coordinate=1))
 56
 57oConstraint0 = mbs.AddObject(CoordinateConstraint(markerNumbers=[mCoordM0, mCoordR0],activeConnector=True))
 58oConstraint1 = mbs.AddObject(CoordinateConstraint(markerNumbers=[mCoordM1, mCoordR1],activeConnector=True))
 59
 60mbs.Assemble()
 61
 62simulationSettings = exu.SimulationSettings() #takes currently set values or default values
 63simulationSettings.timeIntegration.numberOfSteps = 4000
 64simulationSettings.timeIntegration.endTime = 2
 65simulationSettings.timeIntegration.newton.relativeTolerance = 1e-8
 66simulationSettings.timeIntegration.newton.absoluteTolerance = 1e-6
 67simulationSettings.timeIntegration.verboseMode = 1
 68
 69def UFswitchConnector(mbs, t):
 70    if t > 0.3:
 71        mbs.SetObjectParameter(oRJoint, 'activeConnector', False)
 72    if t > 0.1:
 73        mbs.SetObjectParameter(oConstraint0, 'activeConnector', False)
 74        mbs.SetObjectParameter(oConstraint1, 'activeConnector', False)
 75    return True #True, means that everything is alright, False=stop simulation
 76
 77mbs.SetPreStepUserFunction(UFswitchConnector)
 78
 79simulationSettings.timeIntegration.newton.useModifiedNewton = False
 80simulationSettings.timeIntegration.newton.numericalDifferentiation.minimumCoordinateSize = 1
 81simulationSettings.timeIntegration.generalizedAlpha.useNewmark = True
 82simulationSettings.timeIntegration.generalizedAlpha.useIndex2Constraints = True
 83simulationSettings.solutionSettings.solutionInformation = "Rigid pendulum with switching constraints"
 84simulationSettings.displayStatistics = True
 85
 86SC.visualizationSettings.openGL.multiSampling = 1
 87
 88exu.StartRenderer()
 89
 90mbs.SolveDynamic(simulationSettings)
 91print('end time =',mbs.systemData.GetTime()) #time after time integration ...
 92#print('solution =',mbs.systemData.GetODE2Coordinates()) #solution coordinates after time integration ...
 93
 94mbs.WaitForUserToContinue()
 95
 96animate = True
 97if animate:
 98    fileName = 'coordinatesSolution.txt'
 99    solution = LoadSolutionFile('coordinatesSolution.txt')
100    AnimateSolution(mbs, solution, 1, 0.05)
101
102
103#SC.WaitForRenderEngineStopFlag()
104exu.StopRenderer() #safely close rendering window!