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