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