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')