## 15.4. Harmonic Analysis

The harmonic analysis (ANTYPE,HARMIC) solves the time-dependent equations of motion (Equation 15–5) for linear structures undergoing steady-state vibration.

### 15.4.1. Harmonic Analysis Assumptions and Restrictions

1. Valid for structural, fluid, magnetic, and electrical degrees of freedom (DOFs). Thermal DOFs may be present in a coupled field harmonic analysis using structural DOFs.

2. The entire structure has constant or frequency-dependent stiffness, damping, and mass effects.

3. All loads and displacements vary sinusoidally at the same known frequency (although not necessarily in phase).

4. Element loads are assumed to be real (in-phase) only, except for:

1. current density

2. pressures in SURF153, SURF154, SURF156, and SURF159

### 15.4.2. Description of Harmonic Analysis

Consider the general equation of motion for a structural system (Equation 15–5).

 (15–52)

where:

 = structural mass matrix = structural damping matrix = structural stiffness matrix = nodal acceleration vector = nodal velocity vector = nodal displacement vector = applied load vector

As stated above, all points in the structure are moving at the same known frequency, however, not necessarily in phase. Also, it is known that the presence of damping causes phase shifts. Therefore, the displacements may be defined as:

 (15–53)

where:

 = maximum displacement = imposed circular frequency (radians/time) = = imposed frequency (cycles/time) (input as FREQB and FREQE on the HARFRQ command) = time = displacement phase shift (radians)

Note that and may be different at each DOF. The use of complex notation allows a compact and efficient description and solution of the problem. Equation 15–53 can be rewritten as:

 (15–54)

or as:

 (15–55)

where:

 = real displacement vector (input as VALUE on D command, when specified) = imaginary displacement vector (input as VALUE2 on D command, when specified)

The force vector can be specified analogously to the displacement:

 (15–56)
 (15–57)
 (15–58)

where:

 = force amplitude = force phase shift (radians) = real force vector (input as VALUE on F command, when specified) = imaginary force vector (input as on VALUE2 on F command, when specified)

Substituting Equation 15–55 and Equation 15–58 into Equation 15–52 gives:

 (15–59)

The dependence on time () is the same on both sides of the equation and may therefore be removed:

 (15–60)

The solution of this equation is discussed later.

### 15.4.3. Harmonic Analysis Complex Displacement Output

The complex displacement output at each DOF may be given in one of two forms:

1. The same form as and as defined in Equation 15–55 (selected with the command HROUT,ON).

2. The form and (amplitude and phase angle (in degrees)), as defined in Equation 15–54 (selected with the command HROUT,OFF). These two terms are computed at each DOF as:

 (15–61)
 (15–62)

Note that the response lags the excitation by a phase angle of .

### 15.4.4. Nodal and Reaction Load Computation in a Harmonic Analysis

Inertia, damping and static loads on the nodes of each element are computed.

The inertia load is expressed as:

 (15–63)

where:

 = vector of element inertia forces = element mass matrix = element displacement vector

The real and imaginary inertia load parts of the element output are then computed by:

 (15–64)
 (15–65)

where:

 = vector of element inertia forces (real part) = element real displacement vector = vector element inertia forces (imaginary part) = element imaginary displacement vector

The damping load is expressed as:

 (15–66)

where:

 = vector of element damping forces = element damping matrix

The real and imaginary damping loads part of the element output are then computed by:

 (15–67)
 (15–68)

where:

 = vector of element damping forces (real part) = vector of element damping forces (imaginary part)

The real static load is computed the same way as in a static analysis (Solving for Unknowns and Reactions) using the real part of the displacement solution . The imaginary static load is computed also the same way, using the imaginary part . Note that the imaginary part of the element loads (e.g., ) are normally zero, except for current density loads.

The nodal reaction loads are computed as the sum of all three types of loads (inertia, damping, and static) over all elements connected to a given fixed displacement node.

### 15.4.5. Harmonic Analysis Solution

Three methods of solution to Equation 15–60 are available:

#### 15.4.5.1. Full Solution Method

The full solution method (HROPT,FULL) solves Equation 15–60 directly. Equation 15–60 may be expressed as:

 (15–69)

where denotes a complex matrix or vector. Equation 15–69 is solved using the same sparse solver used for a static analysis in Equation Solvers, except that it is done using complex arithmetic.

#### 15.4.5.2. Mode-Superposition Method

The mode-superposition method (HROPT,MSUP) uses the natural frequencies and mode shapes to compute the response to a sinusoidally varying forcing load. This solution method imposes the following additional assumptions and restrictions:

1. Nonzero imposed harmonic displacements are not allowed.

2. There are no element damping matrices. However, various types of system damping are available.

The development of the general mode-superposition procedure is given in Mode-Superposition Method. The equation of motion (Equation 15–52) is converted to modal form, as described in Mode-Superposition Method. Equation 14–122 is:

 (15–70)

where:

 = modal coordinate = natural circular frequency of mode j = fraction of critical damping for mode j = force in modal coordinates

The load vector which is converted to modal coordinates (Equation 14–121) is given by

 (15–71)

where:

 = nodal force vector = load vector scale factor, (input as FACT, LVSCALE command) = load vector from the modal analysis (see Mode-Superposition Method).

For a steady sinusoidal vibration, has the form

 (15–72)

where:

 = complex force amplitude = imposed circular frequency

For Equation 15–70 to be true at all times, must have a similar form as , or

 (15–73)

where:

 = complex amplitude of the modal coordinate for mode j.

Differentiating Equation 15–73, and substituting Equation 15–72 and Equation 15–73 into Equation 15–70,

 (15–74)

Collecting coefficients of and dividing by ()

 (15–75)

solving for ,

 (15–76)

The contribution from each mode is:

 (15–77)

where:

 = contribution of mode j (output if Mcont = ON, on the HROUT command) = mode shape for mode j

Finally, the complex displacements are obtained from Equation 14–104 as

 (15–78)

where:

 = vector of complex displacements

In the case of the QR damped mode extraction method, one substitutes Equation 14–123 for Equation 14–122, so Equation 15–75 becomes:

 (15–79)

Solving the above equation and multiplying by the eigenvectors, one can calculate the complex displacements shown in Equation 15–78.

##### 15.4.5.2.1. Expansion Pass

The expansion pass of the mode-superposition harmonic analysis involves computing element stresses, element nodal forces, reaction forces, etc.

### 15.4.6. Variational Technology Method

A common way to compute the harmonic response of a structure is to compute the normal modes of the undamped structure, and to use a modal superposition method to evaluate the response, after determining the modal damping. Determining the modal damping can be based on modal testing, or by using empirical rules. However, when the structure is non-metallic, the elastic properties can be highly dependent on the frequency and the damping can be high enough that the undamped modes and the damped modes are significantly different, and an approach based on a real, undamped modes is not appropriate.

One alternative to straight modal analysis is to build multiple modal bases, for different property values, and combine them together over the frequency range of the analysis. This technique is complex, error prone, and does not address the problem of determining the modal damping factors. Another alternative is a direct frequency response, updating the elastic properties for every frequency step. This technique give a much better prediction of the frequency response, but is CPU intensive.

The variational technology method (HROPT,VT) is available as the harmonic sweep capability of the VT Accelerator add-on. You can define the material elastic properties as being frequency-dependent (using TB,ELASTIC and TB,SDAMP) and efficiently compute the frequency response over an entire frequency range. For the Variational Technology theory, see Guillaume([334]) and Beley, Broudiscou, et al.([361]).

#### 15.4.6.1. Viscous or Hysteretic Damping

When using the Variational Technology method, the user has a choice between viscous and hysteretic damping.

Viscous Damping

Consider a spring-damper-mass system subjected to a harmonic excitation. The response of the system is given by:

 (15–80)

Due to the damping, the system is not conservative and the energy is dissipated. Using viscous damping, the energy dissipated by the cycle is proportional to the frequency, . In a single DOF spring-mass-damper system, with a viscous damper :

 (15–81)

where:

 = change of energy = viscous damper

Hysteretic Damping

Experience shows that energy dissipated by internal friction in a real system does not depend on frequency, and approximately is a function of :

 (15–82)

where:

 = frequency-dependent damping (input using TB,SDAMP command)

damping is known as structural or hysteretic damping. It can be included in the elastic properties by using a complex Young's modulus:

 (15–83)

where:

 = Young's modulus (input using TB,ELASTIC command)

Using this kind of representation, the equations of motion of the system become:

 (15–84)

where:

 = structural mass matrix = structural stiffness matrix = structural damping matrix = nodal acceleration vector = nodal displacement vector = applied load vector

Note:  The Variational Technology method does not support a material with both viscous and hysteretic damping properties in the same model. In the case, please use the HROPT,FULL command to analyze the material.

### 15.4.7. Automatic Frequency Spacing in a Harmonic Analysis

In a harmonic analysis, the imposed frequencies that involve the most kinetic energy are those near the natural frequencies of the structure. The automatic frequency spacing or "cluster" option (Clust = ON, on the HROUT command) provides an approximate method of choosing suitable imposed frequencies. The nearness of the imposed frequencies to the natural frequencies depends on damping, because the resonance peaks narrow when the damping is reduced. Figure 15.3: Frequency Spacing shows two typical resonance peaks and the imposed frequencies chosen by this method, which are computed from:

 (15–85)
 (15–86)

Equation 15–85 gives frequencies slightly less than the natural circular frequency . Equation 15–86 gives slightly higher frequencies. The spacing parameter is defined as:

 (15–87)

where:

 = modal damping. (If is computed as 0.0, it is redefined to be 0.005 for this equation only). = integer constant (input as NSBSTP, NSUBST command) which may be between 2 and 20. Anything above this range defaults to 10 and anything below this range defaults to 4. = 1, 2, 3, ... N

Each natural frequency, as well as frequencies midway between, are also chosen as imposed frequencies.

### 15.4.8. Logarithm Frequency Spacing in a Harmonic Analysis

The logarithm frequency spacing can be defined over an nth-octave band or a general-frequency band in a harmonic analysis (the HARFRQ command with `LogOpt`).

The octave band is defined by the lower-end frequency:

 (15–88)

and upper-end frequency:

 (15–89)

where the central frequency of the octave band is specified at 16, 31.6, 63, 125, 250, 500, 1000, 2000, 4000, 8000, 16000 Hz for 11 octave bands, and n = 1 (octave band), 2 (1/2 octave band), 3 (1/3 octave band), 6 (1/6 octave band), 12 (1/12 octave band), and 24 (1/24 octave band).

The recursive relationships between the central frequencies of the specified octave bands are:

 (15–90)
 (15–91)

The logarithm frequency increment is given by:

 (15–92)

Where `NSBSTP` is input by the NSUBST command.

The imposed working frequency is calculated by:

 (15–93)

The imposed working frequency is , if `NSBSTP` = 1. For a general logarithm frequency span, the beginning frequency and ending frequency are defined by the FREQB and FREQE parameters of the HARFRQ command respectively.

### 15.4.9. Harmonic Analysis with Rotating Forces on Rotating Structures

If a structure is rotating, forces rotating synchronously or asynchronously with the structure are of interest.

General rotating asynchronous forces are described in General Asynchronous Rotating Force. A specific synchronous force: mass unbalance is shown in Specific Synchronous Forces: Mass Unbalance.

In both cases, the equation solved for harmonic analysis is the same as (Equation 15–60) except for the coefficients of the damping matrix which will be a function of the rotational velocity of the structure (see the CORIOLIS command). will be updated for each excitation frequency step using the following rotational velocity:

 (15–94)

where:

 = rotational velocity of the structure (rd/s) = frequency of excitation (rd/s) = ratio between and (s = 1 for synchronous excitations) (input as `RATIO` in the SYNCHRO command).

The right-hand term of the equation is given below depending on the force considered.

#### 15.4.9.1. General Asynchronous Rotating Force

If the structure is rotating about X axis, then an asynchronous force having its direction in the plane perpendicular to the spin axis is expressed as:

 (15–95)
 (15–96)

where:

 = amplitude of force

Using complex notations, the equations become:

 (15–97)
 (15–98)

where:

 = real force in Y-direction; also, negative of imaginary force in Z-direction. Input as `VALUE` on F, label FY; input as `VALUE2` on F, label FZ = negative of real force in Z-direction; also, negative of imaginary force in Y-direction. Input as `VALUE` on F, label FZ; input as `VALUE2` on F, label FY

The expression of the forces for structures rotating about another direction than X are developed analogously.

#### 15.4.9.2. Specific Synchronous Forces: Mass Unbalance

Consider a structure rotating about the X axis. The mass unbalance m situated at node I with the eccentricity e may be represented as shown in Figure 15.4: Mass Unbalance at Node I

If we only consider the motion in the plane perpendicular to the spin axis (YZ plane), the kinetic energy of the unbalanced mass is written as:

 (15–99)

where:

 = mass unbalance = distance from the mass unbalance to the spin axis = amplitude of the rotational velocity vector of the structure (input as OMEGA or CMOMEGA command). It is equal to the frequency of excitation . = phase of the unbalance = instantaneous velocity along Y = instantaneous velocity along Z

Because the mass unbalance is much smaller than the weight of the structure, the first two terms are neglected. The third term being constant, will have no effect on the final equations.

Applying Lagrange's equations, the force vector is:

 (15–100)
 (15–101)

where:

Using complex notations, it can be written as:

 (15–102)
 (15–103)

Note:  The multiplication of the forces by is done internally at each frequency step.

### 15.4.10. Harmonic Ocean Wave Procedure (HOWP)

Standard harmonic analysis with ocean loading includes the added mass of the water outside of line elements (such as LINK180, BEAM188, BEAM189, PIPE288, and PIPE289). The load vector is calculated based on the loads at a given time, but the standard analysis method usually misses some important effects, as all peak wave loads rarely occur at the same time.

All relevant ocean wave loading effects are accounted for via a specialized variation of the harmonic analysis, the harmonic ocean wave procedure (HOWP).

The HOWP applies to regular waves only (Airy and Wheeler one-component waves, as well as Stokes and Deans Stream Function waves), and works only with the full-solution harmonic analysis method (HROPT,FULL).

The HOWP is accessed via the HROCEAN command. The frequency is obtained automatically, directly from the ocean information (via OCDATA and OCTABLE commands). As with a standard harmonic analysis, a damping matrix must be added separately if one is needed. Ocean loads are calculated with the assumption that the structure is stationary.

For every component of every element, HOWP-specific calculations obtain the real and imaginary loads caused by ocean waves. All other loads (such as those related to pressure, temperature, gravity and buoyancy) are removed. Stated differently, all loads are sinusoidal at the given frequency, but the magnitude and phase angle of every component of every element must be determined. To determine the magnitude and phase angle, a series of static analyses is performed along the phase angle ranging from 0 to 360 degrees. (The number of these analyses is controlled via the HROCEAN command.) The result is a roughly sinusoidal force pattern, as shown in this figure:

From this force history, the maximum and minimum forces are calculated by fitting the highest and lowest points. Then:

 (15–104)
 (15–105)
 (15–106)

so that is the coefficient on the real load vector and is the coefficient on the imaginary load vector.