contactSphereSphereTest.py

You can view and download this file on Github: contactSphereSphereTest.py

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Test with user-defined load function and user-defined spring-damper function (Duffing oscillator)
  5#
  6# Author:   Johannes Gerstmayr
  7# Date:     2024-11-01
  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
 30SC = exu.SystemContainer()
 31mbs = SC.AddSystem()
 32exu.Print('EXUDYN version='+exu.GetVersionString())
 33
 34use2contacts = True
 35yInit = 0.5
 36vyInit = 10
 37radius=0.1
 38mass = 1.6                              #mass in kg
 39contactStiffness = 1e6                  #stiffness of spring-damper in N/m
 40contactDamping = 0*5e-4*contactStiffness  #damping constant in N/(m/s)
 41restitutionCoefficient = 0.7
 42impactModel = 2
 43
 44tEnd = 2     #end time of simulation
 45stepSize = 2e-5*20
 46g = 50*0
 47exu.Print('impact vel=', np.sqrt(2*yInit*g))
 48
 49mbs.CreateGround(referencePosition=[0,-2*radius,0],
 50                graphicsDataList=[graphics.Sphere(radius=radius, color=graphics.color.green, nTiles=64)])
 51#ground node
 52nGround1=mbs.AddNode(NodePointGround(referenceCoordinates = [0,-2*radius,0]))
 53nGround2=mbs.AddNode(NodePointGround(referenceCoordinates = [0,1+2*radius,0]))
 54
 55#add mass point (this is a 3D object with 3 coordinates):
 56massPoint = mbs.CreateRigidBody(referencePosition=[0,yInit,0],
 57                                initialVelocity=[0,vyInit,0],
 58                                inertia=InertiaSphere(mass, radius),
 59                                gravity = [0,-g,0],
 60                                graphicsDataList=[graphics.Sphere(radius=radius, color=graphics.color.orange, nTiles=64)])
 61nMassPoint = mbs.GetObject(massPoint)['nodeNumber']
 62
 63mGround1 = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nGround1))
 64mMass = mbs.AddMarker(MarkerBodyRigid(bodyNumber=massPoint))
 65
 66nData1 = mbs.AddNode(NodeGenericData(initialCoordinates=[0.1,0,0,0],
 67                                    numberOfDataCoordinates=4))
 68oSSC = mbs.AddObject(ObjectContactSphereSphere(markerNumbers=[mGround1, mMass],
 69                                               nodeNumber=nData1,
 70                                               spheresRadii=[radius,radius],
 71                                               contactStiffness = contactStiffness,
 72                                               # contactStiffnessExponent=1.5,
 73                                               contactDamping = contactDamping,
 74                                               impactModel = impactModel,
 75                                               restitutionCoefficient = restitutionCoefficient,
 76                                               ))
 77if use2contacts:
 78    mbs.CreateGround(referencePosition=[0,1+2*radius,0],
 79                    graphicsDataList=[graphics.Sphere(radius=radius, color=graphics.color.green, nTiles=64)])
 80    mGround2 = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nGround2))
 81    nData2 = mbs.AddNode(NodeGenericData(initialCoordinates=[0.1,0,0,0],
 82                                        numberOfDataCoordinates=4))
 83    oSSC = mbs.AddObject(ObjectContactSphereSphere(markerNumbers=[mGround2, mMass],
 84                                                    nodeNumber=nData2,
 85                                                    spheresRadii=[radius,radius],
 86                                                    contactStiffness = contactStiffness,
 87                                                    # contactStiffnessExponent=1.5,
 88                                                    contactDamping = contactDamping,
 89                                                    impactModel = impactModel,
 90                                                    restitutionCoefficient = restitutionCoefficient,
 91                                                    ))
 92
 93
 94sPos=mbs.AddSensor(SensorBody(bodyNumber=massPoint, storeInternal=True,
 95                                outputVariableType=exu.OutputVariableType.Position))
 96sVel=mbs.AddSensor(SensorBody(bodyNumber=massPoint, storeInternal=True,
 97                                outputVariableType=exu.OutputVariableType.Velocity))
 98
 99#exu.Print(mbs)
100mbs.Assemble()
101
102simulationSettings = exu.SimulationSettings()
103simulationSettings.solutionSettings.writeSolutionToFile = False
104simulationSettings.solutionSettings.solutionWritePeriod = 0.02
105simulationSettings.solutionSettings.sensorsWritePeriod = stepSize  #output interval
106simulationSettings.timeIntegration.numberOfSteps = int(tEnd/stepSize)
107simulationSettings.timeIntegration.endTime = tEnd
108#simulationSettings.timeIntegration.discontinuous.iterationTolerance = 1e-8
109# simulationSettings.timeIntegration.discontinuous.useRecommendedStepSize = False
110
111simulationSettings.timeIntegration.newton.useModifiedNewton = True
112simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1
113
114simulationSettings.displayStatistics = True
115simulationSettings.timeIntegration.verboseMode = 1
116
117SC.visualizationSettings.window.renderWindowSize=[1600,2000]
118SC.visualizationSettings.openGL.multiSampling=4
119
120if useGraphics:
121    exu.StartRenderer()              #start graphics visualization
122    mbs.WaitForUserToContinue()    #wait for pressing SPACE bar to continue
123
124#start solver:q
125mbs.SolveDynamic(simulationSettings)
126
127if useGraphics:
128    SC.WaitForRenderEngineStopFlag()#wait for pressing 'Q' to quit
129    exu.StopRenderer()               #safely close rendering window!
130
131#evaluate final (=current) output values
132# u = mbs.GetNodeOutput(n1, exu.OutputVariableType.Position)
133# exu.Print('u     =',u)
134uTotal = mbs.GetNodeOutput(nMassPoint, exu.OutputVariableType.CoordinatesTotal)
135exu.Print('uTotal=',uTotal[1])
136
137# mbs.SolutionViewer()
138
139#plot results:
140if useGraphics:
141    mbs.PlotSensor([sPos,sVel], components=[1,1], closeAll=True)
142
143#+++++++++++++++++++++++++++++++++++++++++++++++++++++
144
145exudynTestGlobals.testError = uTotal[1] - (0.7092489359461815)
146exudynTestGlobals.testResult = uTotal[1]
147
148
149#+++++++++++++++++++++++++++++++++++++++++++++++++++++
150def MaximaAfterCrossings(y):
151    # Identify where the signal crosses from (-) to (+)
152    crossings = np.where((y[:-1] < 0) & (y[1:] > 0))[0] + 1  # Indices where crossing occurs
153
154    # Add the end of the signal as the last segment
155    crossings = np.concatenate(([0], crossings, [len(y)]))
156
157    # Calculate the maximum of each segment
158    maxima = [np.max(y[crossings[i]:crossings[i + 1]]) for i in range(len(crossings) - 1)]
159    minima = [np.min(y[crossings[i]:crossings[i + 1]]) for i in range(len(crossings) - 1)]
160
161    return [maxima[0],-minima[0],maxima[1],-minima[1],maxima[2],-minima[2]]
162    # return [maxima[0],-minima[0]]
163
164# Example usage:
165y = mbs.GetSensorStoredData(sVel)[:,2]
166maxima = MaximaAfterCrossings(y) #max 4 crossings
167exu.Print("Maxima after each crossing:", maxima)
168#exu.Print('relations=',maxima[1]/maxima[0],maxima[2]/maxima[1],maxima[3]/maxima[2])
169exu.Print('relations=',maxima[1]/maxima[0],maxima[2]/maxima[1],maxima[3]/maxima[2],maxima[4]/maxima[3])
170
171#+++++++++++++++++++++++++++++++++++++++++++++++++++++
172#Results:
173#e=0.5, stepSize=1e-5, Gonthier-CarvalhoMartins:
174#Maxima after each crossing: [10.0, 5.001784818335925, 2.501696294255996, 1.251250609056101]
175#e=0.5, stepSize=2e-5, Gonthier-CarvalhoMartins:
176#Maxima after each crossing: [10.0, 5.001975126052718, 2.501837721858895, 1.2513458250009686]
177#e=0.5, stepSize=5e-5, Gonthier-CarvalhoMartins:
178#Maxima after each crossing: [10.0, 5.003239017129689, 2.501566602043766, 1.2509080365758343]
179
180#e=0.9, stepSize=2e-5, Gonthier-CarvalhoMartins:
181#Maxima after each crossing: [10.0, 9.000168963897542, 8.100761764779794, 7.290837898880879]
182#e=0.8, stepSize=2e-5, Gonthier-CarvalhoMartins:
183#Maxima after each crossing: [10.0, 8.000159389915021, 6.400886105172352, 5.120771413440074]
184#e=0.3, stepSize=2e-5, Gonthier-CarvalhoMartins:
185#Maxima after each crossing: [10.0, 3.0423710530759065, 0.9256019663437189]
186#e=0.1, stepSize=2e-5, Gonthier-CarvalhoMartins:
187#Maxima after each crossing: [10.0, 1.0098986986792768, 0.10198953720292783]
188#e=0.025, stepSize=2e-5, Gonthier-CarvalhoMartins:
189#Maxima after each crossing: [10.0, 0.25015634771732226]
190
191#e=0.9, stepSize=2e-5, Hunt-Crossley:
192#Maxima after each crossing: [10.0, 9.090298887177969, 8.263743088547175, 7.512142692359256]
193#e=0.8, stepSize=2e-5, Hunt-Crossley:
194#Maxima after each crossing: [10.0, 8.32903493005253, 6.937608834244499, 5.778355621679387]
195#e=0.5, stepSize=2e-5, Hunt-Crossley:
196#Maxima after each crossing: [10.0, 6.629838391096419, 4.3951552519078065, 2.913688143667944]
197#e=0.3, stepSize=2e-5, Hunt-Crossley:
198#Maxima after each crossing: [10.0, 5.813176318310211, 3.3794625180566293, 1.964542613454227]
199#e=0.1, stepSize=2e-5, Hunt-Crossley:
200#Maxima after each crossing: [10.0, 5.158573889338283, 2.6613482605395395, 1.3728789060291537]
201#+++++++++++++++++++++++++++++++++++++++++++++++++++++