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