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