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