CoordinateSpringDamperTest.py

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

 1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 2# This is an EXUDYN example
 3#
 4# Details:  The same mass point as `CartesianSpringDamperTest`, held by three `CoordinateSpringDamper`
 5#           connectors instead - one per coordinate, the first one loaded. The two results differ in
 6#           the last digits because this one carries an offset.
 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
27L=0.5
28mass = 1.6
29k = 4000
30omega0 = 50 # sqrt(4000/1.6)
31dRel = 0.05
32d = dRel * 2 * 80 #80=sqrt(1.6*4000)
33u0=-0.08
34v0=1
35f = 80
36x0 = f/k
37fFriction = 20 #force in Newton, only depends on direction of velocity
38
39#node for mass point:
40n1=mbs.AddNode(Point(referenceCoordinates = [L,0,0], initialCoordinates = [u0,0,0], initialVelocities= [v0,0,0]))
41nGround=mbs.AddNode(NodePointGround(referenceCoordinates = [L,0,0]))
42
43#add mass points and ground object:
44objectGround = mbs.AddObject(ObjectGround(referencePosition = [0,0,0]))
45massPoint = mbs.AddObject(MassPoint(mass = mass, nodeNumber = n1))
46
47#marker for constraint / springDamper
48groundCoordinateMarker = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= nGround, coordinate = 0))
49nodeCoordinateMarker0  = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= n1, coordinate = 0))
50nodeCoordinateMarker1  = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= n1, coordinate = 1))
51nodeCoordinateMarker2  = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= n1, coordinate = 2))
52
53#Spring-Dampers
54mbs.AddObject(CoordinateSpringDamper(markerNumbers = [groundCoordinateMarker, nodeCoordinateMarker0], 
55                                     stiffness = k, damping = d) ) #changed 2023-01-21: dryFriction=0*fFriction, dryFrictionProportionalZone=0.01)) #offset must be zero, because coordinates just represent the displacements
56mbs.AddObject(CoordinateSpringDamper(markerNumbers = [groundCoordinateMarker, nodeCoordinateMarker1], stiffness = k)) 
57mbs.AddObject(CoordinateSpringDamper(markerNumbers = [groundCoordinateMarker, nodeCoordinateMarker2], stiffness = k)) 
58
59#add loads:
60mbs.AddLoad(LoadCoordinate(markerNumber = nodeCoordinateMarker0, load = f))
61
62mbs.Assemble()
63
64simulationSettings = exu.SimulationSettings()
65tEnd = 1 #1
66steps = 1000    #1000
67simulationSettings.solution.file.writePeriod = 1e-3
68simulationSettings.timeIntegration.numberOfSteps = steps
69simulationSettings.timeIntegration.endTime = tEnd
70simulationSettings.show.statistics = False
71
72simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1 #SHOULD work with 0.9 as well
73simulationSettings.solution.file.write=False
74
75if not testIsActive: 
76    SC.renderer.Start()
77
78mbs.SolveDynamic(simulationSettings)
79
80if not testIsActive: 
81    SC.renderer.DoIdleTasks()
82    SC.renderer.Stop() #safely close rendering window!
83
84u = mbs.GetNodeOutput(n1, exu.OutputVariableType.Position)
85uCoordinateSpringDamper= u[0] - L
86if True:
87    errorCoordinateSpringDamper = uCoordinateSpringDamper - 0.011834933407061654 #until 2022-01-25 (changed jacobians): 0.01183493340619235 #2021-09-27: new JacobianODE2RHS
88    #errorCoordinateSpringDamper = uCoordinateSpringDamper -0.011834933407368853 #2021-02-04: 0.011834933407368853
89else:
90    errorCoordinateSpringDamper = uCoordinateSpringDamper - 0.011834933406690284 #15.12.2019: 0.011834933406690284; beofre 15.12.2019: 0.011834933407047 #for 1000 steps, endtime=1; this is different from CartesianSpringDamper because of offset L (rounding errors around 1e-14)
91
92
93exu.Print('solution CoordinateSpringDamper=',uCoordinateSpringDamper)
94testResult = abs(errorCoordinateSpringDamper)
95exu.Print('solution of CoordinateSpringDamperTest=', testResult)
96exu.sys['testResult'] = testResult