CartesianSpringDamperTest.py

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

 1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 2# This is an EXUDYN example
 3#
 4# Details:  A mass point held by a `CartesianSpringDamper` against a constant force, integrated over
 5#           1 s in 1000 steps: the x-displacement at the end, which agrees with the analytical
 6#           solution to about 3e-6.
 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
37
38#node for mass point:
39n1=mbs.AddNode(Point(referenceCoordinates = [L,0,0], initialCoordinates = [u0,0,0], initialVelocities= [v0,0,0]))
40
41#add mass points and ground object:
42objectGround = mbs.AddObject(ObjectGround(referencePosition = [0,0,0]))
43massPoint = mbs.AddObject(MassPoint(mass = mass, nodeNumber = n1))
44
45#marker for constraint / springDamper
46groundMarker = mbs.AddMarker(MarkerBodyPosition(bodyNumber = objectGround, localPosition= [0, 0, 0]))
47bodyMarker = mbs.AddMarker(MarkerBodyPosition(bodyNumber = massPoint, localPosition= [0, 0, 0]))
48
49mbs.AddObject(CartesianSpringDamper(markerNumbers = [groundMarker, bodyMarker], stiffness = [k,k,k], damping = [d,0,0], offset = [L,0,0]))
50
51#add loads:
52mbs.AddLoad({'loadType': 'ForceVector',  'markerNumber': bodyMarker,  'loadVector': [f, 0, 0]}) 
53
54#exu.Print(mbs)
55mbs.Assemble()
56
57simulationSettings = exu.SimulationSettings()
58simulationSettings.solution.file.write=False
59tEnd = 1
60steps = 1000
61simulationSettings.timeIntegration.numberOfSteps = steps
62simulationSettings.timeIntegration.endTime = tEnd
63simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1 #SHOULD work with 0.9 as well
64simulationSettings.show.statistics = False
65
66mbs.SolveDynamic(simulationSettings)
67
68u = mbs.GetNodeOutput(n1, exu.OutputVariableType.Position)
69uCartesianSpringDamper= u[0] - L
70if True:
71    errorCartesianSpringDamper = uCartesianSpringDamper - 0.011834933407061654 #until 2022-01-25 (changed jacobians): 0.011834933406052683 #2021-09-27: new JacobianODE2RHS
72    #errorCartesianSpringDamper = uCartesianSpringDamper -0.011834933407364412 #2021-02-04: 0.011834933407364412 
73else:
74    errorCartesianSpringDamper = uCartesianSpringDamper - 0.011834933407594783 #15.12.2019: 0.011834933407594783; beofre 15.12.2019: 0.011834933407038783 #for 1000 steps, endtime=1; accurate up to 3e-6 to exact solution
75exu.Print('solution cartesianSpringDamper=',uCartesianSpringDamper)
76
77testResult = abs(errorCartesianSpringDamper)
78exu.Print('solution of CartesianSpringDamperTest=', testResult)
79exu.sys['testResult'] = testResult