contactSphereSphereTestEAPM.py
You can view and download this file on Github: contactSphereSphereTestEAPM.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: Test of the SphereSphere contact module
5#
6# Author: Johannes Gerstmayr and Sebastian Weyrer
7# Date: 2025-01-07
8#
9# 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.
10#
11#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
12
13import exudyn as exu
14from exudyn.utilities import *
15import exudyn.graphics as graphics
16import numpy as np
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
30
31SC = exu.SystemContainer()
32mbs = SC.AddSystem()
33exu.Print('EXUDYN version = ' + exu.GetVersionString())
34
35#+++++++++++++++++++++++++++++++++++++++++++++
36# define behaviour of script
37T = 1 # [s] simulation time
38stepSize = 2e-5 # [s]
39impactModel = 0 # 0 = Adhesive Elasto-Plastic Model; 1 = Hunt-Crossley Model; 2 = mixed Gonthier/EtAl-Carvalho/Martins Model
40restitutionCoefficient = 0.9 # used restitution coefficient if impactModel != 0
41use2Contacts = True # only applies if impactModel != 0, specifies numbers of used contacts for chosen impact model
42
43#+++++++++++++++++++++++++++++++++++++++++++++
44# define parameters of the speheres
45radius = 0.1 # [m] radius of spheres
46mass = 1.6 # [kg] mass of each sphere
47# define reference position of the moving sphere and initial velocity
48if impactModel == 0:
49 xInit = 0
50 yInit = 2 * radius # gap = 0 at beginning of simulation
51 vyInit = 0
52else:
53 xInit = 0
54 yInit = 2 * radius + 0.5 # gap = 0.5m at beginning of simulation
55 vyInit = 10
56
57
58#+++++++++++++++++++++++++++++++++++++++++++++
59# create ground (with marker) that is needed for every test case
60mbs.CreateGround(referencePosition=[0, 0, 0], graphicsDataList=[graphics.Sphere(radius=radius, color=graphics.color.red, nTiles=64)])
61nGround = mbs.AddNode(NodePointGround(referenceCoordinates = [0, 0, 0]))
62mGround = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nGround))
63
64#+++++++++++++++++++++++++++++++++++++++++++++
65# add mass point (with marker) for moving spehere
66massPoint = mbs.CreateRigidBody(referencePosition=[xInit, yInit, 0],
67 initialVelocity=[0, vyInit, 0],
68 inertia=InertiaSphere(mass, radius),
69 gravity=[0, 0, 0],
70 graphicsDataList=[graphics.Sphere(radius=radius, color=graphics.color.green, nTiles=64)])
71mMass = mbs.AddMarker(MarkerBodyRigid(bodyNumber=massPoint))
72nMassPoint = mbs.GetObject(massPoint)['nodeNumber']
73
74#+++++++++++++++++++++++++++++++++++++++++++++
75# if impactModel = 0, simulate loading and unloading of (soil) particles
76if impactModel == 0:
77 tMove = T/2
78 dMove = 0.002
79 def UFforce(mbs, t, itemNumber, u, v, k, d, F0):
80 if t < tMove:
81 setPosition = 2*radius - (dMove/tMove)*t
82 else:
83 setPosition = 2*radius - dMove + (dMove/tMove)*(t-tMove)
84 # print(setPosition)
85 return (u-setPosition)*1e5
86
87 mbs.AddObject(SpringDamper(markerNumbers=[mGround, mMass], referenceLength=0, springForceUserFunction=UFforce))
88 nData = mbs.AddNode(NodeGenericData(initialCoordinates=[0, 0, 0, 0], numberOfDataCoordinates=4))
89 oSSC = mbs.AddObject(ObjectContactSphereSphere(markerNumbers=[mGround, mMass], # m1 is moving spehere
90 nodeNumber=nData,
91 spheresRadii=[radius,radius],
92 impactModel=impactModel,
93 contactStiffness=1e5,
94 contactStiffnessExponent=2,
95 contactDamping=0.001,
96 contactPlasticityRatio=0.6,
97 constantPullOffForce=0.01,
98 adhesionCoefficient=4e4,
99 adhesionExponent=2))
100
101#+++++++++++++++++++++++++++++++++++++++++++++
102# if impactModel != 0 do impact simulation of particles
103else:
104 nData = mbs.AddNode(NodeGenericData(initialCoordinates=[0, 0, 0, 0], numberOfDataCoordinates=4))
105 oSSC = mbs.AddObject(ObjectContactSphereSphere(markerNumbers=[mGround, mMass], # m1 is moving spehere
106 nodeNumber=nData,
107 spheresRadii=[radius,radius],
108 contactStiffness=1e6,
109 contactDamping=0,
110 impactModel=impactModel,
111 restitutionCoefficient=restitutionCoefficient,
112 minimumImpactVelocity=0))
113 if use2Contacts:
114 # add this ground object with a gap of 0.5m above the moving spehere
115 mbs.CreateGround(referencePosition=[0, 1+4*radius, 0], graphicsDataList=[graphics.Sphere(radius=radius, color=graphics.color.red, nTiles=64)])
116 nGround1 = mbs.AddNode(NodePointGround(referenceCoordinates = [0, 1+4*radius, 0]))
117 mGround1 = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nGround1))
118 nData1 = mbs.AddNode(NodeGenericData(initialCoordinates=[0, 0, 0, 0], numberOfDataCoordinates=4))
119 oSSC1 = mbs.AddObject(ObjectContactSphereSphere(markerNumbers=[mGround1, mMass], # m1 is moving spehere
120 nodeNumber=nData1,
121 spheresRadii=[radius,radius],
122 contactStiffness=1e6,
123 contactDamping=0,
124 impactModel=impactModel,
125 restitutionCoefficient=restitutionCoefficient,
126 minimumImpactVelocity=0))
127
128#+++++++++++++++++++++++++++++++++++++++++++++
129# meassure output variables of moving spehere in any test case
130sPos = mbs.AddSensor(SensorBody(bodyNumber=massPoint, outputVariableType=exu.OutputVariableType.Position, writeToFile=False, storeInternal=True))
131sVel = mbs.AddSensor(SensorBody(bodyNumber=massPoint, outputVariableType=exu.OutputVariableType.Velocity, writeToFile=False, storeInternal=True))
132sGap = mbs.AddSensor(SensorObject(objectNumber=oSSC, outputVariableType=exu.OutputVariableType.DisplacementLocal, writeToFile=False, storeInternal=True))
133sContactForce = mbs.AddSensor(SensorObject(objectNumber=oSSC, outputVariableType=exu.OutputVariableType.Force, writeToFile=False, storeInternal=True))
134
135#+++++++++++++++++++++++++++++++++++++++++++++
136mbs.Assemble()
137# exu.Print(mbs)
138
139#+++++++++++++++++++++++++++++++++++++++++++++
140simulationSettings = exu.SimulationSettings()
141simulationSettings.timeIntegration.numberOfSteps = int(T/stepSize)
142simulationSettings.solutionSettings.sensorsWritePeriod = stepSize
143simulationSettings.timeIntegration.endTime = T
144simulationSettings.displayStatistics = True
145simulationSettings.timeIntegration.verboseMode = 1
146
147simulationSettings.timeIntegration.newton.useModifiedNewton = True
148simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1
149
150#+++++++++++++++++++++++++++++++++++++++++++++
151if useGraphics:
152 exu.StartRenderer() # start graphics visualization
153 mbs.WaitForUserToContinue() # wait for pressing SPACE bar to continue
154mbs.SolveDynamic(simulationSettings)
155if useGraphics:
156 SC.WaitForRenderEngineStopFlag() # wait for pressing 'Q' to quit
157 exu.StopRenderer() # safely close rendering window!
158
159
160#+++++++++++++++++++++++++++++++++++++++++++++
161uTotal = mbs.GetNodeOutput(nMassPoint, exu.OutputVariableType.CoordinatesTotal)
162exu.Print('uTotal=',uTotal[1])
163
164exudynTestGlobals.testResult = uTotal[1]
165
166if useGraphics:
167 #+++++++++++++++++++++++++++++++++++++++++++++
168 # plot sensor data for different test cases
169 import matplotlib.pyplot as plt
170 def setMatplotlibSettings():
171 import matplotlib as mpl
172 import matplotlib.font_manager as font_manager
173 import matplotlib.pyplot as plt
174 fsize = 9 # general fontsize
175 tsize = 9 # legend
176 # setting ticks, fontsize
177 tdir = 'out' # 'in' or 'out': direction of ticks
178 major = 5.0 # length major ticks
179 minor = 3.0 # length minor ticks
180 lwidth = 0.8 # width of the frame
181 lhandle = 2 # length legend handle
182 plt.style.use('default')
183 # set font to the locally saved computermodern type
184 plt.rcParams['font.family']='serif'
185 cmfont = font_manager.FontProperties(fname=mpl.get_data_path() + '/fonts/ttf/cmr10.ttf')
186 plt.rcParams['font.serif']=cmfont.get_name()
187 plt.rcParams['mathtext.fontset']='cm'
188 plt.rcParams['axes.unicode_minus']=False
189 plt.rcParams['font.size'] = fsize
190 plt.rcParams['legend.fontsize'] = tsize
191 plt.rcParams['xtick.direction'] = tdir
192 plt.rcParams['ytick.direction'] = tdir
193 plt.rcParams['xtick.major.size'] = major
194 plt.rcParams['xtick.minor.size'] = minor
195 plt.rcParams['ytick.major.size'] = major
196 plt.rcParams['ytick.minor.size'] = minor
197 plt.rcParams['axes.linewidth'] = lwidth
198 plt.rcParams['axes.formatter.use_mathtext'] = True
199 plt.rcParams['legend.handlelength'] = lhandle
200 plt.rcParams['lines.linewidth'] = 1.2
201 return
202 setMatplotlibSettings()
203 mbs.PlotSensor([sPos,sVel], components=[1, 1], closeAll=False)
204 plt.show(block=True) # for figures to stay open at end of plot routines
205 if impactModel == 0:
206 def add_arrow(line, position=None, direction='right', size=15, color=None):
207 if color is None:
208 color = line.get_color()
209 xdata = line.get_xdata()
210 ydata = line.get_ydata()
211 if position is None:
212 position = xdata.mean()
213 start_ind = np.argmin(np.absolute(xdata - position))
214 if direction == 'right':
215 end_ind = start_ind + 1
216 else:
217 end_ind = start_ind - 1
218 line.axes.annotate('', xytext=(xdata[start_ind], ydata[start_ind]),
219 xy=(xdata[end_ind], ydata[end_ind]),
220 arrowprops=dict(arrowstyle="->", color=color), size=size)
221 fig, (ax1, ax2) = plt.subplots(1, 2, figsize=[8, 4])
222 t = mbs.GetSensorStoredData(sGap)[:, 0]
223 penetration = -mbs.GetSensorStoredData(sGap)[:, 1]
224 globalContactForce = mbs.GetSensorStoredData(sContactForce)[:, 2]
225 repulsiveForce = -globalContactForce
226 indexFirstPositivePenetration = np.where(penetration > 0)[0][0]
227 indexLastPositivePenetration = np.where(penetration > 0)[0][-1]
228 line1 = ax1.plot(penetration[indexFirstPositivePenetration:indexLastPositivePenetration+1], repulsiveForce[indexFirstPositivePenetration:indexLastPositivePenetration+1])[0]
229 add_arrow(line1, position=0.001)
230 add_arrow(line1, position=0.0016)
231 ax1.grid()
232 ax1.set_xlabel('negative gap in m')
233 ax1.set_ylabel('force (m1) in N')
234 line2 = ax2.plot(penetration[indexFirstPositivePenetration:indexLastPositivePenetration+1], globalContactForce[indexFirstPositivePenetration:indexLastPositivePenetration+1])[0]
235 add_arrow(line2, position=0.001)
236 add_arrow(line2, position=0.0016)
237 ax2.grid()
238 ax2.set_xlabel('negative gap in m')
239 ax2.set_ylabel('force (m0) in N (= global force)')
240 plt.tight_layout()
241 if impactModel != 0:
242 def MaximaAfterCrossings(y):
243 # Identify where the signal crosses from (-) to (+)
244 crossings = np.where((y[:-1] < 0) & (y[1:] > 0))[0] + 1 # Indices where crossing occurs
245 # Add the end of the signal as the last segment
246 crossings = np.concatenate(([0], crossings, [len(y)]))
247 # Calculate the maximum of each segment
248 maxima = [np.max(y[crossings[i]:crossings[i + 1]]) for i in range(len(crossings) - 1)]
249 minima = [np.min(y[crossings[i]:crossings[i + 1]]) for i in range(len(crossings) - 1)]
250 return [maxima[0],-minima[0],maxima[1],-minima[1]]
251 # return [maxima[0],-minima[0]]
252 y = mbs.GetSensorStoredData(sVel)[:,2]
253 maxima = MaximaAfterCrossings(y) #max 4 crossings
254 exu.Print("Maxima after each crossing:", maxima)
255 exu.Print('relations=',maxima[1]/maxima[0],maxima[2]/maxima[1],maxima[3]/maxima[2])
256
257#+++++++++++++++++++++++++++++++++++++++++++++++++++++
258#Results:
259#e=0.5, stepSize=1e-5, Gonthier-CarvalhoMartins:
260#Maxima after each crossing: [10.0, 5.001784818335925, 2.501696294255996, 1.251250609056101]
261#e=0.5, stepSize=2e-5, Gonthier-CarvalhoMartins:
262#Maxima after each crossing: [10.0, 5.001975126052718, 2.501837721858895, 1.2513458250009686]
263#e=0.5, stepSize=5e-5, Gonthier-CarvalhoMartins:
264#Maxima after each crossing: [10.0, 5.003239017129689, 2.501566602043766, 1.2509080365758343]
265
266#e=0.9, stepSize=2e-5, Gonthier-CarvalhoMartins:
267#Maxima after each crossing: [10.0, 9.000168963897542, 8.100761764779794, 7.290837898880879]
268#e=0.8, stepSize=2e-5, Gonthier-CarvalhoMartins:
269#Maxima after each crossing: [10.0, 8.000159389915021, 6.400886105172352, 5.120771413440074]
270#e=0.3, stepSize=2e-5, Gonthier-CarvalhoMartins:
271#Maxima after each crossing: [10.0, 3.0423710530759065, 0.9256019663437189]
272#e=0.1, stepSize=2e-5, Gonthier-CarvalhoMartins:
273#Maxima after each crossing: [10.0, 1.0098986986792768, 0.10198953720292783]
274#e=0.025, stepSize=2e-5, Gonthier-CarvalhoMartins:
275#Maxima after each crossing: [10.0, 0.25015634771732226]
276
277#e=0.9, stepSize=2e-5, Hunt-Crossley:
278#Maxima after each crossing: [10.0, 9.090298887177969, 8.263743088547175, 7.512142692359256]
279#e=0.8, stepSize=2e-5, Hunt-Crossley:
280#Maxima after each crossing: [10.0, 8.32903493005253, 6.937608834244499, 5.778355621679387]
281#e=0.5, stepSize=2e-5, Hunt-Crossley:
282#Maxima after each crossing: [10.0, 6.629838391096419, 4.3951552519078065, 2.913688143667944]
283#e=0.3, stepSize=2e-5, Hunt-Crossley:
284#Maxima after each crossing: [10.0, 5.813176318310211, 3.3794625180566293, 1.964542613454227]
285#e=0.1, stepSize=2e-5, Hunt-Crossley:
286#Maxima after each crossing: [10.0, 5.158573889338283, 2.6613482605395395, 1.3728789060291537]
287#+++++++++++++++++++++++++++++++++++++++++++++++++++++