contactCurvePolynomial.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Test for ObjectContactCurveCircles with polynomial enhancement; pin moving on smooth circular arc
  5#
  6# Author:   Johannes Gerstmayr
  7# Date:     2025-05-12
  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#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 12import exudyn as exu
 13from exudyn.utilities import *
 14import exudyn.graphics as graphics
 15import numpy as np
 16
 17
 18useGraphics = True #without test
 19SC = exu.SystemContainer()
 20mbs = SC.AddSystem()
 21
 22L = 2
 23a = 0.1
 24b = 0.1
 25rPin = a*0.75
 26
 27oGround = mbs.CreateGround()
 28
 29contactStiffness = 1e6
 30contactDamping = 1e3
 31
 32inertiaBrick=InertiaCuboid(1000, sideLengths=[a,a,b])
 33
 34#sprocket wheel:
 35oPin = mbs.CreateRigidBody(
 36                            referencePosition=[0.99*L-rPin,0,0],
 37                            inertia=inertiaBrick,
 38                            gravity = [0,-g,0],
 39                            create2D=True,
 40                            graphicsDataList=[graphics.Brick(size=[a,a,b], color=graphics.color.dodgerblue),
 41                                              graphics.Sphere(point=[0*0.5*L,0,0], radius=rPin, nTiles=16)
 42                                              ]
 43                            )
 44
 45mBody = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oPin, localPosition=[0*0.5*L,0,0]))
 46sVel = mbs.AddSensor(SensorBody(bodyNumber = oPin, storeInternal=True,
 47                                outputVariableType=exu.OutputVariableType.Velocity))
 48sAcc = mbs.AddSensor(SensorBody(bodyNumber = oPin, storeInternal=True,
 49                                outputVariableType=exu.OutputVariableType.Acceleration))
 50
 51pList = []
 52pList3D = []
 53zPos = 0
 54n = 128 #8 points is already quite smooth
 55pData = np.zeros((n,2))
 56
 57deltaPhi = 1/(n-1)*pi
 58r = L
 59for i in range(1*n):
 60    phi = -0.5*pi+i*deltaPhi
 61    x = r*sin(phi)
 62    y = -r*cos(phi)
 63    pList.append([x,y])
 64    #polynomial enhancement with 2 cubic polys; xi=[0,c]:
 65    c = 2*r*sin(deltaPhi/2)
 66    theta = 2*np.arcsin(c/(2*r))
 67    slopeAng = theta/2
 68    slope = np.tan(slopeAng) #this is the slope we like to add as poly coeff
 69    #print('theta=',theta*180/pi, ', slopeAng=',slopeAng*180/pi)
 70    pData[i,:] = [-slope,slope]
 71
 72
 73nPoints = len(pList)
 74
 75segmentsData = np.zeros((nPoints,4))
 76segmentsData[:,0:2] = pList
 77segmentsData[:,2:4] = np.roll(pList,2) #roll is element wise on rows and columns, therefore 2>shift one row
 78
 79segmentsData = segmentsData[1:,:]
 80nSeg = len(segmentsData)
 81
 82polynomialData = exu.MatrixContainer(pData)
 83
 84mbs.CreateGround(graphicsDataList=[graphics.CheckerBoard(point=[0,0,-b], size=3*L)])
 85
 86mGround = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround, localPosition=[0,0,0]))
 87
 88nGenericData = mbs.AddNode(NodeGenericData(initialCoordinates=[-1,0,0]*nSeg,
 89                                           numberOfDataCoordinates=3*nSeg))
 90
 91mbs.AddObject(ObjectContactCurveCircles(markerNumbers=[mGround]+[mBody],
 92                                        nodeNumber=nGenericData,
 93                                        circlesRadii=[rPin],
 94                                        segmentsData=exu.MatrixContainer(segmentsData),
 95                                        #polynomialData=polynomialData,
 96                                        contactStiffness=contactStiffness, contactDamping=contactDamping,
 97                                        visualization=VObjectContactCurveCircles(show=True, color=graphics.color.blue)
 98                                        ))
 99
100
101
102mbs.Assemble()
103
104stepSize=0.0001
105tEnd = 3
106simulationSettings = exu.SimulationSettings()
107simulationSettings.solutionSettings.writeSolutionToFile = True
108simulationSettings.solutionSettings.solutionWritePeriod = 0.001
109simulationSettings.solutionSettings.sensorsWritePeriod = stepSize  #output interval
110simulationSettings.timeIntegration.numberOfSteps = int(tEnd/stepSize)
111simulationSettings.timeIntegration.endTime = tEnd
112simulationSettings.timeIntegration.simulateInRealtime = True
113#simulationSettings.timeIntegration.realtimeFactor = 0.2
114# simulationSettings.timeIntegration.discontinuous.iterationTolerance = 1e-2
115# simulationSettings.timeIntegration.discontinuous.useRecommendedStepSize = False
116
117#simulationSettings.linearSolverType = exu.LinearSolverType.EigenSparse
118simulationSettings.timeIntegration.newton.useModifiedNewton = True
119#simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1
120
121simulationSettings.timeIntegration.verboseMode = 1
122
123SC.visualizationSettings.general.graphicsUpdateInterval = 0.02
124SC.visualizationSettings.window.renderWindowSize=[1600,2000]
125SC.visualizationSettings.openGL.multiSampling=4
126#SC.visualizationSettings.openGL.facesTransparent=True
127SC.visualizationSettings.openGL.shadow=0.3
128SC.visualizationSettings.loads.show = False
129SC.visualizationSettings.connectors.showContact = True
130SC.visualizationSettings.contact.tilingCurves = 16
131
132SC.renderer.Start()              #start graphics visualization
133SC.renderer.DoIdleTasks()    #wait for pressing SPACE bar to continue
134
135#start solver:
136mbs.SolveDynamic(simulationSettings)
137
138# SC.renderer.DoIdleTasks()#wait for pressing 'Q' to quit
139SC.renderer.Stop()               #safely close rendering window!
140
141if True:
142    mbs.PlotSensor([sVel, sAcc], components=[1,1])
143
144
145#mbs.SolutionViewer()