SpringDamperMesh.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  A rectangular mesh of mass points connected by spring-dampers, hanging from its top row.
  5#           It is solved **statically and dynamically**, and the result is the sum of the two errors;
  6#           the static solver runs with numerical differentiation of the connectors.
  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
 27nBodies = 6
 28nBodies2 = 3
 29
 30for j in range(nBodies2): 
 31    body = mbs.AddObject({'objectType': 'Ground', 'referencePosition': [0,j,0]})
 32    mbs.AddMarker({'markerType': 'BodyPosition',  'bodyNumber': body,  'localPosition': [0.0, 0.0, 0.0], 'bodyFixed': False})
 33    for i in range(nBodies-1): 
 34        #does not work for analytical jacobian in static case: z-coordinate unconstrained
 35        # node = mbs.AddNode({'nodeType': 'Point','referenceCoordinates': [i+1, j, 0.0],
 36        #                     'initialCoordinates': [(i+1)*0.05*0, 0.0, 0.0], 
 37        #                     'initialVelocities': [0., 0., 0.],})
 38        # body = mbs.AddObject({'objectType': 'MassPoint', 'mass': 10, 'nodeNumber': node})
 39        node = mbs.AddNode(NodePoint2D(referenceCoordinates=[i+1, j],
 40                                       initialCoordinates=[(i+1)*0.05*0, 0]))
 41        body = mbs.AddObject(ObjectMassPoint2D(mass= 10, nodeNumber= node))
 42        mbs.AddMarker({'markerType': 'BodyPosition',  'bodyNumber': body,  'localPosition': [0.0, 0.0, 0.0], 'bodyFixed': False})
 43
 44#add spring-dampers:
 45for j in range(nBodies2-1): 
 46    for i in range(nBodies-1): 
 47        mbs.AddObject({'objectType': 'ConnectorSpringDamper', 'stiffness': 4000, 'damping': 10, 'force': 0,
 48                        'referenceLength':1, 'markerNumbers': [j*nBodies + i,j*nBodies + i+1]})
 49        mbs.AddObject({'objectType': 'ConnectorSpringDamper', 'stiffness': 4000, 'damping': 10, 'force': 0,
 50                        'referenceLength':1, 'markerNumbers': [j*nBodies + i,(j+1)*nBodies + i]})
 51        #diagonal spring: l*sqrt(2)
 52        mbs.AddObject({'objectType': 'ConnectorSpringDamper', 'stiffness': 4000, 'damping': 10, 'force': 0,
 53                        'referenceLength':sqrt2, 'markerNumbers': [j*nBodies + i,(j+1)*nBodies + i+1]})
 54
 55for i in range(nBodies-1): 
 56    j = nBodies2-1
 57    mbs.AddObject({'objectType': 'ConnectorSpringDamper', 'stiffness': 4000, 'damping': 10, 'force': 0,
 58                    'referenceLength':1, 'markerNumbers': [j*nBodies + i,j*nBodies + i+1]})
 59for j in range(nBodies2-1): 
 60    i = nBodies-1
 61    mbs.AddObject({'objectType': 'ConnectorSpringDamper', 'stiffness': 4000, 'damping': 10, 'force': 0,
 62                    'referenceLength':1, 'markerNumbers': [j*nBodies + i,(j+1)*nBodies + i]})
 63
 64#add loads:
 65mbs.AddLoad({'loadType': 'ForceVector',  'markerNumber': nBodies*nBodies2-1,  'loadVector': [0, -50*2, 0]})
 66
 67#add constraints for testing:
 68nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[-0.5,0,0])) #ground node for coordinate constraint
 69mGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nGround, coordinate=0)) #Ground node ==> no action
 70
 71mNC1 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = 1, coordinate=1))
 72mNC2 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = 2, coordinate=1))
 73
 74mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mNC1]))
 75
 76mbs.Assemble()
 77
 78if not testIsActive: 
 79    SC.renderer.Start()
 80
 81simulationSettings = exu.SimulationSettings()
 82simulationSettings.timeIntegration.numberOfSteps = 100
 83simulationSettings.timeIntegration.endTime = 1
 84simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.6
 85simulationSettings.timeIntegration.generalizedAlpha.useNewmark = False
 86simulationSettings.timeIntegration.generalizedAlpha.useIndex2Constraints = False
 87simulationSettings.timeIntegration.verboseMode = 0
 88simulationSettings.timeIntegration.newton.useModifiedNewton = True
 89simulationSettings.show.statistics = True
 90simulationSettings.solution.file.write=False
 91
 92SC.visualizationSettings.nodes.defaultSize = 0.05
 93
 94#    if not testIsActive: 
 95#        SC.renderer.DoIdleTasks()
 96
 97mbs.SolveDynamic(simulationSettings)
 98
 99if not testIsActive: 
100    SC.renderer.DoIdleTasks()
101
102u = mbs.GetNodeOutput(nBodies-2, exu.OutputVariableType.Position) #tip node
103exu.Print('dynamic tip displacement (y)=', u[1])
104if True:
105    dynamicError = u[1]-(-0.6385807469187298 )  #2021-02-04: -0.6385807469187298 
106else:
107    dynamicError = u[1]-(-0.6383785907891227)  #2019-12-26: -0.6383785907891227; 2019-12-15: (-0.6349442849103891); before 15.12.2019: (-0.6349442850473246)
108   
109#dynamic tip displacement for                                             s=1000: -0.6386766492418571,s=10000: -0.6386985511667669, s=100000: -0.6387006546281098
110#dynamic tip displacement for index2 Newmark: s=100: -0.6386060431598312, s=1000: -0.638699952155624, s=10000: -0.6387008780175608, s=100000: -0.6387008872717617
111
112
113#simulationSettings.solution.file.name = "staticSolution.txt"
114#simulationSettings.solution.file.append = False
115simulationSettings.staticSolver.newton.numericalDifferentiation.relativeEpsilon = 1e-5
116simulationSettings.staticSolver.newton.relativeTolerance = 1e-6
117simulationSettings.staticSolver.newton.absoluteTolerance = 1e-1
118simulationSettings.staticSolver.newton.numericalDifferentiation.forODE2Connectors = True #be compatible with old solution
119#simulationSettings.staticSolver.verboseMode = 1
120
121mbs.SolveStatic(simulationSettings)
122
123u = mbs.GetNodeOutput(nBodies-2, exu.OutputVariableType.Position) #tip node
124exu.Print('static tip displacement (y)=', u[1])
125staticError = u[1]-(-0.44056224799446486)
126
127totalError = abs(staticError) + abs(dynamicError)
128
129if not testIsActive: 
130    SC.renderer.DoIdleTasks()
131    SC.renderer.Stop() 
132
133testResult = totalError
134exu.Print('solution of SpringDamperMesh=', testResult)
135exu.sys['testResult'] = testResult