Coverage for src / sdynpy / fem / sdynpy_shaker.py: 26%
43 statements
« prev ^ index » next coverage.py v7.13.4, created at 2026-03-11 16:22 +0000
1# -*- coding: utf-8 -*-
2"""
3Tools for modeling shakers, and includes commercial shaker objects.
4"""
5"""
6Copyright 2022 National Technology & Engineering Solutions of Sandia,
7LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S.
8Government retains certain rights in this software.
10This program is free software: you can redistribute it and/or modify
11it under the terms of the GNU General Public License as published by
12the Free Software Foundation, either version 3 of the License, or
13(at your option) any later version.
15This program is distributed in the hope that it will be useful,
16but WITHOUT ANY WARRANTY; without even the implied warranty of
17MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18GNU General Public License for more details.
20You should have received a copy of the GNU General Public License
21along with this program. If not, see <https://www.gnu.org/licenses/>.
22"""
24import numpy as np
25import matplotlib.pyplot as plt
27# def shaker_4dof(m_armature,m_body,m_forcegauge,
28# k_suspension,k_stinger,
29# c_suspension,c_stinger,
30# resistance,inductance,force_factor):
31# '''Creates a four degree of freedom electromechanical model of a shaker
33# Parameters
34# ----------
35# m_armature : float
36# Mass of the armaturesistance
37# m_body : float
38# Mass of the body and trunion of the shaker
39# m_forcegauge : float
40# Mass of the force gauge at the tip of the stinger
41# k_suspension : float
42# Stiffness of the spring between the body and armature
43# k_stinger : float
44# Stiffness of the stinger between the armature and force gauge
45# c_suspension : float
46# Damping of the spring between the body and armature
47# c_stinger : float
48# Damping of the stinger between the armature and force gauge
49# resistance : float
50# Coil resistance of the electronics portion of the shaker
51# inductance : float
52# Coil inductance of the electronics portion of the shaker
53# force_factor : float
54# Force factor BL of the magnet and coil
56# Returns
57# -------
58# M_shaker : ndarray
59# The mass matrix of the shaker
60# C_shaker : ndarray
61# The damping matrix of the shaker
62# K_shaker : ndarray
63# The stiffness matrix of the shaker
65# Notes
66# -----
67# The dof order of the responses are:
68# Displacement of Armature,
69# Displacement of Body,
70# Displacement of Force Gauge,
71# Shaker Current.
73# The dof order of the inputs are:
74# Force into Armature
75# Force into Body
76# Force into Force Gauge
77# Shaker Voltage.
79# The force imparted by the shaker can be computed as the stinger stiffness
80# times the difference in displacement of the armature and force gauge plus
81# the stinger damping times the difference in the velocity of the armature
82# and force gauge.
83# '''
85# # Create the matrices
86# M_shaker = np.array([[m_armature, 0, 0, 0],
87# [ 0, m_body, 0, 0],
88# [ 0, 0, m_forcegauge, 0],
89# [ 0, 0, 0, 0]])
90# C_shaker = np.array([[(c_suspension+c_stinger), -c_suspension, -c_stinger, 0],
91# [ -c_suspension, c_suspension, 0, 0],
92# [ -c_stinger, 0, c_stinger, 0],
93# [ force_factor, -force_factor, 0, inductance]])
94# K_shaker = np.array([[(k_suspension+k_stinger), -k_suspension, -k_stinger, -force_factor],
95# [ -k_suspension, k_suspension, 0, force_factor],
96# [ -k_stinger, 0, k_stinger, 0],
97# [ 0, 0, 0, resistance]])
99# return M_shaker,C_shaker,K_shaker
101# test_shaker_dakota = shaker_4dof( m_armature = 0.227, # [kg] ~ armaturesistance mass
102# m_body = 24.9, # [kg] ~ body mass
103# k_suspension = 2312, # [N/m] ~ suspension stiffness
104# c_suspension = 20, # [N/(m/s)] damping ...
105# k_stinger = 4e6, # [N/m] stiff stinger spring
106# m_forcegauge = 1e-2, # [kg] force gauge mass
107# c_stinger = 10, # [N/(m/s)] damping
108# resistance = 1.3, # [ohm] ~ coil resistance
109# inductance = 100e-6, # [Henry] ~ coil inductance
110# force_factor = 30 # [-] ~ Force Factor of magnet & coil
111# )
114class Shaker4DoF:
115 """A class defining a four degree-of-freedom model of a shaker
117 """
119 def __init__(self, m_armature, m_body, m_forcegauge,
120 k_suspension, k_stinger,
121 c_suspension, c_stinger,
122 resistance, inductance, force_factor):
123 '''Creates a four degree of freedom electromechanical model of a shaker
125 Parameters
126 ----------
127 m_armature : float
128 Mass of the armature
129 m_body : float
130 Mass of the body and trunion of the shaker
131 m_forcegauge : float
132 Mass of the force gauge at the tip of the stinger
133 k_suspension : float
134 Stiffness of the spring between the body and armature
135 k_stinger : float
136 Stiffness of the stinger between the armature and force gauge
137 c_suspension : float
138 Damping of the spring between the body and armature
139 c_stinger : float
140 Damping of the stinger between the armature and force gauge
141 resistance : float
142 Coil resistance of the electronics portion of the shaker
143 inductance : float
144 Coil inductance of the electronics portion of the shaker
145 force_factor : float
146 Force factor BL of the magnet and coil
147 '''
148 self.m_armature = m_armature
149 self.m_body = m_body
150 self.m_forcegauge = m_forcegauge
151 self.k_suspension = k_suspension
152 self.k_stinger = k_stinger
153 self.c_suspension = c_suspension
154 self.c_stinger = c_stinger
155 self.resistance = resistance
156 self.inductance = inductance
157 self.force_factor = force_factor
159 def MCK(self):
160 """
161 Returns mass, damping, and stiffness matrices for the shakers.
163 Returns
164 -------
165 M : np.ndarray
166 The mass matrix of the shaker
167 C : np.ndarray
168 The damping matrix of the shaker
169 K : np.ndarray
170 The stiffness matrix of the shaker
172 Notes
173 -----
174 The dof order of the responses are:
175 Displacement of Armature,
176 Displacement of Body,
177 Displacement of Force Gauge,
178 Shaker Current.
180 The dof order of the inputs are:
181 Force into Armature
182 Force into Body
183 Force into Force Gauge
184 Shaker Voltage.
186 The force imparted by the shaker can be computed as the stinger stiffness
187 times the difference in displacement of the armature and force gauge plus
188 the stinger damping times the difference in the velocity of the armature
189 and force gauge.
191 """
192 M = np.array([[self.m_armature, 0, 0, 0], [0, self.m_body, 0, 0],
193 [0, 0, self.m_forcegauge, 0], [0, 0, 0, 0]])
194 C = np.array([[(self.c_suspension + self.c_stinger), -self.c_suspension, -self.c_stinger, 0],
195 [-self.c_suspension, self.c_suspension, 0, 0],
196 [-self.c_stinger, 0, self.c_stinger, 0],
197 [self.force_factor, -self.force_factor, 0, self.inductance]])
198 K = np.array([[(self.k_suspension + self.k_stinger), -self.k_suspension, -self.k_stinger, -self.force_factor],
199 [-self.k_suspension, self.k_suspension, 0, self.force_factor],
200 [-self.k_stinger, 0, self.k_stinger, 0],
201 [0, 0, 0, self.resistance]])
202 return M, C, K
204 def state_space(self):
205 """
206 Returns the state space formulation of the shaker model of the form
208 x_dot = A@x + B@u
209 y = C@x + D@u
211 Returns
212 -------
213 A : np.ndarray
214 State (system) Matrix.
215 B : np.ndarray
216 Input Matrix
217 C : np.ndarray
218 Output Matrix
219 D : np.ndarray
220 Feedthrough (feedforward) Matrix
222 Notes
223 -----
224 The dof order of the state matrix is
225 0. armature displacement
226 1. body displacement
227 2. forcegauge displacement,
228 3. armature velocity
229 4. body velocity
230 5. forcegauge velocity
231 6. current.
233 The dof order of the input matrix is
234 0. external force on armature
235 1. external force on body
236 2. external force on forcegauge
237 3. shaker voltage
239 The dof order for the output matries is
240 0. Displacement of armature
241 1. Velocity of armature
242 2. Acceleration of armature
243 3. External force on armature
244 4. Displacement of body
245 5. Velocity of body
246 6. Acceleration of body
247 7. External force on body
248 8. Displacement of forcegauge
249 9. Velocity of forcegauge
250 10. Acceleration of forcegauge
251 11. External force on forcegauge
252 12. Shaker voltage
253 13. Shaker current
254 14. Current change per time
255 15. Force on coil due to electronics
256 16. Force in stinger
257 """
258 A = np.array([[0, 0, 0, 1, 0, 0, 0],
259 [0, 0, 0, 0, 1, 0, 0],
260 [0, 0, 0, 0, 0, 1, 0],
261 [(-self.k_suspension - self.k_stinger) / self.m_armature,
262 self.k_suspension / self.m_armature,
263 self.k_stinger / self.m_armature,
264 (-self.c_suspension
265 - self.c_stinger) / self.m_armature,
266 self.c_suspension / self.m_armature,
267 self.c_stinger / self.m_armature,
268 self.force_factor / self.m_armature],
269 [self.k_suspension / self.m_body, -self.k_suspension / self.m_body, 0, self.c_suspension /
270 self.m_body, -self.c_suspension / self.m_body, 0, -self.force_factor / self.m_body],
271 [self.k_stinger / self.m_forcegauge, 0, -self.k_stinger / self.m_forcegauge,
272 self.c_stinger / self.m_forcegauge, 0, -self.c_stinger / self.m_forcegauge, 0],
273 [0, 0, 0, -self.force_factor / self.inductance, self.force_factor / self.inductance, 0, -self.resistance / self.inductance]])
274 B = np.array([[0, 0, 0, 0],
275 [0, 0, 0, 0],
276 [0, 0, 0, 0],
277 [1 / self.m_armature, 0, 0, 0],
278 [0, 1 / self.m_body, 0, 0],
279 [0, 0, 1 / self.m_forcegauge, 0],
280 [0, 0, 0, 1 / self.inductance]])
281 C = np.array([
282 # States: x1, x2, x3, xd1, xd2, xd3, i
283 [1, 0, 0, 0, 0, 0, 0], # Displacement of mass 1
284 [0, 0, 0, 1, 0, 0, 0], # Velocity of mass 1
285 [(-self.k_suspension - self.k_stinger) / self.m_armature, self.k_suspension / self.m_armature,
286 self.k_stinger / self.m_armature, (-self.c_suspension - self.c_stinger) / \
287 self.m_armature, self.c_suspension / self.m_armature,
288 self.c_stinger / self.m_armature, self.force_factor / self.m_armature], # Acceleration of mass 1
289 [0, 0, 0, 0, 0, 0, 0], # External Force on mass 1
290 [0, 1, 0, 0, 0, 0, 0], # Displacement of mass 2
291 [0, 0, 0, 0, 1, 0, 0], # Velocity of mass 2
292 [self.k_suspension / self.m_body, -self.k_suspension / self.m_body, 0, self.c_suspension / self.m_body, - \
293 self.c_suspension / self.m_body, 0, -self.force_factor / self.m_body], # Acceleration of mass 2
294 [0, 0, 0, 0, 0, 0, 0], # External Force on mass 2
295 [0, 0, 1, 0, 0, 0, 0], # Displacement of mass 3
296 [0, 0, 0, 0, 0, 1, 0], # Velocity of mass 3
297 [self.k_stinger / self.m_forcegauge, 0, -self.k_stinger / self.m_forcegauge, self.c_stinger / \
298 self.m_forcegauge, 0, -self.c_stinger / self.m_forcegauge, 0], # Acceleration of mass 3
299 [0, 0, 0, 0, 0, 0, 0], # External Force on mass 3
300 [0, 0, 0, 0, 0, 0, 0], # Voltage
301 [0, 0, 0, 0, 0, 0, 1], # Current
302 [0, 0, 0, -self.force_factor / self.inductance, self.force_factor / self.inductance,
303 0, -self.resistance / self.inductance], # Current change per time
304 [0, 0, 0, 0, 0, 0, self.force_factor], # Force on coil
305 [self.k_stinger, 0, -self.k_stinger, self.c_stinger,
306 0, -self.c_stinger, 0], # Force in stinger
307 ])
309 D = np.array([
310 # Inputs: f_armature, f_body, f_forcegauge, e
311 [0, 0, 0, 0], # Displacement of mass 1
312 [0, 0, 0, 0], # Velocity of mass 1
313 [1 / self.m_armature, 0, 0, 0], # Acceleration of mass 1
314 [1, 0, 0, 0], # External Force on mass 1
315 [0, 0, 0, 0], # Displacement of mass 2
316 [0, 0, 0, 0], # Velocity of mass 2
317 [0, 1 / self.m_body, 0, 0], # Acceleration of mass 2
318 [0, 1, 0, 0], # External Force on mass 2
319 [0, 0, 0, 0], # Displacement of mass 3
320 [0, 0, 0, 0], # Velocity of mass 3
321 [0, 0, 1 / self.m_forcegauge, 0], # Acceleration of mass 3
322 [0, 0, 1, 0], # External Force on mass 3
323 [0, 0, 0, 1], # Voltage
324 [0, 0, 0, 0], # Current
325 [0, 0, 0, 1 / self.inductance], # Current change per time
326 [0, 0, 0, 0], # Force on coil
327 [0, 0, 0, 0], # Force in stinger
328 ])
329 return A, B, C, D
331 def transfer_function(self, frequencies):
332 """
333 Computes the shaker transfer function matrix
335 Returns
336 -------
337 H : np.ndarray
338 A transfer function matrix.
340 """
341 angular_frequencies = 2 * np.pi * frequencies[:, np.newaxis, np.newaxis]
342 M, C, K = self.MCK()
343 H = np.linalg.inv(-angular_frequencies**2 * M + 1j * angular_frequencies * C + K)
344 return H
346 def plot_electrical_impedance(self, frequencies=np.arange(4000) + 1, test_data=None):
347 """
348 Plots the electrical impedance voltage/current for the shaker model
349 and compares to test data if supplied
351 Parameters
352 ----------
353 frequencies : np.array, optional
354 Frequency lines. The default is np.arange(4000).
355 test_data : np.ndarray, optional
356 Test data defined at frequencies. The default is None.
359 """
360 fig, ax = plt.subplots(2, 1, sharex=True)
361 H = self.transfer_function(frequencies)[..., -1, -1]
362 ax[0].plot(frequencies, np.real(1 / H))
363 if test_data is not None:
364 ax[0].plot(frequencies, np.real(test_data))
365 ax[1].plot(frequencies, np.imag(1 / H))
366 if test_data is not None:
367 ax[1].plot(frequencies, np.real(test_data))
369 @classmethod
370 def modal_shop_100lbf(cls):
371 """
372 Returns a shaker model for a Modal Shop Model 2100E11 100lbf modal shaker.
374 Returns
375 -------
376 shaker_model : Shaker4DoF
377 A 4 degree of freedom model of the Modal Shop 2100E11 100lbf modal
378 shaker.
380 Notes
381 -----
382 The shaker uses the following parameters
383 - m_armature = 0.44 kg
384 - m_body = 15 kg
385 - m_forcegauge = 1e-2 kg
386 - k_suspension = 1.1e5 N/m
387 - k_stinger = 9.63e6 N/m
388 - c_suspension = 9.6 Ns/m
389 - c_stinger = 0.42 Ns/m
390 - resistance = 4 Ohm
391 - inductance = 6e-4 Henry
392 - force_factor = 36 (-)
394 """
395 return cls(m_armature=0.44,
396 m_body=15,
397 m_forcegauge=1e-2,
398 k_suspension=1.1e5,
399 k_stinger=9.63e6,
400 c_suspension=9.6,
401 c_stinger=0.42,
402 resistance=4,
403 inductance=6e-4,
404 force_factor=36)
406# shaker = Shaker4DoF.modal_shop_100lbf()
407# shaker.m_forcegauge = 2.6
408# shaker.plot_electrical_impedance()