HydraulicActuatorStaticInitialization.py
You can view and download this file on Github: HydraulicActuatorStaticInitialization.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: A one arm mechanism is actuated by the HydraulicActuatorSimple;
5# This particular example shows how a static computation can be performed with the hydraulic actuator;
6# For static computation, a distance constraint is used to replace the hydraulic actuator;
7# Hereafter, the dynamic simulation is initialized with the static equilibrium; this can be used for flexible booms
8#
9# Author: Johannes Gerstmayr
10# Date: 2023-09-07
11#
12# 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.
13#
14#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
15
16import exudyn as exu
17from exudyn.utilities import *
18
19useGraphics = True #without test
20
21import numpy as np
22from math import sin, cos, sqrt,pi
23
24SC = exu.SystemContainer()
25mbs = SC.AddSystem()
26
27L = 1 #x-dim of arm
28b = 0.1 #y-dim of arm
29
30
31#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
32#one arm mechanism
33background = GraphicsDataCheckerBoard(point=[0,0.5*L*0,-2*b],size=2)
34oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0], visualization=VObjectGround(graphicsData= [background])))
35massRigid = 12*10
36inertiaRigid = massRigid/12*(L)**2
37g = 9.81 # gravity
38
39graphicsList = [GraphicsDataOrthoCubePoint(size= [L,b,0.1*b], color= color4dodgerblue, addEdges=True)]
40
41graphicsList += [GraphicsDataCylinder(pAxis=[-0.5*L,0,-0.7*b], vAxis= [0,0,1.4*b], radius = 0.55*b,
42 color= color4lightgrey, addEdges=True, nTiles=32)]
43#print(graphicsList[2])
44nRigid = mbs.AddNode(Rigid2D(referenceCoordinates=[0.5*L,0,0], initialVelocities=[0,0,0]));
45oRigid = mbs.AddObject(RigidBody2D(physicsMass=massRigid, physicsInertia=inertiaRigid,nodeNumber=nRigid,
46 visualization=VObjectRigidBody2D(graphicsData= graphicsList)))
47
48mR1 = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid, localPosition=[-0.5*L,0.,0.])) #support point
49mR2 = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid, localPosition=[ 0.,0.,0.])) #end point
50
51#add joint
52mG0 = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oGround, localPosition=[0,0,0]))
53mbs.AddObject(RevoluteJoint2D(markerNumbers=[mG0,mR1]))
54
55mbs.AddLoad(Force(markerNumber = mR2, loadVector = [0, -massRigid*g, 0]))
56
57#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
58#add hydraulics actuator:
59mGH = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oGround, localPosition=[0,-0.25*L-0.5*b*0,0.]))
60mRH = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid, localPosition=[-0.25*L,-0.5*b*0,0.]))
61
62
63LH0 = sqrt(2*(0.25*L)**2) #zero length of actuator
64
65#hydraulics parameters:
66V0 = 1. #oil volume (could actually change ...)
67V1 = V0 #oil volume (could actually change ...)
68A=[0.01,0.01] #piston area side 1/2
69Eoil = 1e11
70Av1 = 1 #valve opening (factor)
71Av2 = 0.0 #valve opening (factor)
72Qn = 2e-5 #nominal flow
73pS = 200.*1e5 #system pressure (200bar)
74pT = 0.*1e5 #tank pressure;
75dampingHA = 2e5
76
77
78useHydraulics=True
79staticInitialization=True
80if useHydraulics:
81 #ODE1 for pressures:
82 nODE1 = mbs.AddNode(NodeGenericODE1(referenceCoordinates=[0,0],
83 initialCoordinates=[2e6,2e6], #initialize with 20 bar
84 numberOfODE1Coordinates=2))
85 oHA = mbs.AddObject(HydraulicActuatorSimple(markerNumbers=[mGH, mRH],
86 nodeNumbers=[nODE1],
87 offsetLength=LH0, strokeLength=LH0*0.5,
88 chamberCrossSection0=A[0], chamberCrossSection1=A[1],
89 hoseVolume0=V0, hoseVolume1=V1,
90 valveOpening0=0, valveOpening1=0,
91 oilBulkModulus=Eoil, actuatorDamping=dampingHA, nominalFlow=Qn,
92 systemPressure=pS, tankPressure=pT,
93 useChamberVolumeChange=False,
94 visualization=VHydraulicActuatorSimple(cylinderRadius= 0.6*b, rodRadius= 0.3*b,
95 baseMountLength = 0.4*b, baseMountRadius = 0.4*b,
96 rodMountRadius = 0.3*b, pistonLength = 0.2*b, pistonRadius = 0.55*b,
97 colorCylinder=color4blue, colorPiston=color4lightgrey),
98 ))
99
100 def PreStepUserFunction(mbs, t):
101 LHact = mbs.GetObjectOutput(oHA, variableType=exu.OutputVariableType.Distance)
102 x = (max(0.5, min(1.5,(1-cos(t*pi*2*0.5))) ) - 0.5)*0.1+LH0
103 #if t>2: x=LH0
104
105 Av0 = (x-LHact)*2 #valve position control ==> penalize set value LH0
106 #print('Av0=',Av0)
107 Av1 = -Av0
108 mbs.SetObjectParameter(oHA, "valveOpening0", Av0)
109 mbs.SetObjectParameter(oHA, "valveOpening1", Av1)
110 return True
111
112
113 sForce = mbs.AddSensor(SensorObject(objectNumber=oHA, storeInternal=True, outputVariableType=exu.OutputVariableType.Force))
114 sDistance = mbs.AddSensor(SensorObject(objectNumber=oHA, storeInternal=True, outputVariableType=exu.OutputVariableType.Distance))
115 sVelocity = mbs.AddSensor(SensorObject(objectNumber=oHA, storeInternal=True, outputVariableType=exu.OutputVariableType.Velocity))
116 sPressures = mbs.AddSensor(SensorNode(nodeNumber=nODE1, storeInternal=True, outputVariableType=exu.OutputVariableType.Coordinates))
117
118#compute reference length of distance constraint (this is LH0 in this case, but could be else):
119mGHposition = mbs.GetMarkerOutput(mGH, variableType=exu.OutputVariableType.Position,
120 configuration=exu.ConfigurationType.Reference)
121mRHposition = mbs.GetMarkerOutput(mRH, variableType=exu.OutputVariableType.Position,
122 configuration=exu.ConfigurationType.Reference)
123
124dLH0 = NormL2(mGHposition - mRHposition)
125# print('LH0=', LH0)
126# print('dLH0=', dLH0)
127
128#use distance constraint to compute static equlibrium in static case
129oDC = mbs.AddObject(DistanceConstraint(markerNumbers=[mGH, mRH],
130 distance=dLH0))
131
132mbs.Assemble()
133
134#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
135
136simulationSettings = exu.SimulationSettings() #takes currently set values or default values
137
138
139tEnd = 1
140stepSize = 1e-3
141simulationSettings.timeIntegration.numberOfSteps = int(tEnd/stepSize)
142simulationSettings.timeIntegration.endTime = tEnd
143simulationSettings.timeIntegration.startTime = 0
144simulationSettings.timeIntegration.newton.relativeTolerance = 1e-8*100 #10000
145simulationSettings.timeIntegration.newton.absoluteTolerance = 1e-10
146simulationSettings.timeIntegration.verboseMode = 1
147# simulationSettings.timeIntegration.simulateInRealtime = True #to see what happens ...
148
149simulationSettings.timeIntegration.newton.useModifiedNewton = True
150simulationSettings.timeIntegration.newton.numericalDifferentiation.minimumCoordinateSize = 1
151simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.5
152simulationSettings.displayStatistics = True
153
154simulationSettings.solutionSettings.solutionInformation = 'Hydraulics user function test'
155
156SC.visualizationSettings.openGL.multiSampling = 4
157SC.visualizationSettings.openGL.lineWidth = 2
158
159if useGraphics:
160 exu.StartRenderer()
161 # mbs.WaitForUserToContinue()
162
163simulationSettings.staticSolver.constrainODE1coordinates = True #True: set pressures to initial values
164if staticInitialization:
165 exu.SolveStatic(mbs, simulationSettings, updateInitialValues=True) #results are new initial values
166 force = mbs.GetObjectOutput(oDC, variableType=exu.OutputVariableType.Force)
167 print('initial force=', force)
168
169mbs.SetObjectParameter(oDC, 'activeConnector', False)
170if useHydraulics:
171 if staticInitialization:
172 p0 = 2e6 + force/A[0]
173 p1 = 2e6
174
175 #now we would like to reset the pressures:
176 #1) chance initial in NodeGenericODE1 => then mbs.Assemble() => this would destroy the previously computed initial values
177 #2) change the initial values in the system vector
178
179 sysODE1 = mbs.systemData.GetODE1Coordinates(configuration=exu.ConfigurationType.Initial)
180 nODE1index = mbs.GetNodeODE1Index(nODE1)
181 print('sysODE1=',sysODE1)
182 print('p0,p1=',p0,p1)
183 sysODE1[nODE1index] = p0
184 sysODE1[nODE1index+1] = p1
185
186
187 #now write the updated system variables:
188 mbs.systemData.SetODE1Coordinates(coordinates=sysODE1, configuration=exu.ConfigurationType.Initial)
189
190 #mbs.SetObjectParameter(oHA, '')
191 mbs.SetPreStepUserFunction(PreStepUserFunction)
192 exu.SolveDynamic(mbs, simulationSettings, showHints=False)
193
194if useGraphics:
195 SC.WaitForRenderEngineStopFlag()
196 exu.StopRenderer() #safely close rendering window!
197
198if useHydraulics:
199 exu.Print('hydraulics C++:')
200 exu.Print('pressures=', mbs.GetSensorValues(sPressures))
201 exu.Print('velocity=', mbs.GetSensorValues(sVelocity))
202 #for stepSize=1e-6: error about 1e-5 compared to user function implementation; with initialVelocities=[0,0,2] and tEnd=0.4
203 # hydraulics C++:
204 # pressures= [6441296.09086297 3008420.04232005]
205 # velocity= [-0.0050061 0.20338669 0. ]
206
207 # from exudyn.plot import PlotSensor
208 # PlotSensor(mbs, sensorNumbers=sForce, components=exudyn.plot.componentNorm, labels=['connector force norm'], yLabel='force (N)', closeAll=True)
209 # PlotSensor(mbs, sensorNumbers=sDistance, components=0)
210 mbs.PlotSensor(sensorNumbers=[sPressures]*2, components=[0,1], labels=['p0', 'p1'], yLabel='pressure (N/m^2)')
211
212 #PlotSensor(mbs, sensorNumbers=p01, components=0, labels=['differential hydraulic force'], yLabel='hydraulic force (N)')
213
214 #compute error for test suite:
215 sol2 = mbs.systemData.GetODE2Coordinates();
216 sol1 = mbs.systemData.GetODE1Coordinates();
217 u = np.linalg.norm(sol2);
218 u += np.linalg.norm(sol1)*1e-6;
219 exu.Print('solution of hydraulicActuatorSimpleTest =',u)