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