postNewtonStepContactTest.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Test with postNewtonUserFunction and recommendedStepSize modeling a simplistic 1-mass- penalty contact problem;
  5#           Uses step reduction to resolve contact switching point
  6#
  7# Author:   Johannes Gerstmayr
  8# Date:     2021-03-20
  9#
 10# 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.
 11#
 12#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 13
 14import exudyn as exu
 15from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
 16import exudyn.graphics as graphics
 17
 18testIsActive = exu.sys.get('testIsActive', False)
 19
 20SC = exu.SystemContainer()
 21mbs = SC.AddSystem()
 22exu.Print('EXUDYN version='+exu.config.Version())
 23
 24#define parameters of mass point
 25L=0.5
 26r = 0.05
 27g=9.81
 28mass = 0.25         #mass in kg
 29spring = 20000        #stiffness of spring-damper in N/m
 30damper = 0*0.01*spring          #damping constant in N/(m/s)
 31load0 = -mass*g     #in negative y-direction
 32
 33doRefSol = False
 34tEnd = 0.5     #end time of simulation
 35h = 5e-3
 36if doRefSol:
 37    h=1e-5
 38
 39#data coordinate: 0=no contact, 1=contact
 40nData=mbs.AddNode(NodeGenericData(initialCoordinates=[0], numberOfDataCoordinates=1))
 41
 42#node for 3D mass point:
 43n1=mbs.AddNode(Point(referenceCoordinates = [0,L,0],
 44                     initialCoordinates = [0,r-L+0.01,0],
 45                     initialVelocities = [0,-1,0]))
 46
 47#user function for spring force
 48def springForce(mbs, t, itemIndex, u, v, k, d, offset): #changed 2023-01-21:, mu, muPropZone):
 49    p = L+u-r
 50    #print(p)
 51    data = mbs.systemData.GetDataCoordinates()
 52    if data[0] == 1:
 53        return k*p
 54    else:
 55        return 0
 56
 57def PostNewtonUserFunction(mbs, t):
 58    p0 = mbs.GetNodeOutput(n1, exu.OutputVariableType.Position, configuration=exu.ConfigurationType.StartOfStep)[1] - r
 59    p = mbs.GetNodeOutput(n1, exu.OutputVariableType.Position)[1] - r
 60    #v0 = mbs.GetNodeOutput(n1, exu.OutputVariableType.Velocity, configuration=exu.ConfigurationType.StartOfStep)[1]
 61    #a0 = mbs.GetNodeOutput(n1, exu.OutputVariableType.Acceleration, configuration=exu.ConfigurationType.StartOfStep)[1]
 62    h = mbs.sys['dynamicSolver'].it.currentStepSize #grab current step size from dynamic solver
 63    data = mbs.systemData.GetDataCoordinates()
 64    data0 = mbs.systemData.GetDataCoordinates(configuration=exu.ConfigurationType.StartOfStep)
 65
 66    #data[0] = 0 #no contact; 0 corresponds to the only one data coordinate in the system, attributed to contact
 67    recommendedStepSize = -1
 68    error = 0
 69    #check if previous assumption was wrong ==> set error, reduce step size and set new contact state
 70    if p < 0:
 71        if data0[0] == 0:
 72            error = abs(p)
 73            #recommendedStepSize = 1e-6 #simple alternative
 74            #x = abs(v0)*h #this is the estimated distance (acc=0) per step
 75            #x = 0.5*abs(a0)*h**2
 76            #recommendedStepSize = min(h,abs(h*(x-error)/x)) #assuming almost constant velocity during step
 77            if (p0 > 0):
 78                recommendedStepSize = h*(abs(p0))/(abs(p0)+abs(p))
 79            else:
 80                recommendedStepSize = 0.25*h #simple alternative
 81
 82
 83            data[0] = 1 #contact
 84        #else:
 85        #    recommendedStepSize = 1e-4
 86        #    error = abs(h-1e-4)
 87    else:
 88        if data0[0] == 1:
 89            error = abs(p)
 90            #recommendedStepSize = 1e-6 #simple alternative
 91            if (p0 > 0):
 92                recommendedStepSize = h*(abs(p0))/(abs(p0)+abs(p))
 93            else:
 94                recommendedStepSize = 0.25*h #simple alternative
 95            data[0] = 0 #contact off
 96
 97    #print("t=", round(t,6), ", p=", round(p,6), ", p0=", round(p0,6), #", a0=", round(a0,6), 
 98    #      ", h=", round(h,6), ", hRec=", 
 99    #      round(recommendedStepSize,6), ", tRec=", round(t-h+recommendedStepSize,6), 
100    #      ", c0=", data0[0], ", c=", data[0], ", e=", error)
101
102    mbs.systemData.SetDataCoordinates(data)
103    return [error,recommendedStepSize]
104
105mbs.SetPostNewtonUserFunction(PostNewtonUserFunction)
106
107#ground node
108d=0.01
109gGround = graphics.Brick([0,-d*0.5,0],[2*L,d,d],color=graphics.color.grey)
110oGround=mbs.AddObject(ObjectGround(visualization=VObjectGround(graphicsData=[gGround])))
111
112nGround=mbs.AddNode(NodePointGround(referenceCoordinates = [0,0,0]))
113
114#add mass point (this is a 3D object with 3 coordinates):
115gSphere = graphics.Sphere([0,0,0], r, color=graphics.color.red, nTiles=20)
116massPoint = mbs.AddObject(MassPoint(mass = mass, nodeNumber = n1,
117                                    visualization=VMassPoint(graphicsData=[gSphere])))
118
119#marker for ground (=fixed):
120groundMarker=mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= nGround, coordinate = 0))
121#marker for springDamper for first (x-)coordinate:
122nodeMarker  =mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= n1, coordinate = 1)) #y-coordinate
123
124#Spring-Damper between two marker coordinates
125mbs.AddObject(CoordinateSpringDamper(markerNumbers = [groundMarker, nodeMarker], 
126                                     stiffness = spring, damping = damper, 
127                                     springForceUserFunction = springForce,
128                                     visualization=VCoordinateSpringDamper(show=False))) 
129
130#add load:
131loadC = mbs.AddLoad(LoadCoordinate(markerNumber = nodeMarker, 
132                           load = load0))
133
134
135if not testIsActive:
136    sPos = mbs.AddSensor(SensorNode(nodeNumber=n1, storeInternal=True,#fileName="solution/sensorPos.txt"
137                             outputVariableType=exu.OutputVariableType.Position))
138    sVel = mbs.AddSensor(SensorNode(nodeNumber=n1, storeInternal=True,#fileName="solution/sensorVel.txt"
139                                    outputVariableType=exu.OutputVariableType.Velocity))
140    sAcc = mbs.AddSensor(SensorNode(nodeNumber=n1, storeInternal=True,#fileName="solution/sensorAcc.txt"
141                                    outputVariableType=exu.OutputVariableType.Acceleration))
142    #dummy, for PlotSensor
143    #these files are created, if doRefSol=True:
144    sPosRef = mbs.AddSensor(SensorNode(nodeNumber=n1, outputVariableType=exu.OutputVariableType.Position, 
145                                       storeInternal=not doRefSol,fileName="solution/sensorPosRef.txt",
146                                       writeToFile=doRefSol)) #set True to compute reference solution
147    sVelRef = mbs.AddSensor(SensorNode(nodeNumber=n1, outputVariableType=exu.OutputVariableType.Velocity, 
148                                       storeInternal=not doRefSol,fileName="solution/sensorVelRef.txt",
149                                       writeToFile=doRefSol)) #set True to compute reference solution
150    sAccRef = mbs.AddSensor(SensorNode(nodeNumber=n1, outputVariableType=exu.OutputVariableType.Acceleration, 
151                                       storeInternal=not doRefSol,fileName="solution/sensorAccRef.txt",
152                                       writeToFile=doRefSol)) #set True to compute reference solution
153
154#exu.Print(mbs)
155mbs.Assemble()
156
157simulationSettings = exu.SimulationSettings()
158simulationSettings.timeIntegration.newton.useModifiedNewton = False #Just for the test; modified Newton is usually faster
159simulationSettings.solution.file.write = False
160simulationSettings.solution.sensors.writePeriod = 1e-5
161simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
162simulationSettings.timeIntegration.endTime = tEnd
163simulationSettings.timeIntegration.minimumStepSize = 1e-10
164simulationSettings.timeIntegration.stepInformation = 3 #do not show step increase
165
166#important settings for contact:
167simulationSettings.timeIntegration.discontinuous.iterationTolerance = 1e-8 #this is the accepted penetration before reducing step size
168simulationSettings.timeIntegration.discontinuous.maxIterations = 1 #immediately perform step reduction
169simulationSettings.timeIntegration.discontinuous.ignoreMaxIterations = False #repeat step in case of failure
170simulationSettings.timeIntegration.adaptiveStepRecoverySteps = 0 #number of steps to wait until step size is increased again
171simulationSettings.timeIntegration.adaptiveStepIncrease = 10    #after successful step, increase again rapidly
172
173
174simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1
175
176simulationSettings.show.statistics = True
177simulationSettings.timeIntegration.verboseMode = 1
178#simulationSettings.timeIntegration.realtime.active = True
179
180if not testIsActive:
181    SC.renderer.Start()              #start graphics visualization
182    #SC.renderer.DoIdleTasks()    #wait for pressing SPACE bar to continue
183
184#start solver:
185mbs.SolveDynamic(solverType=exu.DynamicSolverType.TrapezoidalIndex2, simulationSettings=simulationSettings)
186#mbs.SolveDynamic(solverType=exu.DynamicSolverType.RK67, simulationSettings=simulationSettings)
187
188if not testIsActive:
189    #SC.renderer.DoIdleTasks()#wait for pressing 'Q' to quit
190    SC.renderer.Stop()               #safely close rendering window!
191
192u = mbs.GetNodeOutput(n1, exu.OutputVariableType.Position)
193exu.Print('postNewtonStepContactTest=',u[1])
194
195exu.sys['testResult'] = u[1]
196
197
198if not testIsActive:
199    
200    import matplotlib.pyplot as plt
201    plt.close('all')
202
203    mbs.PlotSensor(sensorNumbers=[sPos, sPosRef], components=[1,1], figureName='Pos')
204    mbs.PlotSensor(sensorNumbers=[sVel, sVelRef], components=[1,1], figureName='Vel')
205    mbs.PlotSensor(sensorNumbers=[sAcc, sAccRef], components=[1,1], figureName='Acc')