cartesianSpringDamper.py

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

 1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 2# This is an EXUDYN example
 3#
 4# Details:  Example with CartesianSpringDamper and reference solution, to create simple system with mass point and spring-damper
 5#
 6# Model:    Linear oscillator with mass point and spring-damper
 7#
 8# Author:   Johannes Gerstmayr
 9# Date:     2019-08-15
10#
11# 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.
12#
13# *clean example*
14#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
15
16## import libaries
17import exudyn as exu
18from exudyn.itemInterface import *
19from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
20import exudyn.graphics as graphics #only import if it does not conflict
21
22## set up system 'mbs'
23SC = exu.SystemContainer()
24mbs = SC.AddSystem()
25
26## define overall parameters for linear spring-damper
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
36#static result = f/k = 0.01;
37x0 = f/k
38
39## create ground object
40objectGround = mbs.CreateGround(referencePosition = [0,0,0])
41
42## create mass point
43massPoint = mbs.CreateMassPoint(referencePosition=[L,0,0],
44                    initialDisplacement=[u0,0,0],
45                    initialVelocity=[v0,0,0],
46                    physicsMass=mass)
47
48## create spring damper  between objectGround and massPoint
49mbs.CreateCartesianSpringDamper(bodyNumbers=[objectGround, massPoint],
50                                stiffness = [k,k,k],
51                                damping   = [d,0,0],
52                                offset    = [L,0,0])
53
54## create force vector [f,0,0]
55mbs.CreateForce(bodyNumber=massPoint, loadVector= [f,0,0])
56
57## assemble and solve system
58mbs.Assemble()
59
60simulationSettings = exu.SimulationSettings()
61
62tEnd = 1
63steps = 1000000
64simulationSettings.timeIntegration.numberOfSteps = steps
65simulationSettings.timeIntegration.endTime = tEnd
66simulationSettings.timeIntegration.newton.useModifiedNewton = True
67simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1 #SHOULD work with 0.9 as well
68
69## solve and read displacement at end time
70mbs.SolveDynamic(simulationSettings)
71uCSD = mbs.GetObjectOutputBody(massPoint, exu.OutputVariableType.Displacement)[0]
72
73## compute exact (analytical) solution:
74import numpy as np
75import matplotlib.pyplot as plt
76
77# omega0 = np.sqrt(k/mass)
78# dRel = d/(2*np.sqrt(k*mass))
79
80omega = omega0*np.sqrt(1-dRel**2)
81C1 = u0-x0 #static solution needs to be considered!
82C2 = (v0+omega0*dRel*C1) / omega #C1 used instead of classical solution with u0, because x0 != 0 !!!
83
84refSol = np.zeros((steps+1,2))
85for i in range(0,steps+1):
86    t = tEnd*i/steps
87    refSol[i,0] = t
88    refSol[i,1] = np.exp(-omega0*dRel*t)*(C1*np.cos(omega*t) + C2*np.sin(omega*t))+x0
89
90print('refSol=',refSol[steps,1])
91print('error exact-numerical=',refSol[steps,1] - uCSD)
92
93## compare Exudyn with analytical solution:
94mbs.PlotSensor(['coordinatesSolution.txt', refSol],
95                components=[0,0], yLabel='displacement',
96                labels=['Exudyn','analytical'])