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