hydraulicActuatorSimpleTest.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  A one arm mechanism is actuated by the HydraulicActuatorSimple;
  5#           The actuator contains internal dynamics based on GenericODE1 node
  6#
  7# Author:   Johannes Gerstmayr
  8# Date:     2022-06-16
  9#
 10# 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.
 11#
 12#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 13
 14import exudyn as exu
 15from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
 16import exudyn.graphics as graphics
 17
 18testIsActive = exu.sys.get('testIsActive', False)
 19
 20import numpy as np
 21from math import sin, cos, sqrt,pi
 22
 23SC = exu.SystemContainer()
 24mbs = SC.AddSystem()
 25
 26L = 1    #x-dim of arm 
 27b = 0.1  #y-dim of arm 
 28
 29
 30#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 31#one arm mechanism
 32background = graphics.CheckerBoard(point=[0,0.5*L*0,-2*b],size=2)
 33oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0], visualization=VObjectGround(graphicsData= [background])))
 34massRigid = 12*10
 35inertiaRigid = massRigid/12*(L)**2
 36g = 9.81    # gravity
 37
 38graphicsList = [graphics.Brick(size= [L,b,0.1*b], color= graphics.color.dodgerblue, addEdges=True)]
 39
 40graphicsList += [graphics.Cylinder(pAxis=[-0.5*L,0,-0.7*b], vAxis= [0,0,1.4*b], radius = 0.55*b, 
 41                                     color= graphics.color.lightgrey, addEdges=True, nTiles=32)]
 42#print(graphicsList[2])
 43nRigid = mbs.AddNode(Rigid2D(referenceCoordinates=[0.5*L,0,0], initialVelocities=[0,0,0]));
 44oRigid = mbs.AddObject(RigidBody2D(mass=massRigid, inertia=inertiaRigid,nodeNumber=nRigid,
 45                                   visualization=VObjectRigidBody2D(graphicsData= graphicsList)))
 46
 47mR1 = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid, localPosition=[-0.5*L,0.,0.])) #support point
 48mR2 = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid, localPosition=[ 0.,0.,0.])) #end point
 49
 50#add joint
 51mG0 = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oGround, localPosition=[0,0,0]))
 52mbs.AddObject(RevoluteJoint2D(markerNumbers=[mG0,mR1]))
 53
 54mbs.AddLoad(Force(markerNumber = mR2, loadVector = [0, -massRigid*g, 0]))
 55
 56#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
 57#add hydraulics actuator:
 58mGH = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oGround, localPosition=[0,-0.25*L-0.5*b*0,0.]))
 59mRH = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid,  localPosition=[-0.25*L,-0.5*b*0,0.]))
 60
 61
 62LH0 = sqrt(2*(0.25*L)**2) #zero length of actuator
 63
 64#hydraulics parameters:
 65V0 = 1. #oil volume (could actually change ...)
 66V1 = V0 #oil volume (could actually change ...)
 67A=[0.01,0.01] #piston area side 1/2
 68Eoil = 1e11
 69Av1 = 1 #valve opening (factor)
 70Av2 = 0.0 #valve opening (factor)
 71Qn = 2e-5 #nominal flow
 72pS = 200.*1e5 #system pressure (200bar)
 73pT = 0.*1e5   #tank pressure; 
 74dampingHA = 2e5
 75
 76
 77#ODE1 for pressures:
 78nODE1 = mbs.AddNode(NodeGenericODE1(referenceCoordinates=[0,0],
 79                                    initialCoordinates=[2e6,2e6], #initialize with 20 bar
 80                                    numberOfODE1Coordinates=2))
 81
 82oHA = mbs.AddObject(HydraulicActuatorSimple(markerNumbers=[mGH, mRH],
 83                                            nodeNumbers=[nODE1],
 84                                            offsetLength=LH0, strokeLength=LH0*0.5, 
 85                                            chamberCrossSection0=A[0], chamberCrossSection1=A[1],
 86                                            hoseVolume0=V0, hoseVolume1=V1,
 87                                            valveOpening0=0, valveOpening1=0,
 88                                            oilBulkModulus=Eoil, actuatorDamping=dampingHA, nominalFlow=Qn, 
 89                                            systemPressure=pS, tankPressure=pT,
 90                                            useChamberVolumeChange=False,
 91                                            visualization=VHydraulicActuatorSimple(cylinderRadius= 0.6*b, rodRadius= 0.3*b,
 92                                                                                   baseMountLength = 0.4*b, baseMountRadius = 0.4*b,
 93                                                                                   rodMountRadius = 0.3*b, pistonLength = 0.2*b, pistonRadius = 0.55*b,
 94                                                                                   colorCylinder=graphics.color.blue, colorPiston=graphics.color.lightgrey),
 95                                            )) 
 96
 97
 98def PreStepUserFunction(mbs, t):
 99    LHact = mbs.GetObjectOutput(oHA, variableType=exu.OutputVariableType.Distance)
100    x = (max(0.5, min(1.5,(1-cos(t*pi*2*0.5))) ) - 0.5)*0.1+LH0
101    #if t>2: x=LH0
102
103    Av0 = (x-LHact)*2 #valve position control ==> penalize set value LH0
104    #print('Av0=',Av0)
105    Av1 = -Av0
106    mbs.SetObjectParameter(oHA, "valveOpening0", Av0)
107    mbs.SetObjectParameter(oHA, "valveOpening1", Av1)
108    return True
109
110mbs.SetPreStepUserFunction(PreStepUserFunction)
111
112
113sForce = mbs.AddSensor(SensorObject(objectNumber=oHA, storeInternal=True, outputVariableType=exu.OutputVariableType.Force))
114sDistance = mbs.AddSensor(SensorObject(objectNumber=oHA, storeInternal=True, outputVariableType=exu.OutputVariableType.Distance))
115sVelocity = mbs.AddSensor(SensorObject(objectNumber=oHA, storeInternal=True, outputVariableType=exu.OutputVariableType.Velocity))
116sPressures = mbs.AddSensor(SensorNode(nodeNumber=nODE1, storeInternal=True, outputVariableType=exu.OutputVariableType.Coordinates))
117
118mbs.Assemble()
119
120#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
121
122simulationSettings = exu.SimulationSettings() #takes currently set values or default values
123
124
125tEnd = 0.4
126stepSize = 1e-3
127simulationSettings.timeIntegration.numberOfSteps = int(tEnd/stepSize)
128simulationSettings.timeIntegration.endTime = tEnd
129simulationSettings.timeIntegration.startTime = 0
130simulationSettings.timeIntegration.newton.relativeTolerance = 1e-8*100 #10000
131simulationSettings.timeIntegration.newton.absoluteTolerance = 1e-10
132simulationSettings.timeIntegration.verboseMode = 1
133# simulationSettings.timeIntegration.realtime.active = True #to see what happens ...
134
135simulationSettings.timeIntegration.newton.useModifiedNewton = True
136simulationSettings.timeIntegration.newton.numericalDifferentiation.minimumCoordinateSize = 1
137simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.5
138simulationSettings.show.statistics = True
139
140simulationSettings.solution.file.information = 'Hydraulics user function test'
141
142SC.visualizationSettings.openGL.multiSampling = 4
143SC.visualizationSettings.openGL.lineWidth = 2
144
145if not testIsActive:
146    SC.renderer.Start()
147    SC.renderer.DoIdleTasks()
148
149mbs.SolveDynamic(simulationSettings, showHints=False)
150
151if not testIsActive:
152    SC.renderer.DoIdleTasks()
153    SC.renderer.Stop() #safely close rendering window!
154
155exu.Print('hydraulics C++:')
156exu.Print('pressures=', mbs.GetSensorValues(sPressures))
157exu.Print('velocity=', mbs.GetSensorValues(sVelocity))
158#for stepSize=1e-6: error about 1e-5 compared to user function implementation; with initialVelocities=[0,0,2] and tEnd=0.4
159# hydraulics C++:
160# pressures= [6441296.09086297 3008420.04232005]
161# velocity= [-0.0050061   0.20338669  0.        ]
162
163# 
164# mbs.PlotSensor(sensorNumbers=sForce, components=exudyn.plot.componentNorm, labels=['connector force norm'], yLabel='force (N)', closeAll=True)
165# mbs.PlotSensor(sensorNumbers=sDistance, components=0)
166# mbs.PlotSensor(sensorNumbers=[sPressures]*2, components=[0,1], labels=['p1', 'p2'], yLabel='pressure (N/m^2)')
167
168#mbs.PlotSensor(sensorNumbers=p01, components=0, labels=['differential hydraulic force'], yLabel='hydraulic force (N)')
169
170#compute error for test suite:
171sol2 = mbs.systemData.GetODE2Coordinates(); 
172sol1 = mbs.systemData.GetODE1Coordinates(); 
173u = np.linalg.norm(sol2); 
174u += np.linalg.norm(sol1)*1e-6;
175exu.Print('solution of hydraulicActuatorSimpleTest =',u)
176
177exu.sys['testResult'] = u