SlidingJoint2DTest.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  A rigid body sliding along a chain of `ANCFCable2D` elements with an
  5#           `ObjectJointSliding2D`, which moves the contact from one element to the next as it goes:
  6#           the vertical position of its centre of mass after a short dynamic step.
  7#           The model compares against a reference value written into it, so its
  8#           result is that difference and its reference solution is 0 (#2632).
  9#
 10# Author:   Johannes Gerstmayr
 11# Date:     2019-11-01, reworked 2026-09-24
 12#
 13# 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.
 14#
 15#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 16
 17import exudyn as exu
 18from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
 19import exudyn.graphics as graphics
 20
 21testIsActive = exu.sys.get('testIsActive', False)
 22exu.sys['testTolerance'] = 4e-13 #the tolerance RunAllModelUnitTests used for these ten
 23
 24SC = exu.SystemContainer()
 25mbs = SC.AddSystem()
 26
 27
 28#background
 29rect = [-2.5,-2,2.5,1] #xmin,ymin,xmax,ymax
 30background0 = {'type':'Line', 'color':[0.1,0.1,0.8,1], 'data':[rect[0],rect[1],0, rect[2],rect[1],0, rect[2],rect[3],0, rect[0],rect[3],0, rect[0],rect[1],0]} #background
 31background1 = {'type':'Line', 'color':[0.1,0.1,0.8,1], 'data':[0,-1,0, 2,-1,0]} #background
 32oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0], visualization=VObjectGround(graphicsData= [background0])))
 33
 34
 35#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 36#cable:
 37mypi = 3.141592653589793
 38
 39L=2                     # length of ANCF element in m
 40#L=mypi                 # length of ANCF element in m
 41E=2.07e11               # Young's modulus of ANCF element in N/m^2
 42rho=7800                # density of ANCF element in kg/m^3
 43b=0.001                 # width of rectangular ANCF element in m
 44h=0.001                 # height of rectangular ANCF element in m
 45A=b*h                   # cross sectional area of ANCF element in m^2
 46I=b*h**3/12             # second moment of area of ANCF element in m^4
 47f=3*E*I/L**2            # tip load applied to ANCF element in N
 48g=9.81
 49
 50exu.Print("load f="+str(f))
 51exu.Print("EI="+str(E*I))
 52
 53nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0])) #ground node for coordinate constraint
 54mGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nGround, coordinate=0)) #Ground node ==> no action
 55
 56cableList=[]        #for cable elements
 57nodeList=[]  #for nodes of cable
 58markerList=[]       #for nodes
 59nc0 = mbs.AddNode(Point2DS1(referenceCoordinates=[0,0,1,0]))
 60nodeList+=[nc0]
 61nElements = 3
 62lElem = L / nElements
 63for i in range(nElements):
 64    nLast = mbs.AddNode(Point2DS1(referenceCoordinates=[lElem*(i+1),0,1,0]))
 65    nodeList+=[nLast]
 66    elem=mbs.AddObject(Cable2D(length=lElem, massPerLength=rho*A, 
 67                               bendingStiffness=E*I, axialStiffness=E*A, 
 68                               nodeNumbers=[int(nc0)+i,int(nc0)+i+1]))
 69    cableList+=[elem]
 70    mBody = mbs.AddMarker(MarkerBodyMass(bodyNumber = elem))
 71    mbs.AddLoad(Gravity(markerNumber=mBody, loadVector=[0,-g,0]))
 72
 73mANCF0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nc0, coordinate=0))
 74mANCF1 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nc0, coordinate=1))
 75mANCF2 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nc0, coordinate=3))
 76
 77mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF0]))
 78mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF1]))
 79mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF2]))
 80
 81
 82a = 0.1     #y-dim/2 of gondula
 83b = 0.001    #x-dim/2 of gondula
 84massRigid = 12*0.01
 85inertiaRigid = massRigid/12*(2*a)**2
 86
 87slidingCoordinateInit = lElem*1.5 #0.75*L
 88initialLocalMarker = 1 #second element
 89if nElements<2:
 90    slidingCoordinateInit /= 3.
 91    initialLocalMarker = 0
 92
 93addRigidBody = True
 94nRigid = -1
 95if addRigidBody:
 96    #rigid body which slides:
 97    graphicsRigid = {'type':'Line', 'color':[0.1,0.1,0.8,1], 'data':[-b,-a,0, b,-a,0, b,a,0, -b,a,0, -b,-a,0]} #drawing of rigid body
 98    nRigid = mbs.AddNode(Rigid2D(referenceCoordinates=[slidingCoordinateInit,-a,0], initialVelocities=[0,0,0]));
 99    oRigid = mbs.AddObject(RigidBody2D(mass=massRigid, inertia=inertiaRigid,nodeNumber=nRigid,visualization=VObjectRigidBody2D(graphicsData= [graphicsRigid])))
100
101    markerRigidTop = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid, localPosition=[0.,a,0.])) #support point
102    mR2 = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid, localPosition=[ 0.,0.,0.])) #center of mass (for load)
103
104    mbs.AddLoad(Force(markerNumber = mR2, loadVector = [massRigid*g*0.1, -massRigid*g, 0]))
105
106
107#slidingJoint:
108addSlidingJoint = True
109if addSlidingJoint:
110    cableMarkerList = []#list of Cable2DCoordinates markers
111    offsetList = []     #list of offsets counted from first cable element; needed in sliding joint
112    offset = 0          #first cable element has offset 0
113    for item in cableList: #create markers for cable elements
114        m = mbs.AddMarker(MarkerBodyCable2DCoordinates(bodyNumber = item))
115        cableMarkerList += [m]
116        offsetList += [offset]
117        offset += lElem
118
119    #mGroundSJ = mbs.AddMarker(MarkerBodyPosition(bodyNumber = oGround, localPosition=[0.*lElem+0.75*L,0.,0.])) 
120    nodeDataSJ = mbs.AddNode(NodeGenericData(initialCoordinates=[initialLocalMarker,slidingCoordinateInit],numberOfDataCoordinates=2)) #initial index in cable list
121    slidingJoint = mbs.AddObject(ObjectJointSliding2D(name='slider', markerNumbers=[markerRigidTop,cableMarkerList[initialLocalMarker]], 
122                                                      slidingMarkerNumbers=cableMarkerList, slidingMarkerOffsets=offsetList, 
123                                                      nodeNumber=nodeDataSJ, useClassicalFormulation = False))
124
125
126mbs.Assemble()
127
128simulationSettings = exu.SimulationSettings() #takes currently set values or default values
129simulationSettings.solution.file.write=False
130
131fact = 200
132simulationSettings.timeIntegration.numberOfSteps = 1*fact
133simulationSettings.timeIntegration.endTime = 0.001*fact*0.5
134simulationSettings.solution.file.write = True
135simulationSettings.solution.file.writePeriod = simulationSettings.timeIntegration.endTime/fact
136simulationSettings.timeIntegration.verboseMode = 1
137
138simulationSettings.timeIntegration.newton.relativeTolerance = 1e-8*100 #10000
139simulationSettings.timeIntegration.newton.absoluteTolerance = 1e-10*100
140
141simulationSettings.timeIntegration.newton.useModifiedNewton = False
142simulationSettings.timeIntegration.newton.maxModifiedNewtonIterations = 5
143simulationSettings.timeIntegration.newton.numericalDifferentiation.addReferenceCoordinatesToEpsilon = False
144simulationSettings.timeIntegration.newton.numericalDifferentiation.minimumCoordinateSize = 1.e-3
145simulationSettings.timeIntegration.newton.numericalDifferentiation.relativeEpsilon = 1e-8 #6.055454452393343e-06*0.0001 #eps^(1/3)
146simulationSettings.timeIntegration.newton.modifiedNewtonContractivity = 1e8
147simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.6 #0.6 works well 
148simulationSettings.show.statistics = False
149
150if not testIsActive: 
151    SC.renderer.Start()
152
153mbs.SolveDynamic(simulationSettings)
154
155
156error = 0
157if nRigid != -1:
158    u = mbs.GetNodeOutput(nRigid, exu.OutputVariableType.Position) #tip node
159    if True: 
160        error = u[1] - (-0.14920183348514676 ) #2021-09-27: new JacobianODE2RHS
161        #error = u[1] -(-0.14920182499994944 ) #2021-02-06: -0.14920182499994944 (1e-9 different from old solver) 
162    elif True:
163        error = u[1] - (-0.14920182666080586) #2021-02-04: -0.14920182666080586
164    else:
165        error = u[1] - (-0.14920151345936586) #2019-12-26: -0.14920151345936586; 15.12.2019: (-0.1489879442762764); before 15.12.2019: (-0.14898795617249422) #2019-11-22; #20.10.2019: (-0.1489879501348149); 17.10.2019:-0.14898795468724652; old? :(-0.14898795002032308) #old error before projected sliding joint: (-0.14898792622401774) #y-position of COM of sliding body
166                    
167    exu.Print('value SlidingJoint2DTest=',u[1])
168    #exu.Print('error SlidingJoint2DTest=',error)
169
170if not testIsActive: 
171    SC.renderer.DoIdleTasks()
172    SC.renderer.Stop() #safely close rendering window!
173
174testResult = abs(error)
175exu.Print('solution of SlidingJoint2DTest=', testResult)
176exu.sys['testResult'] = testResult