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#+++++++++++++++++++++++++++++++++++++++++++++++++++++