pyyeti.ode.solvepsd

pyyeti.ode.solvepsd(fs, forcepsd, t_frc, freq, drmlist, rbduf=1.0, elduf=1.0, **kwargs)

Solve equations of motion in frequency domain with uncorrelated PSD forces.

See also pyyeti.cla.DR_Results.solvepsd() for a very similar routine, but one that is designed for use within the pyYeti “coupled loads analysis paradigm” (where the classes defined in pyyeti.cla are used).

Parameters:

• fs (class instance) – An instance of SolveUnc or FreqDirect (or similar … must have .fsolve method)

• forcepsd (2d array_like) – The matrix of force psds; each row is a force PSD

• t_frc (2d array_like) – Transform to put forcepsd into the coordinates of the equations of motion: t_frc @ forcepsd. Commonly, t_frc is simply the transpose of a row-partition of the mode shape matrix (phi) and the conversion of forcepsd is from physical space to modal space. In that case, the row-partition is from the full set down to just the forced DOF. However, t_frc can also have force mappings (as from the TLOAD capability in Nastran); in that case, t_frc = phi.T @ mapping_vectors. In any case, the number of columns in t_frc is the number of rows in forcepsd: t_frc.shape[1] == forcepsd.shape[0]

• freq (1d array_like) – Frequency vector at which solution will be computed; len(freq) = cols(forcepsd)

• drmlist (list_like) – List of lists (or similar) of any number of data recovery matrix quadruples (in the order typically used to write equations of motion):

  [
      [drma1, drmv1, drmd1, drmf1],
      [drma2, drmv2, drmd2, drmf2],
      ...
  ]
  

  To not use a particular drm, set it to None. For example, to perform these 3 types of data recovery:

  acce = atm*a
  disp = dtm*d
  loads = ltma*a + ltmv*v + ltmd*d + ltmf*f
  

  drmlist would be:

  [[atm, None, None, None],
   [None, None, dtm, None],
   [ltma, ltmv, ltmd, ltmf]]
  

• rbduf (scalar; optional) – Rigid-body uncertainty factor

• elduf (scalar; optional) – Dynamic uncertainty factor

• **kwargs (keyword arguments for fs.fsolve; optional) – Currently, there are two arguments available:

  argument	brief description
  incrb	specifies how to handle rigid-body responses
  rf_disp_only	specifies how to handle residual-flexibility modes

  See, for example, SolveUnc.fsolve().

Returns:

• rms (list) – List of vectors (corresponding to drmlist) of rms values of all data recovery rows; # of rows of each vector = # rows in each drm pair

• psd (list) – List of matrices (corresponding to drmlist) of PSD responses for all data recovery rows:

  # rows in each PSD = # rows in DRM
  # cols in each PSD = len(freq)
  

Notes

This routine first calls fs.fsolve to solve the modal equations of motion. Then, it scales the responses by the corresponding PSD input force. All PSD responses are summed together for the overall response. For example:

resp_psd = 0
for i in range(forcepsd.shape[0]):
    # solve for unit frequency response function:
    unitforce = np.ones((1, len(freq)))
    genforce = t_frc[:, i:i+1] @ unitforce
    sol = fs.fsolve(genforce, freq, incrb="av")
    frf = (drma @ sol.a
           + drmv @ sol.v
           + drmd @ sol.d
           + drmf[:, [i]] @ unitforce)
    resp_psd += forcepsd[i] * abs(frf)**2

In that example, the data recovery uses all four drms. Also, the looping over the drmlist is not included for simplicity.

Examples

>>> from pyyeti import ode
>>> import numpy as np
>>> m = np.array([10., 30., 30., 30.])        # diag mass
>>> k = np.array([0., 6.e5, 6.e5, 6.e5])      # diag stiffness
>>> zeta = np.array([0., .05, 1., 2.])        # % damping
>>> b = 2.*zeta*np.sqrt(k/m)*m                # diag damping
>>> freq = np.arange(.1, 35, .1)              # frequency
>>> forcepsd = 10000*np.ones((4, freq.size))  # PSD forces
>>> fs = ode.SolveUnc(m, b, k)
>>> atm = np.eye(4)    # recover modal accels
>>> t_frc = np.eye(4)  # assume forces already modal
>>> drms = [[atm, None, None, None]]
>>> rms, psd = ode.solvepsd(fs, forcepsd, t_frc, freq,
...                         drms)

The rigid-body results should be 100.0 g**2/Hz flat; rms = np.sqrt(100*34.8)

>>> np.allclose(100., psd[0][0])
True
>>> np.allclose(np.sqrt(3480.), rms[0][0])
True

Plot the four accelerations PSDs:

>>> import matplotlib.pyplot as plt
>>> fig = plt.figure('Example', figsize=[8, 8], clear=True,
...                  layout='constrained')
>>> labels = ['Rigid-body', 'Underdamped',
...           'Critically Damped', 'Overdamped']
>>> for j, name in zip(range(4), labels):
...     _ = plt.subplot(4, 1, j+1)
...     _ = plt.plot(freq, psd[0][j])
...     _ = plt.title(name)
...     _ = plt.ylabel(r'Accel PSD ($g^2$/Hz)')
...     _ = plt.xlabel('Frequency (Hz)')