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 *
20
21## set up system 'mbs'
22SC = exu.SystemContainer()
23mbs = SC.AddSystem()
24
25## define overall parameters for linear spring-damper
26L=0.5
27mass = 1.6
28k = 4000
29omega0 = 50 # sqrt(4000/1.6)
30dRel = 0.05
31d = dRel * 2 * 80 #80=sqrt(1.6*4000)
32u0=-0.08
33v0=1
34f = 80
35#static result = f/k = 0.01;
36x0 = f/k
37
38## create ground object
39objectGround = mbs.CreateGround(referencePosition = [0,0,0])
40
41## create mass point
42massPoint = mbs.CreateMassPoint(referencePosition=[L,0,0],
43                    initialDisplacement=[u0,0,0],
44                    initialVelocity=[v0,0,0],
45                    physicsMass=mass)
46
47## create spring damper  between objectGround and massPoint
48mbs.CreateCartesianSpringDamper(bodyList=[objectGround, massPoint],
49                                stiffness = [k,k,k],
50                                damping   = [d,0,0],
51                                offset    = [L,0,0])
52
53## create force vector [f,0,0]
54mbs.CreateForce(bodyNumber=massPoint, loadVector= [f,0,0])
55
56## assemble and solve system
57mbs.Assemble()
58
59simulationSettings = exu.SimulationSettings()
60
61tEnd = 1
62steps = 1000000
63simulationSettings.timeIntegration.numberOfSteps = steps
64simulationSettings.timeIntegration.endTime = tEnd
65simulationSettings.timeIntegration.newton.useModifiedNewton = True
66simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1 #SHOULD work with 0.9 as well
67
68## solve and read displacement at end time
69mbs.SolveDynamic(simulationSettings)
70uCSD = mbs.GetObjectOutputBody(massPoint, exu.OutputVariableType.Displacement)[0]
71
72## compute exact (analytical) solution:
73import numpy as np
74import matplotlib.pyplot as plt
75
76# omega0 = np.sqrt(k/mass)
77# dRel = d/(2*np.sqrt(k*mass))
78
79omega = omega0*np.sqrt(1-dRel**2)
80C1 = u0-x0 #static solution needs to be considered!
81C2 = (v0+omega0*dRel*C1) / omega #C1 used instead of classical solution with u0, because x0 != 0 !!!
82
83refSol = np.zeros((steps+1,2))
84for i in range(0,steps+1):
85    t = tEnd*i/steps
86    refSol[i,0] = t
87    refSol[i,1] = np.exp(-omega0*dRel*t)*(C1*np.cos(omega*t) + C2*np.sin(omega*t))+x0
88
89print('refSol=',refSol[steps,1])
90print('error exact-numerical=',refSol[steps,1] - uCSD)
91
92## compare Exudyn with analytical solution:
93mbs.PlotSensor(['coordinatesSolution.txt', refSol],
94                components=[0,0], yLabel='displacement',
95                labels=['Exudyn','analytical'])