contactCoordinateTest.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Test with ObjectContactCoordinate, which can be used to achieve accurate contact simulation
  5#           Uses step reduction to resolve contact switching point
  6#           Similar to postNewtonStepContactTest but with stiffer contact
  7#
  8# Author:   Johannes Gerstmayr
  9# Date:     2021-08-12
 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.utilities import * #includes itemInterface and rigidBodyUtilities
 17import exudyn.graphics as graphics #only import if it does not conflict
 18
 19useGraphics = True #without test
 20#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 21#you can erase the following lines and all exudynTestGlobals related operations if this is not intended to be used as TestModel:
 22try: #only if called from test suite
 23    from modelUnitTests import exudynTestGlobals #for globally storing test results
 24    useGraphics = exudynTestGlobals.useGraphics
 25except:
 26    class ExudynTestGlobals:
 27        pass
 28    exudynTestGlobals = ExudynTestGlobals()
 29#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 30
 31SC = exu.SystemContainer()
 32mbs = SC.AddSystem()
 33
 34withUserFunction = False #compare to user function based test
 35
 36#define parameters of mass point
 37L=0.5
 38r = 0.05
 39g=9.81
 40mass = 0.25         #mass in kg
 41spring = 20000*100        #stiffness of spring-damper in N/m
 42damper = 0.0001*spring          #damping constant in N/(m/s)
 43load0 = -mass*g     #in negative y-direction
 44
 45doRefSol = False
 46tEnd = 0.25        #end time of simulation
 47h = 2e-3*0.1
 48if doRefSol:
 49    h=1e-5
 50
 51#data coordinate: contains gap for ObjectContactCoordinate, for user function: 1=no contact, 0=contact
 52nData=mbs.AddNode(NodeGenericData(initialCoordinates=[1], numberOfDataCoordinates=1))
 53
 54#node for 3D mass point:
 55n1=mbs.AddNode(Point(referenceCoordinates = [0,0,0],
 56                     initialCoordinates = [0,r+(L*0+0.05),0],
 57                     initialVelocities = [0,-0.1*0,0]))
 58
 59#user function for spring force
 60def springForce(mbs, t, itemIndex, u, v, k, d, offset, mu, muPropZone):
 61    p = 0*L+u-r
 62    data = mbs.systemData.GetDataCoordinates()
 63    #print("p=", p, ", contact=", data[0])
 64    if data[0] == 0:
 65        return k*p + d*v
 66    else:
 67        return 0
 68
 69def PostNewtonUserFunction(mbs, t):
 70    p0 = mbs.GetNodeOutput(n1, exu.OutputVariableType.Position, configuration=exu.ConfigurationType.StartOfStep)[1] - r
 71    p = mbs.GetNodeOutput(n1, exu.OutputVariableType.Position)[1] - r
 72    #v0 = mbs.GetNodeOutput(n1, exu.OutputVariableType.Velocity, configuration=exu.ConfigurationType.StartOfStep)[1]
 73    #a0 = mbs.GetNodeOutput(n1, exu.OutputVariableType.Acceleration, configuration=exu.ConfigurationType.StartOfStep)[1]
 74    h = mbs.sys['dynamicSolver'].it.currentStepSize #grab current step size from dynamic solver
 75    data = mbs.systemData.GetDataCoordinates()
 76    data0 = mbs.systemData.GetDataCoordinates(configuration=exu.ConfigurationType.StartOfStep)
 77
 78    #data[0] = 0 #no contact; 0 corresponds to the only one data coordinate in the system, attributed to contact
 79    recommendedStepSize = -1
 80    error = 0
 81    #check if previous assumption was wrong ==> set error, reduce step size and set new contact state
 82    if p < 0:
 83        if data0[0] == 1:
 84            error = abs(p)
 85
 86            if (p0 > 0):
 87                recommendedStepSize = h*(abs(p0))/(abs(p0)+abs(p))
 88            else:
 89                recommendedStepSize = 0.25*h #simple alternative
 90
 91            data[0] = 0 #contact
 92
 93    else:
 94        if data0[0] == 0:
 95            error = abs(p)
 96            #recommendedStepSize = 1e-6 #simple alternative
 97            if (p0 > 0):
 98                recommendedStepSize = h*(abs(p0))/(abs(p0)+abs(p))
 99            else:
100                recommendedStepSize = 0.25*h #simple alternative
101            data[0] = 1 #contact off
102
103    mbs.systemData.SetDataCoordinates(data)
104    return [error,recommendedStepSize]
105
106sMode = ""
107if withUserFunction:
108    mbs.SetPostNewtonUserFunction(PostNewtonUserFunction)
109    sMode = "User"
110
111#ground node
112d=0.01
113gGround = graphics.Brick([0,-d*0.5,0],[2*L,d,d],color=graphics.color.grey)
114oGround=mbs.AddObject(ObjectGround(visualization=VObjectGround(graphicsData=[gGround])))
115
116nGround=mbs.AddNode(NodePointGround(referenceCoordinates = [0,0,0]))
117
118#add mass point (this is a 3D object with 3 coordinates):
119gSphere = graphics.Sphere([0,0,0], r, color=graphics.color.red, nTiles=20)
120massPoint = mbs.AddObject(MassPoint(physicsMass = mass, nodeNumber = n1,
121                                    visualization=VMassPoint(graphicsData=[gSphere])))
122
123#marker for ground (=fixed):
124groundMarker=mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= nGround, coordinate = 0))
125#marker for springDamper for first (x-)coordinate:
126nodeMarker  =mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= n1, coordinate = 1)) #y-coordinate
127
128#Spring-Damper between two marker coordinates
129if withUserFunction:
130    sensorFileName='solution/sensorPos'+sMode+'.txt'
131    mbs.AddObject(CoordinateSpringDamper(markerNumbers = [groundMarker, nodeMarker],
132                                         stiffness = spring, damping = damper,
133                                         springForceUserFunction = springForce,
134                                         visualization=VCoordinateSpringDamper(show=False)))
135else:
136    sensorFileName=''
137    mbs.AddObject(ObjectContactCoordinate(markerNumbers = [groundMarker, nodeMarker],
138                                          nodeNumber = nData,
139                                          contactStiffness = spring, contactDamping = damper,
140                                          offset = r,
141                                          visualization=VObjectContactCoordinate(show=False)))
142
143#add load:
144loadC = mbs.AddLoad(LoadCoordinate(markerNumber = nodeMarker,
145                           load = load0))
146
147
148if useGraphics:
149    sPos = mbs.AddSensor(SensorNode(nodeNumber=n1, outputVariableType=exu.OutputVariableType.Position,
150                                    storeInternal=True,fileName=sensorFileName
151                                    ))
152
153mbs.Assemble()
154
155simulationSettings = exu.SimulationSettings()
156simulationSettings.solutionSettings.writeSolutionToFile = False
157simulationSettings.solutionSettings.sensorsWritePeriod = 1e-10
158simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
159simulationSettings.timeIntegration.endTime = tEnd
160simulationSettings.timeIntegration.minimumStepSize = 1e-10
161simulationSettings.timeIntegration.stepInformation = 3 #do not show step increase
162
163#important settings for contact:
164simulationSettings.timeIntegration.verboseMode = 1
165simulationSettings.timeIntegration.newton.useModifiedNewton = False #True=does not work yet
166simulationSettings.timeIntegration.discontinuous.iterationTolerance = 1e-8 #this is the accepted penetration before reducing step size
167if not withUserFunction:
168    simulationSettings.timeIntegration.discontinuous.iterationTolerance = 1e-3 #this is the accepted contact force error before reducing step size
169
170simulationSettings.timeIntegration.discontinuous.maxIterations = 2 #1=immediately perform step reduction
171simulationSettings.timeIntegration.discontinuous.ignoreMaxIterations = False #repeat step in case of failure
172simulationSettings.timeIntegration.adaptiveStepRecoverySteps = 0 #number of steps to wait until step size is increased again
173simulationSettings.timeIntegration.adaptiveStepIncrease = 10    #after successful step, increase again rapidly
174
175simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1 #for index 3 solver, this would be the best case
176
177simulationSettings.displayStatistics = True
178#simulationSettings.timeIntegration.simulateInRealtime = True
179
180if useGraphics:
181    exu.StartRenderer()              #start graphics visualization
182    mbs.WaitForUserToContinue()    #wait for pressing SPACE bar to continue
183
184#start solver:
185mbs.SolveDynamic(solverType=exu.DynamicSolverType.TrapezoidalIndex2, simulationSettings=simulationSettings) #chose index2, can handle adaptive steps
186#mbs.SolveDynamic(solverType=exu.DynamicSolverType.RK67, simulationSettings=simulationSettings)
187
188if useGraphics:
189    exu.StopRenderer()               #safely close rendering window!
190
191u = mbs.GetNodeOutput(n1, exu.OutputVariableType.Position)
192exu.Print('contactCoordinateTest=',u[1])
193
194exudynTestGlobals.testError = u[1] - (0.055313199503736685) #2021-08-13: 0.055313199503736685 (may change significantly for other disc. solver strategies)
195exudynTestGlobals.testResult = u[1]
196
197#+++++++++++++++++++++++++++++++++++++++++++++++++++++
198
199if useGraphics and True: #to run this, run model first with withUserFunction=True
200
201    mbs.PlotSensor(sensorNumbers=[sPos, 'solution/sensorPosUser.txt'], components=1,
202               labels=['internal contact','user function'])