-
Notifications
You must be signed in to change notification settings - Fork 21
/
springDamperUserFunctionTest.py
139 lines (105 loc) · 5.08 KB
/
springDamperUserFunctionTest.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# This is an EXUDYN example
#
# Details: Test with user-defined load function and user-defined spring-damper function (Duffing oscillator)
#
# Author: Johannes Gerstmayr
# Date: 2019-11-15
#
# 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.
#
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
import exudyn as exu
from exudyn.itemInterface import *
import numpy as np
useGraphics = True #without test
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#you can erase the following lines and all exudynTestGlobals related operations if this is not intended to be used as TestModel:
try: #only if called from test suite
from modelUnitTests import exudynTestGlobals #for globally storing test results
useGraphics = exudynTestGlobals.useGraphics
except:
class ExudynTestGlobals:
pass
exudynTestGlobals = ExudynTestGlobals()
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SC = exu.SystemContainer()
mbs = SC.AddSystem()
exu.Print('EXUDYN version='+exu.GetVersionString())
L=0.5
mass = 1.6 #mass in kg
spring = 4000 #stiffness of spring-damper in N/m
damper = 4 #damping constant in N/(m/s)
load0 = 80
omega0=np.sqrt(spring/mass)
f0 = 0.*omega0/(2*np.pi)
f1 = 1.*omega0/(2*np.pi)
exu.Print('resonance frequency = '+str(omega0))
tEnd = 50 #end time of simulation
steps = 5000 #number of steps
#tEnd = 0.05 #end time of simulation
#steps = 5 #number of steps
#user function for spring force
def springForce(mbs, t, itemIndex, u, v, k, d, offset): #changed 2023-01-21:, mu, muPropZone):
return 0.1*k*u+k*u**3+v*d
#linear frequency sweep in time interval [0, t1] and frequency interval [f0,f1];
def Sweep(t, t1, f0, f1):
k = (f1-f0)/t1
return np.sin(2*np.pi*(f0+k*0.5*t)*t) #take care of factor 0.5 in k*0.5*t, in order to obtain correct frequencies!!!
#user function for load
def userLoad(mbs, t, load):
#return load*np.sin(0.5*omega0*t) #gives resonance
#exu.Print(t)
return load*Sweep(t, tEnd, f0, f1)
#return load*Sweep(t, tEnd, f1, f0) #backward sweep
#node for 3D mass point:
n1=mbs.AddNode(Point(referenceCoordinates = [L,0,0]))
#ground node
nGround=mbs.AddNode(NodePointGround(referenceCoordinates = [0,0,0]))
#add mass point (this is a 3D object with 3 coordinates):
massPoint = mbs.AddObject(MassPoint(physicsMass = mass, nodeNumber = n1))
#marker for ground (=fixed):
groundMarker=mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= nGround, coordinate = 0))
#marker for springDamper for first (x-)coordinate:
nodeMarker =mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= n1, coordinate = 0))
#Spring-Damper between two marker coordinates
mbs.AddObject(CoordinateSpringDamper(markerNumbers = [groundMarker, nodeMarker],
stiffness = spring, damping = damper,
springForceUserFunction = springForce))
#add load:
loadC = mbs.AddLoad(LoadCoordinate(markerNumber = nodeMarker,
load = load0, loadUserFunction=userLoad))
writeSensorFile = False
if useGraphics:
writeSensorFile = True
sLoad=mbs.AddSensor(SensorLoad(loadNumber=loadC, writeToFile = writeSensorFile,
storeInternal=True,#fileName="solution/userFunctionLoad.txt"
))
#mbs.AddSensor(SensorNode(nodeNumber=n1, writeToFile = writeSensorFile, fileName="solution/userFunctionNode.txt"))
sCoords=mbs.AddSensor(SensorNode(nodeNumber=n1, writeToFile = writeSensorFile,
outputVariableType=exu.OutputVariableType.Coordinates,
storeInternal=True,#fileName="solution/userFunctionNode.txt"
))
#exu.Print(mbs)
mbs.Assemble()
simulationSettings = exu.SimulationSettings()
simulationSettings.solutionSettings.solutionWritePeriod = 2e-3 #output interval
simulationSettings.timeIntegration.numberOfSteps = steps
simulationSettings.timeIntegration.endTime = tEnd
simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1
simulationSettings.displayStatistics = True
simulationSettings.timeIntegration.verboseMode = 1
#exu.StartRenderer() #start graphics visualization
#mbs.WaitForUserToContinue() #wait for pressing SPACE bar to continue
#start solver:
mbs.SolveDynamic(simulationSettings)
#SC.WaitForRenderEngineStopFlag()#wait for pressing 'Q' to quit
#exu.StopRenderer() #safely close rendering window!
#evaluate final (=current) output values
u = mbs.GetNodeOutput(n1, exu.OutputVariableType.Position)
exu.Print('displacement=',u[0])
exudynTestGlobals.testError = u[0] - (0.5062872273010898) #2019-12-18: 0.5062872273010898; #2019-12-15: 0.5062872272996835; 2019-12-13:0.5062872273014417; 2019-12-01: 0.5152217339585201
exudynTestGlobals.testResult = u[0]
#+++++++++++++++++++++++++++++++++++++++++++++++++++++
if useGraphics:
mbs.PlotSensor(sCoords, components=[0], closeAll=True)