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