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