Description of Atlas I

A spacecraft consists of mission-specific equipment called the payload and a collection of subsystem components called a bus. Typical subsystem components may include attitude control, propulsion, communications, electrical power, etc. (Sarafin, 1995, p. 451). The spacecraft bus may be able to support different payloads for different missions. As shown in figure [fig18.1], the main components of the Atlas I space launch vehicle from bottom to top are an expendable booster, an expendable second stage called Centaur, and the payload/spacecraft (Isakowitz, 1995). The spacecraft is covered by a fairing. An adapter, also known as a launch vehicle adapter, structurally links the spacecraft and Atlas I. The spacecraft developer usually designs the adapter.

Selected weights and thrust data for the Atlas I in table [tab18.1] are from Isakowitz (1995).

Propellant weights in the booster and second stage, and the approximate structural weight are listed in table [tab18.2].

In terms of percentages

1. 92.6 percent of the total weight is propellant,

2. 6.33 percent of the total weight is structure, and

3. 1.04 percent of the total weight is payload.

Clearly, most of the weight is propellant.

Rigid body load factor

First, we assume Atlas I is a rigid body and determine the load factor from dynamic equilibrium. This analysis is shown in figure [fig18.2], where it is determined that the load factor is 1.213. However, this rigid body analysis omits information on frequencies. To determine frequencies of the structure we must consider flexible body dynamics, which is discussed beginning in the next section.

Steps in flexible body dynamics

1. Equations of motion relative to the equilibrium state.

2. Free vibration problem; normal modes.

3. Coordinate transformation from physical coordinates to modal coordinates.

4. Transform initial conditions from physical coordinates to modal coordinates. Solve for the time history of the modal coordinates.

5. Transform back to physical coordinates.

6. Calculate the payload acceleration time history.

Eigenvalue problems for real symmetric matrices

In general an nxn real matrix \([A]\) acts on an nx1 vector \(\{\phi\}\) by changing both its magnitude and direction. For example consider \[\begin{aligned} \label{18.4a} [A]=\left[\begin{array}{@{}ll@{}}2 & 1 \\1 & 2\end{array}\right] \quad\{\phi\}=\left[\begin{array}{@{}l@{}}2 \\0\end{array}\right] \quad[A]\{\phi\}=\left[\begin{array}{@{}l@{}}4 \\2\end{array}\right].\tag{a}\end{aligned}\] Vectors \(\{\phi\}\) and \([A]\{\phi\}\) are depicted in figure [fig18.5](a). Matrix \([A]\) may act on certain vectors by changing only their magnitude, and leaving their direction unchanged; e.g., \[\begin{aligned} \label{18.4b} \left[\begin{array}{@{}ll@{}}2 & 1 \\1 & 2\end{array}\right]\left[\begin{array}{@{}l@{}} 1\\1\end{array}\right]=\left[\begin{array}{@{}l@{}}3 \\3\end{array}\right]=3\left[\begin{array}{@{}l@{}}1 \\1\end{array}\right].\tag{b}\end{aligned}\] Vector \(\left[\begin{array}{@{}ll@{}}1 & 1\end{array}\right]^{T}\) is an eigenvector of matrix \([A]\) and the number “3” is an eigenvalue. Consider the vector \(\left[\begin{array}{@{}ll@{}}-1 & 1\end{array}\right]^{T}\). Then, the matrix \([A]\) operating on this second vector is \[\begin{aligned} \label{18.4c} \left[\begin{array}{@{}ll@{}}2 & 1 \\1 & 2\end{array}\right]\left[\begin{array}{@{}c@{}}-1 \\1\end{array}\right]=\left[\begin{array}{@{}c@{}}-1 \\1\end{array}\right].\tag{c}\end{aligned}\] Thus, vector \(\left[\begin{array}{@{}ll@{}}-1 & 1\end{array}\right]^{T}\) is an eigenvector and the number “1” is the eigenvalue. The eigenvectors of matrix \([A]\) are depicted in figure [fig18.5](b).

The eigenpairs of matrix \([A]\) are written as \[\begin{aligned} \label{18.4d} \left(\lambda_{1},\left\{\phi_{1}\right\}\right)=\left(3,\left[\begin{array}{@{}l@{}}1 \\1\end{array}\right]\right) \quad \text{and } \quad\left(\lambda_{2},\left\{\phi_{2}\right\}\right)=\left(1,\left[\begin{array}{@{}c@{}}-1 \\1\end{array}\right]\right).\tag{d}\end{aligned}\]

The standard matrix eigenvalue problem is \[\begin{aligned} \label{eq18.38} [A]\{\phi\}=\lambda\{\phi\},\end{aligned}\] where \([A]\) is a real nxn matrix and \(\{\phi\}\) is an nx1 vector. It is also written as \[\begin{aligned} \label{18.4e} ([A]-\lambda[I])\{\phi\}=0.\tag{e}\end{aligned}\] where \([I]\) is the nxn identity matrix. (i.e., a diagonal matrix with each diagonal element equal to one). One solution to eq. ([eq18.38]) is \(\{\phi\}=0_{n X 1}\), which is called the trivial solution. Non-trivial solutions for \(\{\phi\}\) require \[\begin{aligned} \label{eq18.39} \operatorname{\it Det}([A]-\lambda[I])=p(\lambda)=0,\end{aligned}\] where \(p(\lambda)\) is a polynomial of degree \(n\) having \(n\) roots. The polynomial \(p(\lambda)\) is called the characteristic equation. The roots of \(p(\lambda)=0\) are the eigenvalues of matrix \([A]\). The eigenvalue problem for the r-th mode is \[\begin{aligned} \label{eq18.40} [A]\left\{\phi_{r}\right\}=\lambda_{r}\left\{\phi_{r}\right\}.\end{aligned}\] The eigenpairs are denoted by \(\left(\lambda_{r},\left\{\phi_{r}\right\}\right)\). If \(\left\{\phi_{r}\right\}\) is a solution to eq. ([eq18.40])), then \(c_{r}\left\{\phi_{r}\right\}\) is also a solution where \(c_{r}\) is an arbitrary constant. Thus, the solution for the eigenvector is not unique. Scaling, or normalization, is the process to make the amplitude of the eigenvector unique. One method to normalize the amplitude is to require \(\left\{\phi_{r}\right\}\) be a unit vector in the \(n\)-dimensional space. That is, we specify \[\begin{aligned} \label{eq18.41} \left\{\phi_{r}\right\}^{T}\left\{\phi_{r}\right\}=1.\end{aligned}\] Consider the eigenvalue problem for s-th mode \[\begin{aligned} \label{eq18.42} [A]\left\{\phi_{s}\right\}=\lambda_{s}\left\{\phi_{s}\right\}.\end{aligned}\]

Pre-multiply the r-th mode in eq. ([eq18.40]) by \(\left\{\phi_{s}\right\}^{T}\) to get \[\begin{aligned} \label{eq18.43} \left\{\phi_{s}\right\}^{T}[A]\left\{\phi_{r}\right\}=\lambda_{r}\left\{\phi_{s}\right\}^{T}\left\{\phi_{r}\right\}.\end{aligned}\] Pre-multiply the s-th mode in eq. ([eq18.42]) by \(\left\{\phi_{r}\right\}^{T}\): \[\begin{aligned} \label{eq18.44} \left\{\phi_{r}\right\}^{T}[A]\left\{\phi_{s}\right\}=\lambda_{s}\left\{\phi_{r}\right\}^{T}\left\{\phi_{s}\right\}.\end{aligned}\] Take the transpose of eq. ([eq18.43]) and use the property of symmetry \([A]^{T}=[A]\) to find \[\begin{aligned} \label{eq18.45} \left\{\phi_{r}\right\}^{T}[A]\left\{\phi_{s}\right\}=\lambda_{r}\left\{\phi_{r}\right\}^{T}\left\{\phi_{s}\right\}.\end{aligned}\] Subtract eq. ([eq18.45]) from eq. ([eq18.44]) to find \[\begin{aligned} \label{eq18.46} 0=\left(\lambda_{r}-\lambda_{s}\right)\left\{\phi_{r}\right\}^{T}\left\{\phi_{s}\right\}.\end{aligned}\] For \(\lambda_{r} \neq \lambda_{s}\) \[\begin{aligned} \label{eq18.47} \left\{\phi_{r}\right\}^{T}\left\{\phi_{s}\right\}=0 \quad r \neq s \quad r, s=1,2, \ldots, n.\end{aligned}\] The eigenvectors for distinct eigenvalues are orthogonal. The eigenvectors are said to be orthonormal if \[\begin{aligned} \label{eq18.48} \left\{\phi_{r}\right\}^{T}\left\{\phi_{s}\right\}= \begin{cases}1 & r=s \\ 0 & r \neq s\end{cases}.\end{aligned}\]

Since the characteristic equation \(p(\lambda)=0\) has real coefficients, the roots \(\lambda\) are real or occur as complex conjugate pairs. It can be shown that roots are real (Boresi, 1965, p. 34).

Let \([\Phi]\) denote the nxn modal matrix whose columns are the orthonormal eigenvectors. The \(r\)-th column contains the eigenvector \(\left\{\phi_{r}\right\}\) associated with eigenvalue \(\lambda_{r}\), \(r=1,2, \ldots, n\). Since the eigenvectors are orthonormal, the modal matrix has the property \[\begin{aligned} \label{eq18.49} [\Phi]^{T}[\Phi]=[\Phi][\Phi]^{T}=[I].\end{aligned}\] The modal matrix is called an orthogonal matrix, and it has the properties \[\begin{aligned} \label{eq18.50} [\Phi]^{-1}=[\Phi]^{T} \quad \operatorname{\it Det}[\Phi]=1.\end{aligned}\] Using the modal matrix the standard eigenvalue problem can be written as \[\begin{aligned} \label{eq18.51} [A][\Phi]=[\Phi][D],\end{aligned}\] where \([D]\) is a diagonal matrix of the eigenvalues: \[\begin{aligned} \label{eq18.52} [D]=\operatorname{diag}\left[\lambda_{1}, \lambda_{2}, \ldots, \lambda_{n}\right].\end{aligned}\] Post-multiply eq. ([eq18.51]) by \([\Phi]^{T}\) to get \[\begin{aligned} \label{eq18.53} [A] \underbrace{[\Phi][\Phi]^{T}}_{[I]}=[\Phi][D][\Phi]^{T}.\end{aligned}\] Hence, matrix \([A]\) has the decomposition \[\begin{aligned} \label{eq18.54} [A]=[\Phi][D][\Phi]^{T}.\end{aligned}\] An nxn real, symmetric matrix is positive definite if all eigenvalues \(\lambda_{r}>0\), \(r=1,2, \ldots, n\). It is positive semidefinite if the eigenvalues \(\lambda_{r} \geq 0\), that is, some of its eigenvalues may equal zero. For structures having rigid body freedoms the stiffness matrix [K] has as many zero eigenvalues as it has rigid body modes.

The mode superposition method leads to the generalized eigenvalue problem \[\begin{aligned} \label{eq18.55} [K]\{\phi\}=\lambda[M]\{\phi\},\end{aligned}\] where the nxn matrices [K] and [M] are real and symmetric, \(\{\phi\}\) is the nx1 eigenvector, or modal vector, and the eigenvalue \(\lambda\) is equal to the square of the natural frequency \(\omega\). Refer to eq. ([eq18.17]). The generalized eigenvalue problem ([eq18.55]) can be reduced to the standard form in eq. ([eq18.38]) if the mass matrix is positive definite (Bathe, p. 573). Then all the properties of the eigenvalues, eigenvectors, and characteristic polynomials can be derived from the properties of the standard eigenvalue problem. In particular, for the generalized eigenvalue problem

• all \(n\) of the \(\lambda_{r}\)’s are real (no complex conjugates), and

• for \(r \neq s\), \(\left\{\phi_{r}\right\}^{T}[M]\left\{\phi_{s}\right\}=0\) (i.e., the eigenvectors are M-orthogonal).

Vanishing of the determinate \(\operatorname{\it Det}([K]-\lambda[M])=p(\lambda)=0\) is equivalent to finding the roots of a polynomial of degree \(n\). If \(n \geq 5\), iterative solution procedures have to be employed (Graig, 1981). There are many efficient subroutines available in various scientific subroutine libraries to numerically solve eigenvalue problems. In Mathematica the function “Eigensystem[K,M]” finds numerical eigenvalues and eigenvectors of the generalized eigenvalue problem. In MATLAB the statement “[V,D] = EIG(K,M)” finds a diagonal matrix D of generalized eigenvalues and a full matrix V whose columns are the corresponding eigenvectors so that \([K][V]=[M][V][D]\). For the numerical solutions to the practice exercises listed in article 1.10 it is assumed the reader has access to a subroutine to find eigenvalues and eigenvectors.

Hamilton’s principle and Lagrange’s equations

Consider a particle of mass M in motion along the \(x_1\)-\(x_2\) plane as shown in figure [fig18.6]. The position vector of the particle at time \(t\) is denoted by \(\smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{R}}}(t)\), where \[\begin{aligned} \label{eq18.56} \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{R}}}(t)=x_{1}(t) \hat{i}_{1}+x_{2}(t) \hat{i}_{2},\end{aligned}\] where \(\hat{i}_{1}\) and \(\hat{i}_{2}\) are unit vectors in the \(x_{1}\)-, and \(x_2\)-directions, respectively. From Newton’s second law the equation of motion for the particle is solved for its motion \(\smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{R}}}(t)\), \(t>t_{0}\), given its position \(\smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{R}}}\left(t_{0}\right)\) and the velocity \(\frac{d \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{R}}}}{d t}\left(t_{0}\right)\) at \(t=t_{0}\). In contrast to Newton’s method, Hamilton’s approach to this moving particle is to determine the motion \(\smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{R}}}(t)\) given the positions \(\smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{R}}}\left(t_{0}\right)\) and \(\smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{R}}}\left(t_{1}\right)\), \(t_{1}>t_{0}\). All positions in the interval \(t_{0} \leq t \leq t_{1}\) describe a time-dependent path illustrated by a curve in figure [fig18.6]. The path labeled \(\smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{R}}}(t)\) represents the actual motion. The path labeled \(\smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{R}}}^{*}(t)\) represents a motion from \(\smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{R}}}\left(t_{0}\right)\) to \(\smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{R}}}\left(t_{1}\right)\) in the vicinity of the actual motion.

Begin with D’Alembert’s principle and introduce the inertia force in addition the applied force \(\smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{F}}}\) and write the dynamic equation of equilibrium as follows: \[\begin{aligned} \label{eq18.57} \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{F}}}-M \frac{d^{2} \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{R}}}}{d t^{2}}=0.\end{aligned}\] At a fixed time introduce an infinitesimal virtual displacement denoted by the function \(\delta \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{R}}}(t)\). We can imagine the motion is stopped at time \(t\) when the virtual displacement is imposed on the particle with both the applied force and inertia force acting on the particle. The virtual changes in the coordinates imposed on the particle at a fixed time are denoted by \(\delta x_{1}(t)\) and \(\delta x_{2}(t)\), so that \[\begin{aligned} \label{eq18.58} \delta \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{R}}}(t)=\delta x_{1}(t) \hat{i}_{1}+\delta x_{2}(t) \hat{i}_{2}.\end{aligned}\] The varied path \(\smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{R}}}^{*}(t)\) in the vicinity of the actual path \(\smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{R}}}(t)\) is \[\begin{aligned} \label{eq18.59} \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{R}}}^{*}(t)=\smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{R}}}(t)+\delta \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{R}}}(t).\end{aligned}\] Since the varied path and actual path coincide at \(t = t_0\) and \(t=t_1\) the virtual displacement vanishes at these times: \[\begin{aligned} \label{eq18.60} \delta \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{R}}}\left(t_{0}\right)=\delta \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{R}}}\left(t_{1}\right)=0 \hat{i}_{1}+0 \hat{i}_{2}.\end{aligned}\] The virtual work of the applied force and the inertia force acting through the virtual displacement is given by the scalar product of eq. ([eq18.57]) with the virtual displacement in eq. ([eq18.58]): \[\begin{aligned} \label{eq18.61} \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{F}}} \bullet \delta \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{R}}}-M \frac{d^{2} \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{R}}}}{d t^{2}} \bullet \delta \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{R}}}=0.\end{aligned}\] Consider the relation \[\begin{aligned} \label{eq18.62} \frac{d}{d t}\left(\frac{d \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{R}}}}{d t} \bullet \delta \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{R}}}\right)=\frac{d^{2} \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{R}}}}{d t^{2}} \bullet \delta \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{R}}}+\frac{d \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{R}}}}{d t} \bullet \frac{d}{d t}(\delta \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{R}}}).\end{aligned}\] The virtual change in velocity is defined by \[\begin{aligned} \label{eq18.63} \delta\left(\frac{d \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{R}}}}{d t}\right) \equiv \frac{d \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{R}}}^{*}}{d t}-\frac{d \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{R}}}}{d t}=\frac{d}{d t}\left(\smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{R}}}^{*}-\smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{R}}}\right)=\frac{d}{d t}(\delta \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{R}}}).\end{aligned}\] Equation ([eq18.63]) shows we can interchange the virtual operator \(\delta(\quad)\) and the differential operator \(d(\quad)/d t\). Use the result from eq. ([eq18.63]) to write eq. ([eq18.62]) as \[\begin{aligned} \label{eq18.64} \frac{d}{d t}\left(\frac{d \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{R}}}}{d t} \bullet \delta \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{R}}}\right)=\frac{d^{2} \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{R}}}}{d t^{2}} \bullet \delta \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{R}}}+\frac{d \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{R}}}}{d t} \bullet \delta\left(\frac{d \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{R}}}}{d t}\right).\end{aligned}\] Equation ([eq18.64]) is incorporated into the virtual work eq. ([eq18.61]) to get \[\begin{aligned} \label{eq18.65} \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{F}}} \bullet \delta \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{R}}}+M \frac{d \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{R}}}}{d t} \bullet \delta\left(\frac{d \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{R}}}}{d t}\right)-\frac{d}{d t}\left(M \frac{d \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{R}}}}{d t} \bullet \delta \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{R}}}\right)=0.\end{aligned}\] The virtual work of the applied force is written as \(\delta W=\smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{F}}} \bullet \delta \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{R}}}\), and let the velocity of the particle be denoted by \(\smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{v}}}=\frac{d \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{R}}}}{d t}\), and its virtual velocity by \(\delta \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{v}}}=\delta\left(\frac{d \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{R}}}}{d t}\right)\). Then, eq. ([eq18.65]) for the virtual work is \[\begin{aligned} \label{eq18.66} \delta W+M \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{\nu}}} \bullet \delta \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{\nu}}}-\frac{d}{d t}(M \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{\nu}}} \bullet \delta \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{R}}})=0.\end{aligned}\] The kinetic energy of the particle is \(T=\frac{1}{2} M \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{\nu}}} \bullet \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{\nu}}}\), and its kinetic energy on the varied path is \[\begin{aligned} \label{eq18.67} T+\Delta T=\frac{1}{2} M(\smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{\nu}}}+\delta \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{\nu}}}) \bullet(\smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{\nu}}}+\delta \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{\nu}}})=\underbrace{\frac{1}{2} M \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{\nu}}} \bullet \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{\nu}}}}_{T}+\underbrace{M \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{\nu}}} \bullet \delta \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{\nu}}}}_{\delta T}+\underbrace{\frac{1}{2} M \delta \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{\nu}}} \bullet \delta \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{\nu}}}}_{\delta^{2} T}.\end{aligned}\] For infinitesimal virtual velocities \(\Delta T \approx \delta T\). Integrate eq. ([eq18.66]) with respect to time from \(t_0\) to \(t_1\) to get \[\begin{aligned} \label{eq18.68} \int_{t_{0}}^{t_{1}}(\delta W+\delta T) d t-(M \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{v}}} \cdot \delta \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{R}}})\Big|_{t_{0}} ^{t_{1}}=0.\end{aligned}\] The virtual displacement \(\delta \smash{\mathord{\buildrel{\lower 3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{R}}}\) vanishes at times \(t_0\) and \(t_1\). Then eq. ([eq18.68]) becomes \[\begin{aligned} \label{eq18.69} \int_{t_{0}}^{t_{1}}(\delta W+\delta T) d t=0.\end{aligned}\] Equation ([eq18.69]) is Hamilton’s principle. Note that \(\delta W\) is not the total virtual work. It is the virtual work of the non-inertial forces.

The dynamic response of an elastic thin-walled bar

Starting with Hamilton’s principle we develop a weak form for the dynamic response of a thin-walled bar suitable for a finite element analysis. Assume the bar is homogenous with a uniform cross section along the \(z\)-axis, \(0 \leq z \leq L\), where L denotes the length of the bar. Since the bar is a one-dimensional continuum, Hamilton’s principle ([eq18.69]) is expressed as \[\begin{aligned} \label{eq18.80} \int_{t_{0}}^{t_{1}}\left[\int_{0}^{L}(\delta \bar{T}+\delta \bar{W}) d z\right] d t=0,\end{aligned}\] where \(\delta \bar{T}\) is the virtual kinetic energy per unit axial length, and \(\delta \bar{W}\) is the virtual work of the non-inertial forces per unit axial length. From eqs. ([eq3.20]) and ([eq3.22]) on page , and eq. ([eq3.26]) on page , the displacements of a material point on the contour of the cross section are \[\begin{gathered} \label{eq18.81} u_{s}(s,~z,~t)=-u(z,~t) \sin \theta(s)+v(z,~t) \cos \theta(s)+r_{n}(s) \phi_{z}(z,~t), \\ \label{eq18.82} u_{\zeta}(s,~z,~t)=u(z,~t) \cos \theta(s)+v(z,~t) \sin \theta(s)-r_{t}(s) \phi_{z}(z,~t),\ \text{and} \\ \label{eq18.83} u_{z}(s,~z,~t)=w(z,~t)+y(s) \phi_{x}(z,~t)+x(s) \phi_{y}(z,~t).\end{gathered}\] The virtual kinetic energy term in Hamilton’s principle ([eq18.80]) is written as \[\begin{aligned} \label{eq18.84} \int_{t_{0}}^{t_{1}}\left(\int_{0}^{L} \delta \bar{T} d z\right) d t=\int_{t_{0}}^{t_{1}}\left[\int_{0}^{L}\left(\int_{c} \delta T t d s\right) d z\right] d t,\end{aligned}\] where the kinetic energy per unit volume, or the kinetic energy density, is \[\begin{aligned} \label{eq18.85} T=\frac{1}{2} \rho\left[\left(\dot{u}_{s}\right)^{2}+\left(\dot{u}_{\zeta}\right)^{2}+\left(\dot{u}_{z}\right)^{2}\right],\end{aligned}\] The mass density is denoted by \(\rho\). The dot over the displacement function is a shorthand notation for the partial derivative of the function with respect to time; e.g., \[\begin{aligned} \label{eq18.86} \dot{u}_{s}=\frac{\partial u_{s}}{\partial t}.\end{aligned}\] The virtual kinetic energy density at a fixed time is \[\begin{aligned} \label{eq18.87} \delta T=\frac{\partial T}{\partial \dot{u}_{s}} \delta \dot{u}_{s}+\frac{\partial T}{\partial \dot{u}_{\zeta}} \delta \dot{u}_{\zeta}+\frac{\partial T}{\partial \dot{u}_{z}} \delta \dot{u}_{z}=\rho\left(\dot{u}_{s} \delta \dot{u}_{s}+\dot{u}_{\zeta} \delta \dot{u}_{\zeta}+\dot{u}_{z} \delta \dot{u}_{z}\right),\end{aligned}\] where \(\delta \dot{u}_{s}\), \(\delta \dot{u}_{\zeta},\text{ and }\delta \dot{u}_{z}\) are the infinitesimal virtual velocity components. Interchange the order of the temporal and spatial integrals in eq. ([eq18.84]) to write the virtual kinetic energy term as \[\begin{aligned} \label{eq18.88} \int_{t_{0}}^{t_{1}}\left[\int_{0}^{L}\left(\int_{c}\delta T t d s\right) d z\right] d t=\int_{0}^{L}\left[\int_{c}\left(\int_{t_{0}}^{t_{1}} \delta T d t\right) t d s\right] d z.\end{aligned}\] The interior integral in eq. ([eq18.88]) with respect to time is integrated by parts as follows: \[\label{eq18.89} \begin{split}\int_{t_{0}}^{t_{1}} \delta T d t &=\rho \int_{t_{0}}^{t_{1}}\left(\dot{u}_{s} \delta \dot{u}_{s}+\dot{u}_{\zeta} \delta \dot{u}_{\zeta}+\dot{u}_{z} \delta \dot{u}_{z}\right) d t \\&=-\rho \int_{t_{1}}^{t_1}\left(\ddot{u}_{s} \delta u_{s}+\ddot{u}_{\zeta} \delta u_{\zeta}+\ddot{u}_{s} \delta u_{s}\right) d t+\left.\rho\left(\dot{u}_{s} \delta u_{s}+\dot{u}_{\zeta} \delta u_{\zeta}+\dot{u}_{z} \delta u_{z}\right)\right|_{t_{0}} ^{t_{1}}\end{split}.\] By Hamilton’s principle the virtual displacements \(\left\{\delta u_{s}, \delta u_{\zeta}, \delta u_{z}\right\}\) vanish at times \(t_{0}\) and \(t_{1}\). Return to eq. ([eq18.84]) and identify \[\begin{aligned} \label{eq18.90} \int_{t_{0}}^{t_{1}}\left[\int_{0}^{t}\left(\int_c \delta T t d s\right) d z\right] d t=\int_{t_{0}}^{t_{1}}\left[\int_{0}^{L} \delta \bar{T} d z\right] d t,\end{aligned}\] where the virtual kinetic energy per unit axial length is \[\begin{aligned} \label{eq18.91} \delta \bar{T}=\int_{c} \delta T t d s=-\rho \int_{c}\left(\ddot{u}_{s} \delta u_{s}+\ddot{u}_{\zeta} \delta u_{\zeta}+\ddot{u}_{s} \delta u_{s}\right) t d s.\end{aligned}\] We assume the density of mass \(\rho\) is independent of the contour coordinate \(s\). The individual terms in the integral for \(\delta \bar{T}\) are as follows: \[\begin{gathered} \label{eq18.92} \rho \int_{c}\left(\ddot{u}_{s} \delta u_{s}\right) t d s=\rho \int_{c}\left[-\ddot{u} \sin \theta+\ddot{v} \cos \theta+r_{n} \ddot{\phi}_{z}\right]\left[-\delta u \sin \theta+\delta v \cos \theta+r_{n} \delta \phi_{z}\right] t d s, \\\label{eq18.93} \rho \int_{c}\left(\ddot{u}_{\zeta} \delta u_{\zeta}\right) t d s=\rho \int_{c}\left[\ddot{u} \cos \theta+\ddot{v} \sin \theta-r_{t} \ddot{\phi}_{z}\right]\left[\delta u \cos \theta+\delta v \sin \theta-r_{t} \delta \phi_{z}\right] t d s,\ \text{and} \\\label{eq18.94} \rho \int_{c}\left(\ddot{u}_{z} \delta u_{z}\right) t d s=\rho \int_{c}\left[\ddot{w}+y \ddot{\phi}_{x}+x \ddot{\phi}_{y}\right]\left[\delta w+y \delta \phi_{x}+x \delta \phi_{y}\right] t d s.\end{gathered}\] Sum eqs. ([eq18.92]) to ([eq18.94]) and collect terms in the integral multiplying the virtual displacements and virtual rotations. In the process of evaluating and simplifying the result for \(\delta \bar{T}\), the following terms are identified or defined:

• From eq. ([eq3.9]) on page the relationship between normal and tangential coordinates of the contour and the Cartesian coordinates \(\{x(s), y(s)\}\) of the contour is \[\begin{aligned} \label{eq18.95} r_{n} \cos \theta-r_{t} \sin \theta=x(s)-x_{s c} \quad r_{n} \sin \theta+r_{t} \cos \theta=y(s)-y_{s c}.\end{aligned}\]

• From eqs. ([eq4.2]) and ([eq4.3]) on page the geometric properties of the cross section are identified as \[\begin{aligned} \label{eq18.96} A=\int_{c} t d s \quad \int_{c} x t d s=0 \quad \int_{c} y t d s=0 \quad\left(I_{x x}, I_{y y}, I_{x y}\right)=\int_{c}\left(y^{2}, x^{2}, x y\right) t d s.\end{aligned}\]

• Terms associated with the mass are defined by \[\begin{gathered} \label{eq18.97} m=\rho \int_{c} t d s=\rho A, \\ \rho\left(I_{x x}, I_{y y}, I_{x y}\right)=\rho A\left(I_{x x}/A, I_{y y}/A, I_{x y}/A\right)=m\left(r_{x}^{2}, r_{y}^{2}, r_{x y}^{2}\right),\ \text{and}\label{eq18.98}\\ I_{z}=\rho \int_{c}\left(r_{n}^{2}+r_{t}^{2}\right) t d s=\rho A\left[\frac{1}{A} \int_{c}\left(r_{n}^{2}+r_{t}^{2}\right) t d s\right] = m r_{z}^{2}.\label{eq18.99}\end{gathered}\]

The mass per unit length is denoted by \(m\), rotary inertial moments due to bending by \(\rho(I_{x x}, I_{y y}, I_{x y})\), and the polar moment of inertia about the \(z\)-axis through the shear center by \(I_{z}\). The final result for virtual kinetic energy per unit length in eq. ([eq18.91]) is \[\begin{aligned} \label{eq18.100} \delta \bar{T} =&-\left\{m\left(\ddot{u}+y_{s c} \ddot{\phi}_{z}\right) \delta u+m\left(\ddot{v}-x_{s c} \ddot{\phi}_{z}\right) \delta v+m \ddot{w} \delta w+m\left(r_{x}^{2} \ddot{\phi}_{x}+r_{x y}^{2} \ddot{\phi}_{y}\right) \delta \phi_{x}+m\left(r_{x y}^{2} \ddot{\phi}_{x}+r_{y}^{2} \ddot{\phi}_{y}\right) \delta \phi_{y}\right.\nonumber\\ &+\left.m\left(r_{z}^{2} \ddot{\phi}_{z}+y_{s c} \ddot{u}-x_{s c} \ddot{v}\right) \delta \phi_{z}\right\}.\end{aligned}\]

The virtual work per unit axial length is determined from conservative forces and non-conservative forces. It is expressed as \[\begin{aligned} \label{eq18.101} \delta \bar{W}=-\delta \bar{U}+\delta \bar{W}_{\text {n.c.}},\end{aligned}\] where \(\bar{U}\) is the potential energy of the conservative forces and \(\delta \bar{W}_{\text {n.c. }}\) is the virtual work of the non-conservative forces acting at point \(z\) in the bar. The potential energy per unit length represents the energy due to elastic deformation of the bar and it is given by \[\begin{aligned} \label{eq18.102} \bar{U}=\bar{U}_{\varepsilon}+\bar{U}_{\gamma}.\end{aligned}\] The strain energy due extension and bending is denoted by \(\bar{U}_{\varepsilon}\), and the strain energy due to transverse shear and torsion by \(\bar{U}_{\gamma}\). From eq. ([eq5.47]) on page \[\begin{aligned} \label{eq18.103} \bar{U}_{\varepsilon}=\frac{1}{2}\left[E A\left(w^{\prime}\right)^{2}+E I_{x x}\left(\phi_{x}^{\prime}\right)^{2}+E I_{y y}\left(\phi_{y}^{\prime}\right)^{2}+2 E I_{x y}\left(\phi_{x}^{\prime}\right)\left(\phi_{y}^{\prime}\right)\right],\end{aligned}\] and from eq. ([eq5.82]) on page \[\begin{aligned} \label{eq18.104} \bar{U}_{\gamma}=\frac{1}{2}\left[s_{x x} \psi_{x}^{2}+2 s_{x y} \psi_{x} \psi_{y}+s_{y y} \psi_{y}^{2}+G J\left(\phi_{z}^{\prime}\right)^{2}\right].\end{aligned}\] Note that the partial derivatives of the axial displacement and bending rotations are denoted by a prime superscript in eq. ([eq18.103]). That is, \[\begin{aligned} \label{eq18.105} w^{\prime}=\frac{\partial w}{\partial z} \quad \phi_{x}^{\prime}=\frac{\partial \phi_{x}}{\partial z} \quad \phi_{y}^{\prime}=\frac{\partial \phi_{y}}{\partial z}.\end{aligned}\] The transverse shear strains in eq. ([eq18.104]) are related to the displacements and rotations by \[\begin{aligned} \label{eq18.106} \psi_{x}=u^{\prime}+\phi_{y} \quad \Psi_{y}=v^{\prime}+\phi_{x}.\end{aligned}\] Refer to figure [fig3.6] on page . The virtual strains are denoted by \(\delta w^{\prime}\), \(\delta \phi_{x}^{\prime}\), \(\delta \phi_{y}^{\prime}\), \(\delta \psi_{x}\), \(\delta \Psi_{y},\text{ and }\delta \phi_{z}\), where \[\begin{aligned} \label{eq18.107} \delta \psi_{x}=\delta u^{\prime}+\delta \phi_{y},\text{ and }\delta \psi_{y}=\delta v^{\prime}+\delta \phi_{x}.\end{aligned}\] At a fixed value of the coordinate \(z\) the virtual strain energies per unit axial length are \[\begin{aligned} \label{eq18.108} \delta \bar{U}_{\varepsilon}&=\frac{\partial \bar{U}_{\varepsilon}}{\partial w^{\prime}} \delta w^{\prime}+\frac{\partial \bar{U}_{\varepsilon}}{\partial \phi_{x}^{\prime}} \delta \phi_{x}^{\prime}+\frac{\partial \bar{U}_{\varepsilon}}{\partial \phi_{y}^{\prime}} \delta \phi_{y}^{\prime}=N \delta w^{\prime}+M_{x} \delta \phi_{x}^{\prime}+M_{y} \delta \phi_{y}^{\prime},\ \text{and} \\\label{eq18.109} \delta \bar{U}_{\gamma} &=\frac{\partial \bar{U}_{\gamma}}{\partial \psi_{x}} \delta \psi_{x}+\frac{\partial \bar{U}_{\gamma}}{\partial \psi_{y}} \delta \psi_{y}+\frac{\partial \bar{U}_{\gamma}}{\partial \phi_{z}^{\prime}} \delta \phi_{z}^{\prime}=V_{x} \delta \psi_{x}+V_{y} \delta \psi_{y}+M_{z} \delta \phi_{z}^{\prime}.\end{aligned}\] The generalized forces corresponding to the virtual strains in eqs. ([eq18.108]) and ([eq18.109]) are given by \[\begin{aligned} \label{eq18.110} N=E A w^{\prime} \quad M_{x}=E I_{x x} \phi_{x}^{\prime}+E I_{x y} \phi_{y}^{\prime} \quad M_{y}=E I_{x y} \phi_{x}^{\prime}+E I_{y y} \phi_{y}^{\prime},\ \text{and}\\ \label{eq18.111} V_{x}=s_{x x} \psi_{x}+s_{x y} \psi_{y} \quad V_{y}=s_{x y} \psi_{x}+s_{y y} \psi_{y} \quad M_{z}=G J \phi_{z}^{\prime}.\end{aligned}\] The total virtual strain energy is the sum of eqs. ([eq18.108]) and ([eq18.109]). That is, \[\begin{aligned} \label{eq18.112} \delta \bar{U}=\delta \bar{U}_{\varepsilon}+\delta \bar{U}_{\gamma}=N \delta w^{\prime}+M_{x} \delta \phi_{x}^{\prime}+M_{y} \delta \phi_{y}^{\prime}+V_{x} \delta \psi_{x}+V_{y} \delta \psi_{y}+M_{z} \delta \phi_{z}^{\prime}.\end{aligned}\] Substitute eq. ([eq18.101]) into Hamilton’s principle ([eq18.80]) to get \[\begin{aligned} \label{eq18.113} \int_{t_{0}}^{t_{1}}\left[\int_{0}^{L}(\delta \bar{T}+\delta \bar{W}) d z\right] d t=-\int_{t_0}^{t_1}\left[\int_{0}^{L}\left(-\delta \bar{T}+\delta \bar{U}-\delta \bar{W}_{\mathrm{n}, \mathrm{c}.}^{L}\right) d z\right] d t=0.\end{aligned}\] Equation ([eq18.113]) is satisfied for \(t_{0} \leq t \leq t_{1}\) by \[\begin{aligned} \label{eq18.114} \left[\int_{0}^{L}\left(-\delta \bar{T}+\delta \bar{U}-\delta \bar{W}_{\mathrm{n}, \mathrm{c}.}\right) d z\right]=0.\end{aligned}\] Substitute eq. ([eq18.100]) for \(\delta \bar{T}\) and eq. ([eq18.112]) for \(\delta \bar{U}\) into eq. ([eq18.114]) to get the weak form for the elastodynamics of the thin-walled bar as \[\begin{aligned} \label{eq18.115} \begin{split} &\int_{0}^{L}\left\{\left[\left(m\left(\ddot{u}+y_{s c} \ddot{\phi}_{z}\right) \delta u+m\left(\ddot{v}-x_{s c} \ddot{\phi}_{z}\right) \delta v+m \ddot{w} \delta w+m\left(r_{x}^{2} \ddot{\phi}_{x}+r_{x y}^{2} \ddot{\phi}_{y}\right) \delta \phi_{x}+m\left(r_{x y}^{2} \ddot{\phi}_{x}+r_{y}^{2} \ddot{\phi}_{y}\right) \delta \phi_{y}+\right.\right.\right. \\ &\left.\left.\quad m\left(r_{z}^{2} \ddot{\phi}_{z}+y_{s c} \ddot{u}-x_{s c} \ddot{v}\right) \delta \phi_{z}\right]+\left[N \delta w^{\prime}+M_{x} \delta \phi_{x}^{\prime}+M_{y} \delta \phi_{y}^{\prime}+V_{x} \delta \psi_{x}+V_{y} \delta \psi_{y}+M_{z} \delta \phi_{z}^{\prime}\right]\right\} d z=\int_{0}^{L} \delta \bar{W}_{\text {n.c. }} d z\end{split}.\end{aligned}\]

-1The non-conservative virtual work is expressed in terms of the prescribed external loads acting on the bar shown in figure [fig18.8]. The generalized forces acting on the cross sections at \(z=0\) and \(z=L\) are denoted by \(Q_{i}(t)\), \(i=1,2, \ldots, 12\). The prescribed distribution of force intensities (F/L) and distribution of moments intensities(F-L/L) are functions of coordinate \(z\) and time \(t\). Distributions \(f_{x}(z,~t)\) and \(f_{y}(z,~t)\) are resolved at the shear center, and distribution \(f_{z}(z,~t)\) is resolved at the centroid. Distributions of moments \(m_{x}(z,~t)\) and \(m_{y}(z,~t)\) are resolved at the centroid, and distribution \(m_{z}(z,~t)\) is resolved at the shear center. Also, refer to figure [fig3.8](b) on page .

The virtual displacements and rotations are independent variables. We can take \(\delta u\) and its derivative \(\delta u'\) not equal to zero, and the remaining virtual displacements and rotations equal to zero in eq. ([eq18.115]). Then the weak form governing virtual displacement \(\delta u\) is \[\begin{aligned} \label{eq18.116} \int_{0}^{L}\left[\left(m\left(\ddot{u}+y_{s c} \ddot{\phi}_{z}\right) \delta u+V_{x} \delta u^{\prime}\right] d z=\int_{0}^{L} f_{x}(z, t) \delta u(z, t) d z+Q_{1}(t) \delta u(0, t)+Q_{7}(t) \delta u(L, t)\right.\!.\end{aligned}\] From eq. ([eq18.115]) and eq. ([eq18.107]) the weak form governing \(\delta v\) and its derivative is \[\begin{aligned} \label{eq18.117} \int_{0}^{L}\left[m\left(\ddot{v}-x_{s c} \ddot{\phi}_{z}\right) \delta v+V_{y} \delta v^{\prime}\right] d z=\int_{0}^{L} f_{y}(z, t) \delta v(z, t) d z+Q_{2}(t) \delta v(0, t)+Q_{8}(t) \delta v(L, t).\end{aligned}\] The weak form governing \(\delta w\) and its derivative is \[\begin{aligned} \label{eq18.118} \int_{0}^{L}\left[m \ddot{w} \delta w+N \delta w^{\prime}\right] d z=\int_{0}^{L} f_{z}(z, t) \delta w(z, t) d z+Q_{3}(t) \delta w(0, t)+Q_{9}(t) \delta w(L, t).\end{aligned}\] The remaining weak forms governing the virtual rotations in the order \(\delta \phi_{x}\), \(\delta \phi_{y},\text{ and }\delta \phi_{z}\) are as follows: \[\begin{aligned} &\int_{0}^{L}\left[m\left(r_{x}^{2} \ddot{\phi}_{x}+r_{x y}^{2} \ddot{\phi}_{y}\right) \delta \phi_{x}+M_{x} \delta \phi_{x}^{\prime}+V_{y} \delta \phi_{x}\right] d z\nonumber\\ &\qquad=\int_{0}^{L} m_{x}(z, t) \delta \phi_{x}(z, t) d z+Q_{4}(t) \delta \phi_{x}(0, t)+Q_{10}(t) \delta \phi_{x}(L, t),\label{eq18.119}\\ &\int_{0}^{L}\left[m\left(r_{x y}^{2} \ddot{\phi}_{x}+r_{y}^{2} \ddot{\phi}_{y}\right) \delta \phi_{y}+M_{y} \delta \phi_{y}^{\prime}+V_{x} \delta \phi_{y}\right] d z\nonumber\\* &\qquad=\int_{0}^{L} m_{y}(z, t) \delta \phi_{y}(z, t) d z+Q_{5}(t) \delta \phi_{y}(0, t)+Q_{11}(t) \delta \phi_{y}(L, t),\ \text{and}\label{eq18.120}\\ &\int_{0}^{L}\left[m\left(r_{z}^{2} \ddot{\phi}_{z}+y_{s c} \ddot{u}-x_{s c} \ddot{v}\right) \delta \phi_{z}+M_{z} \delta \phi_{z}^{\prime}\right] d z\nonumber\\ &\qquad=\int_{0}^{L} m_{z}(z, t) \delta \phi_{z}(z, t) d z+Q_{6}(t) \delta \phi_{z}(0, t)+Q_{12}(t) \delta \phi_{z}(L, t).\label{eq18.121}\end{aligned}\]

Finite element model for the dynamics of an axial bar

The weak form governing the response of the bar is given by eq. ([eq18.118]). Following the discussion in article [sec17.2] on page consider the mesh \[\begin{aligned} \label{eq18.122} 0=z_{1}<z_{2}<z_{3}<\ldots<z_{M}<z_{M+1}=L,\end{aligned}\] where M denotes the number of elements, and M+1 is the number of nodes. The \(k\)th element is denoted by \[\begin{aligned} \label{eq18.123} \Omega_{k} \equiv\left\{z \mid z_{k} \leq z \leq z_{k+1}\right\} \quad k=1,2, \ldots, M.\end{aligned}\] Each element is mapped into a standard element denoted by \[\begin{aligned} \label{eq18.124} \Omega_{\mathrm{st}} \equiv\{\zeta \mid-1<\zeta<1\}.\end{aligned}\] The standard element is mapped onto the \(k\)th element by \[\begin{aligned} \label{eq18.125} z=\eta_{1}(\zeta) z_{k}+\eta_{2}(\zeta) z_{k+1}, \text{where}\ \eta_{1}(\zeta)=(1-\zeta)/2\ \text{and}\ \eta_{2}(\zeta)=(1+\zeta)/2.\end{aligned}\] The inverse mapping is \[\begin{aligned} \label{eq18.126} \zeta=\frac{2 z-z_{z}-z_{k+1}}{z_{k+1}-z_{k}}=\frac{2 z-z_{z}-z_{k+1}}{h_{k}} \quad z \in \Omega_{k}.\end{aligned}\] The length of the element is denoted by \(h_{k}\) where \(h_{k}=z_{k+1}-z_{k}\). Let \(w^{(k)}(z, t)\) denote the axial displacement in \(\Omega_k\), which expressed in terms of the nodal displacements and the shape functions is \[\begin{aligned} \label{eq18.127} w^{(k)}(\zeta, t)=\eta_{1}(\zeta) w_{k}(t)+\eta_{2}(\zeta) w_{k+1}(t) \quad \zeta \in \Omega_{\mathrm{st}}.\end{aligned}\]\[\begin{aligned} \label{eq18.128} \delta w^{(k)}=b_{k} \eta_{1}(\zeta)+b_{k+1} \eta_{2}(\zeta)\end{aligned}\] The external loads acting on element \(\Omega_k\) are shown in figure [fig18.9].

Let \(\eta_{1}^{\prime}=\frac{d \eta_{1}}{d \zeta}\). The finite element representation of eq. ([eq18.118]) is \[\begin{aligned} \label{eq18.129} &\sum_{k=1}^{M} b_{k}\left\{\int_{-1}^1\left[m\left(\eta_{1} \ddot{w}_{k}+\eta_{2} \ddot{w}_{k+1}\right) \eta_{1}+E A \frac{2}{h_{k}}\left(\eta_{1}^{\prime} w_{k}+\eta_{2}^{\prime} w_{k+1}\right) \frac{2}{h_{k}} \eta_{1}^{\prime}\right] \frac{h_{k}}{2} d \zeta\right\} \nonumber\\ &\quad\qquad+ b_{k+1}\left\{\int_{-1}^{1}\left[m\left(\eta_{1} \ddot{w}_{k}+\eta_{2} \ddot{w}_{k+1}\right) \eta_{2}+E A \frac{2}{h_{k}}\left(\eta_{1}^{\prime} w_{k}+\eta_{2}^{\prime} w_{k+1}\right) \frac{2}{h_{k}} \eta_{2}^{\prime}\right] \frac{h_{k}}{2} d \zeta\right\}\nonumber\\ &\qquad=\sum_{k=1}^{M} b_{k}\left[Q_{k}(t)+F_{k}(t)\right]+b_{k+1}\left[Q_{k+1}(t)+F_{k+1}(t)\right] .\end{aligned}\]

External loads at a common node between elements are defined by \[\begin{aligned} Q_{k}(t)=Q_{9}^{(k-1)}(t)+Q_{3}^{(k)}(t) \quad k=2,3, \ldots, M-1,\label{eq18.130}\end{aligned}\] which are depicted in figure [fig18.10].

At the first node \({z_1}\) and the last node \({z_{M+1}}\) the external forces are \[\begin{aligned} \label{eq18.131} Q_{1}(t)=Q_{3}^{(1)}(t),\text{ and }Q_{M}(t)=Q_{9}^{(M)}.\end{aligned}\] The distributed loads lead to forces \[\begin{aligned} \label{eq18.132} F_{k}(t)=\int_{-1}^{1} f_{z}\left[\eta_{1} z_{k}+\eta_{2} z_{k+1}, t\right] \eta_{1} \frac{h_{k}}{2} d \zeta \quad F_{k+1}(t)=\int_{-1}^{1} f_{z}\left[\eta_{1} z_{k}+\eta_{2} z_{k+1}, t\right] \eta_{2} \frac{h_{k}}{2} d \zeta.\end{aligned}\] For \(b_{k} \neq 0,\text{ and }b_{k+1}=0\), eq. ([eq18.129]) leads to \[\begin{aligned} \label{eq18.133} m h_{k}\left(\frac{1}{3} \ddot{w}_{k}+\frac{1}{6} \ddot{w}_{k+1}\right)+\frac{E A}{h_{k}}\left(w_{k}-w_{k+1}\right)=Q_{k}(t)+F_{k}(t).\end{aligned}\] For \(b_{k}=0\) and \(b_{k+1} \neq 0\) eq. ([eq18.129]) leads to \[m h_{k}\left(\frac{1}{6} \ddot{w}_{k}+\frac{1}{3} \ddot{w}_{k+1}\right)+\frac{E A}{h_{k}}\left(-w_{k}+w_{k+1}\right)=Q_{k+1}(t)+F_{k+1}(t).\] Finally we write the matrix form for element \(\Omega_k\) as \[\begin{aligned} \label{eq18.134} \frac{m h_{k}}{6}\left[\begin{array}{@{}ll@{}}2 & 1 \\1 & 2\end{array}\right]\left[\begin{array}{@{}c@{}}\ddot{w}_{k} \\\ddot{w}_{k+1}\end{array}\right]+\frac{E A}{h_{k}}\left[\begin{array}{@{}cc@{}}1 & -1 \\-1 & 1\end{array}\right]\left[\begin{array}{@{}c@{}}w_{k} \\w_{k+1}\end{array}\right]=\left[\begin{array}{@{}c@{}}Q_{k}(t)+F_{k}(t) \\Q_{k+1}(t)+F_{k+1}(t)\end{array}\right]\end{aligned}\] Note that the mass matrix in eq. ([eq18.134]) is not diagonal, but it is symmetric. This mass matrix is called a consistent mass matrix, since the spatial distribution of the velocity function is consistent with the assumed spatial distribution of the axial displacement function ([eq18.127]).

[ex18.2]Consider a three-degree-of freedom model of a uniform bar shown in figure [fig18.11] that undergoes unrestrained axial motion. The mass per unit length \(m\) and the axial stiffness \(E A\) are constant. The bar is modeled with two finite elements \(\Omega_1\) and \(\Omega_2\), each with lengths \(h_{1}=h_{2}=L\).

For element \(\Omega_1\) the mass and stiffness matrices are \[\begin{aligned} \label{ex18.2a} M^{(1)}=\frac{m L}{6} \begin{array}{@{}l@{}} \hspace*{3pt}{\arraycolsep=2pt\begin{array}{@{}cc@{}} \ddot{w}_{1} & \ddot{w}_{2} \end{array}}\\ \left[\begin{array}{@{}ll@{}}2 & 1 \\1 & 2\end{array}\right] \end{array} \quad \text{and} \quad K^{(1)}=\frac{E A}{L} \begin{array}{@{}l@{}} \hspace*{5pt}{\arraycolsep=6pt\begin{array}{@{}cc@{}} w_{1} & w_{2} \end{array}}\\ \left[\begin{array}{@{}cc@{}}1 & -1 \\-1 & 1\end{array}\right] \end{array}.\end{aligned}\] For element \(\Omega_2\) the mass and stiffness matrices are \[\begin{aligned} \label{ex18.2b} M^{(2)}=\frac{m L}{6} \begin{array}{@{}l@{}} \hspace*{3pt}{\arraycolsep=2pt\begin{array}{@{}cc@{}} \ddot{w}_{2} & \ddot{w}_{3} \end{array}}\\ \left[\begin{array}{@{}ll@{}}2 & 1 \\1 & 2\end{array}\right]\end{array} \quad \text{and} \quad K^{(2)}=\frac{E A}{L} \begin{array}{@{}l@{}} \hspace*{6pt}{\arraycolsep=5pt\begin{array}{@{}cc@{}} w_{2} & w_{3} \end{array}}\\ \left[\begin{array}{@{}cc@{}}1 & -1 \\ -1 & 1\end{array}\right] \end{array}.\end{aligned}\] The matrix equation of motion after assembly is \[\begin{aligned} \label{ex18.2c} \frac{m L}{6}\left[\begin{array}{@{}ccc@{}}2 & 1 & 0 \\1 & 4 & 1 \\0 & 1 & 2\end{array}\right]\left[\begin{array}{@{}c@{}}\ddot{w}_{1} \\\ddot{w}_{2} \\\ddot{w}_{3}\end{array}\right]+\frac{E A}{L}\left[\begin{array}{@{}ccc@{}}1 & -1 & 0 \\-1 & 2 & -1 \\0 & -1 & 1\end{array}\right]\left[\begin{array}{@{}c@{}}w_{1} \\w_{2} \\w_{3}\end{array}\right]=\left[\begin{array}{@{}c@{}}Q_{1}(t) \\0 \\0\end{array}\right].\end{aligned}\] For the bar at rest at time \(t\) = 0 the initial axial displacements are \[\begin{aligned} \label{ex18.2d} w_{1}(0)=w_{2}(0)=w_{3}(0)=0.\end{aligned}\] Consider the impulsive response of the bar \(Q_{1}(t)=\hat{Q} \delta(t)\), where \(\hat{Q}\) is the magnitude of the impulse and \(\delta(t)\) is the Dirac delta function. The Dirac delta function is defined by \[\begin{aligned} \label{ex18.2e} \delta(t)=0\ \text{for}\ t \neq 0,\text{ and }\int_{-\infty}^{\infty} \delta(t) d t=1.\end{aligned}\]

r63pt

The Dirac delta function is depicted in figure [fig18.12], in which we take the limit as \(\varepsilon \rightarrow 0\). The amplitude of \(\delta(0)\) is undefined, but the area under the function is well defined. Note the dimensional units of the Dirac delta function are reciprocal seconds, and the dimensional units of impulse \(\hat{Q}\) are force-seconds. The impulsive force causes an initial velocity of the bar. To determine the initial velocity of the bar we integrate the equations of motion (c) with respect to time over an infinitesimal interval at \(t\) = 0 from \(t=-0\) to \(t=+0\). During this time interval there is no displacement of the bar and the velocities \(\dot{w}_{i}(-0)=0\), \(i=1,2,3\). Integration of eq. (c) over the infinitesimal time interval is \[\begin{aligned} \label{ex18.2f} \frac{m L}{6}\left[\begin{array}{lll}2 & 1 & 0 \\1 & 4 & 1 \\0 & 1 & 2\end{array}\right]\left[\begin{array}{@{}c@{}}\dot{w}_{1}(+0) \\\dot{w}_{2}(+0) \\\dot{w}_{2}(+0)\end{array}\right]=\hat{Q} \int_{-0}^{+0}\left[\begin{array}{@{}c@{}}\delta(t) \\0 \\0\end{array}\right] d t=\hat{Q}\left[\begin{array}{@{}l@{}}1 \\0 \\0\end{array}\right].\end{aligned}\] Solve eq. (f) for the initial velocities to get \[\begin{aligned} \label{ex18.2g} \left[\begin{array}{@{}c@{}}\dot{w}_{1}(+0) \\\dot{w}_{2}(+0) \\\dot{w}_{2}(+0)\end{array}\right]=\frac{\hat{Q}}{m L}\left[\begin{array}{@{}c@{}}7/2 \\-1 \\1/2\end{array}\right].\end{aligned}\]

Dynamic bending of a bar with two axes of symmetry

If the cross section is symmetric with respect to both the \(x\)- and \(y\)-axes through the centroid, then \(x_{s c}=y_{s c}=0\), \(I_{x y}=0\), \(r_{x y}=0,\text{ and }s_{x y}=0\). In this case of double symmetry, transverse bending is decoupled from torsion in both the inertia and stiffness terms. That is, the inertia axis and elastic axis coincide. However, the motions of the lateral displacement \(v(z, t)\) and the rotation \(\phi_{x}(z, t)\) are linked because of the presence of transverse shear deformation \(\psi_{y}\). The governing weak forms ([eq18.117]) and ([eq18.119]) are \[\begin{aligned} &\int_{0}^{L}\left[m(\ddot{v}) \delta v+s_{y y}\left(v^{\prime}+\phi_{x}\right) \delta v^{\prime}\right] d z=\int_{0}^{L} f_{y}(z, t) \delta v(z) d z+Q_{2}(t) \delta v(0, t)+Q_{8}(t) \delta v(L, t).\\ &\int_{0}^{L}\left[m\left(r_{x}^{2} \ddot{\phi}_{x}\right) \delta \phi_{x}+E I_{x x} \phi_{x}^{\prime} \delta \phi_{x}^{\prime}+s_{y y}\left(v^{\prime}+\phi_{x}\right) \delta \phi_{x}\right] d z\\ &\quad=\int_{0}^{L} m_{x}(z, t) \delta \phi_{x}(z, t) d z+Q_{4}(t) \delta \phi_{x}(0, t)+Q_{10}(t) \delta \phi_{x}(L, t)\end{aligned}\] The previous two equations are combined to the matrix form \[\begin{aligned} \label{eq18.137} \int_{0}^{L}\left\{\left[\delta \phi_{x} \delta v\right]\left[\begin{array}{@{}cc@{}}m r_{x}^{2} & 0 \\0 & m\end{array}\right]\left[\begin{array}{@{}c@{}}\ddot{\phi}_x \\\ddot{v}\end{array}\right]+\delta \phi_{x}^{\prime} E I_{x x} \phi_{x}^{\prime}+\left(\delta v^{\prime}+\delta \phi_{x}\right) s_{y y}\left(v^{\prime}+\phi_{x}\right)\right\} d z=\delta W_{\mathrm{n}. c},\end{aligned}\] where the virtual work of the non-conservative forces is \[\begin{aligned} \label{eq18.138} \delta W_{\text {n.c.}}=\int_{0}^{L}\left[\delta \phi_{x}~\delta v\right]\left[\begin{array}{@{}c@{}}m_{z}(z, t) \\f_{y}(z, t)\end{array}\right] d z+\left[\delta \phi_{x}(0, t) \delta v(0, t)\right]\left[\begin{array}{@{}l@{}}Q_{2}(t) \\Q_{4}(t)\end{array}\right]+\left[\begin{array}{@{}l@{}}\delta \phi_{x}(L, t)~\delta v(L, t)\end{array}\right]\left[\begin{array}{@{}c@{}}Q_{8}(t) \\Q_{10}(t)\end{array}\right].\end{aligned}\]

Callstack:
    at (Under_Construction/Aerospace_Structures_(Johnson)/18:_Introduction_to_flexible_body_dynamics), /content/body/div/div[8]/div[2]/p[2]/span[4]/span[1], line 1, column 2

Callstack:
    at (Under_Construction/Aerospace_Structures_(Johnson)/18:_Introduction_to_flexible_body_dynamics), /content/body/div/div[8]/div[2]/p[2]/span[4]/span[2], line 1, column 2

Callstack:
    at (Under_Construction/Aerospace_Structures_(Johnson)/18:_Introduction_to_flexible_body_dynamics), /content/body/div/div[8]/div[2]/p[2]/span[4]/span[3], line 1, column 2

Callstack:
    at (Under_Construction/Aerospace_Structures_(Johnson)/18:_Introduction_to_flexible_body_dynamics), /content/body/div/div[8]/div[2]/p[2]/span[4]/span[4], line 1, column 2

Callstack:
    at (Under_Construction/Aerospace_Structures_(Johnson)/18:_Introduction_to_flexible_body_dynamics), /content/body/div/div[8]/div[2]/p[2]/span[4]/span[5], line 1, column 3

Callstack:
    at (Under_Construction/Aerospace_Structures_(Johnson)/18:_Introduction_to_flexible_body_dynamics), /content/body/div/div[8]/div[2]/p[2]/span[4]/span[6], line 1, column 2

Callstack:
    at (Under_Construction/Aerospace_Structures_(Johnson)/18:_Introduction_to_flexible_body_dynamics), /content/body/div/div[8]/div[2]/p[2]/span[4]/span[7], line 1, column 2

Callstack:
    at (Under_Construction/Aerospace_Structures_(Johnson)/18:_Introduction_to_flexible_body_dynamics), /content/body/div/div[8]/div[2]/p[2]/span[4]/span[8], line 1, column 2

Callstack:
    at (Under_Construction/Aerospace_Structures_(Johnson)/18:_Introduction_to_flexible_body_dynamics), /content/body/div/div[8]/div[2]/p[2]/span[4]/span[9], line 1, column 2

Callstack:
    at (Under_Construction/Aerospace_Structures_(Johnson)/18:_Introduction_to_flexible_body_dynamics), /content/body/div/div[8]/div[2]/p[2]/span[4]/span[10], line 1, column 2

Callstack:
    at (Under_Construction/Aerospace_Structures_(Johnson)/18:_Introduction_to_flexible_body_dynamics), /content/body/div/div[8]/div[2]/p[2]/span[4]/span[11], line 1, column 2

Callstack:
    at (Under_Construction/Aerospace_Structures_(Johnson)/18:_Introduction_to_flexible_body_dynamics), /content/body/div/div[8]/div[2]/p[2]/span[4]/span[12], line 1, column 2

Callstack:
    at (Under_Construction/Aerospace_Structures_(Johnson)/18:_Introduction_to_flexible_body_dynamics), /content/body/div/div[8]/div[2]/p[2]/span[4]/span[13], line 1, column 2

Callstack:
    at (Under_Construction/Aerospace_Structures_(Johnson)/18:_Introduction_to_flexible_body_dynamics), /content/body/div/div[8]/div[2]/p[2]/span[4]/span[14], line 1, column 2

Callstack:
    at (Under_Construction/Aerospace_Structures_(Johnson)/18:_Introduction_to_flexible_body_dynamics), /content/body/div/div[8]/div[2]/p[2]/span[4]/span[15], line 1, column 2

Callstack:
    at (Under_Construction/Aerospace_Structures_(Johnson)/18:_Introduction_to_flexible_body_dynamics), /content/body/div/div[8]/div[2]/p[2]/span[4]/span[16], line 1, column 2

Callstack:
    at (Under_Construction/Aerospace_Structures_(Johnson)/18:_Introduction_to_flexible_body_dynamics), /content/body/div/div[8]/div[2]/p[2]/span[4]/span[17], line 1, column 2

Callstack:
    at (Under_Construction/Aerospace_Structures_(Johnson)/18:_Introduction_to_flexible_body_dynamics), /content/body/div/div[8]/div[2]/p[2]/span[4]/span[18], line 1, column 2

Callstack:
    at (Under_Construction/Aerospace_Structures_(Johnson)/18:_Introduction_to_flexible_body_dynamics), /content/body/div/div[8]/div[2]/p[2]/span[4]/span[19], line 1, column 2

Callstack:
    at (Under_Construction/Aerospace_Structures_(Johnson)/18:_Introduction_to_flexible_body_dynamics), /content/body/div/div[8]/div[2]/p[2]/span[4]/span[20], line 1, column 2

Callstack:
    at (Under_Construction/Aerospace_Structures_(Johnson)/18:_Introduction_to_flexible_body_dynamics), /content/body/div/div[8]/div[2]/p[2]/span[4]/span[21], line 1, column 2

Callstack:
    at (Under_Construction/Aerospace_Structures_(Johnson)/18:_Introduction_to_flexible_body_dynamics), /content/body/div/div[8]/div[2]/p[2]/span[4]/span[22], line 1, column 2

Callstack:
    at (Under_Construction/Aerospace_Structures_(Johnson)/18:_Introduction_to_flexible_body_dynamics), /content/body/div/div[8]/div[2]/p[2]/span[4]/span[23], line 1, column 2

Callstack:
    at (Under_Construction/Aerospace_Structures_(Johnson)/18:_Introduction_to_flexible_body_dynamics), /content/body/div/div[8]/div[2]/p[2]/span[4]/span[24], line 1, column 2

Callstack:
    at (Under_Construction/Aerospace_Structures_(Johnson)/18:_Introduction_to_flexible_body_dynamics), /content/body/div/div[8]/div[2]/p[2]/span[4]/span[25], line 1, column 2

Callstack:
    at (Under_Construction/Aerospace_Structures_(Johnson)/18:_Introduction_to_flexible_body_dynamics), /content/body/div/div[8]/div[2]/p[2]/span[4]/span[26], line 1, column 2

Callstack:
    at (Under_Construction/Aerospace_Structures_(Johnson)/18:_Introduction_to_flexible_body_dynamics), /content/body/div/div[8]/div[2]/p[2]/span[4]/span[27], line 1, column 2

Callstack:
    at (Under_Construction/Aerospace_Structures_(Johnson)/18:_Introduction_to_flexible_body_dynamics), /content/body/div/div[8]/div[2]/p[2]/span[4]/span[28], line 1, column 2

Callstack:
    at (Under_Construction/Aerospace_Structures_(Johnson)/18:_Introduction_to_flexible_body_dynamics), /content/body/div/div[8]/div[2]/p[2]/span[4]/span[29], line 1, column 2

Callstack:
    at (Under_Construction/Aerospace_Structures_(Johnson)/18:_Introduction_to_flexible_body_dynamics), /content/body/div/div[8]/div[2]/p[2]/span[4]/span[30], line 1, column 2

Callstack:
    at (Under_Construction/Aerospace_Structures_(Johnson)/18:_Introduction_to_flexible_body_dynamics), /content/body/div/div[8]/div[2]/p[2]/span[4]/span[31], line 1, column 2

Callstack:
    at (Under_Construction/Aerospace_Structures_(Johnson)/18:_Introduction_to_flexible_body_dynamics), /content/body/div/div[8]/div[2]/p[2]/span[4]/span[32], line 1, column 2

Callstack:
    at (Under_Construction/Aerospace_Structures_(Johnson)/18:_Introduction_to_flexible_body_dynamics), /content/body/div/div[8]/div[2]/p[2]/span[4]/span[33], line 1, column 2

Callstack:
    at (Under_Construction/Aerospace_Structures_(Johnson)/18:_Introduction_to_flexible_body_dynamics), /content/body/div/div[8]/div[2]/p[2]/span[4]/span[34], line 1, column 2

Callstack:
    at (Under_Construction/Aerospace_Structures_(Johnson)/18:_Introduction_to_flexible_body_dynamics), /content/body/div/div[8]/div[2]/p[2]/span[4]/span[35], line 1, column 2

Callstack:
    at (Under_Construction/Aerospace_Structures_(Johnson)/18:_Introduction_to_flexible_body_dynamics), /content/body/div/div[8]/div[2]/p[2]/span[4]/span[36], line 1, column 2

Vibrations of a coplanar frame

A bar in a coplanar frame is subject to extension and bending. In the local x-y-z coordinates of a bar, we assume that it is symmetric with the respect to the longitudinal x-z plane and the y-z plane as in article 1.8. The dynamic response of the bar for extension and bending in its y-z plane is governed by eq. ([eq18.115]) on page , which reduces to \[\begin{aligned} \label{eq18.174} &\int_{0}^{L}\left(m \ddot{v} \delta v+m \ddot{w} \delta w+m r_{x}^{2} \ddot{\phi}_{x} \delta \phi_{x}+N \delta\left(\frac{\partial w}{\partial z}\right)+M_{x} \delta\left(\frac{\partial \phi_{x}}{\partial z}\right)+V_{y} \delta \psi_{y}\right) d z\nonumber\\ &\quad =\left.\left[Q_{2} \delta v+Q_{3} \delta w+Q_{4} \delta \phi_{x}\right]\right|_{z=0}+\left.\left[Q_{8} \delta v+Q_{9} \delta w+Q_{10} \delta \phi_{x}\right]\right|_{z=L}\nonumber\\ &\qquad +\int_{0}^{L} f_{y}(z, t) \delta v(z) d z+\int f_{z}(z, t) \delta w(z) dz.\end{aligned}\] We use Euler-Bernoulli theory so that transverse shear strain \(\Psi_{y}=0\), which means the cross-sectional rotation\(.\phi_{x}=-\left(\frac{\partial v}{\partial z}\right)\). Also, neglect the rotary inertia \(m r^2_{x}\ddot{\phi}_x\) in eq. ([eq18.174]). Ignoring rotary inertia and transverse shear deformation has negligible effect on the first few modes of vibration for slender bars (Langhaar, p. 288) . Let one finite element model the bar so that \(L=h_{k}\) and \(z=\eta_{1}(\zeta) L\). The frame element is shown in its reference configuration in figure [fig18.20](a) and in rotated configuration in figure [fig18.20](b). The 6X1 displacement vector of the frame element shown in figure [fig18.20](a) is \[\begin{aligned} \label{eq18.175} \{\bar{u}(t)\}=\left[\begin{array}{@{}llllll@{}}w_{3} & v_{2} & \phi_{4} & w_{9} & v_{8} & \phi_{10}\end{array}\right]^{T}.\end{aligned}\]

The displacements are expressed as \[\begin{aligned} \label{eq18.176} \left[\begin{array}{@{}l@{}}v(\zeta, t) \\w(\zeta, t)\end{array}\right]=[N(\zeta)]\{\bar{u}(t)\},\end{aligned}\] where the interpolation matrix is \[\begin{aligned} \label{eq18.177} [N(\zeta)]=\left[\begin{array}{@{}cccccc@{}}0 & \phi_{1}(\zeta) & \phi_{2}(\zeta) & 0 & \phi_{3}(\zeta) & \phi_{4}(\zeta) \\\eta_{1}(\zeta) & 0 & 0 & \eta_{2}(\zeta) & 0 & 0\end{array}\right].\end{aligned}\] The interpolation functions for the lateral displacement \(\phi_{i}(\zeta)\), \(i=1,2,3,4\), are the Hermite cubic interpolation functions given by eq. ([eq17.140]) on page . The virtual displacements have the same interpolation as the displacements, \[\begin{aligned} \label{eq18.178} \left[\begin{array}{@{}l@{}}\delta v(\zeta, t) \\\delta w(\zeta, t)\end{array}\right]=[N(\zeta)]\{b\},\end{aligned}\] where the 6X1vector \(\{b\}\) given in the following equation is independent of the nodal displacements \[\begin{aligned} \label{eq18.179} \{b\}=\left[\begin{array}{@{}llllll@{}}b_{3} & b_{2} & b_{4} & b_{9} & b_{8} & b_{10}\end{array}\right]^{T}.\end{aligned}\] The weak form for the element is \[\begin{aligned} \label{eq18.180} &\int_{-1}^{1}\left([\delta v b w] \left[\begin{array}{@{}cc@{}}m & 0 \\0 & m \end{array}\right] \left[\begin{array}{@{}c@{}}\ddot{v} \\\ddot{v}\end{array}\right] +\left[\left(\frac{2}{L}\right)^{2} \delta v^{\prime \prime}\left(\frac{2}{L}\right) \delta w^{\prime}\right] \left[\begin{array}{@{}cc@{}}E I_{x x} & 0 \\0 & E A\end{array}\right] \left[\begin{array}{@{}l@{}}\left(\frac{2}{L}\right)^{2} v^{\prime \prime} \\[4pt] \left(\frac{2}{L}\right) w^{\prime}\end{array}\right]\right)\left[\frac{L}{2} d \zeta\right.\nonumber\\ &\quad=\{b\}^{T}\{Q\}+\int_{-1}^{1}[\delta v \delta w]\left[\begin{array}{@{}c@{}}f_{y}\left(\eta_{1}(L, t)\right) \\[4pt] f_{x}\left(\eta_{1}(L, t)\right)\end{array}\right] \frac{L}{2} d \zeta,\end{aligned}\] where the prime superscript means the derivative with respect to \(\zeta\), and the generalized force vector is \[\begin{aligned} \label{eq18.181} \{Q\}=\left[\begin{array}{@{}llllll@{}}Q_{3} & Q_{2} & Q_{4} & Q_{9} & Q_{8} & Q_{10}\end{array}\right]^{T}.\end{aligned}\] The displacement derivatives with respect to \(\zeta\) appearing in the stiffness terms of eq. ([eq18.180]) are written in the matrix form \[\begin{aligned} \label{eq18.182} \left[\begin{array}{@{}c@{}}\left(\frac{2}{L}\right)^{2} v^{\prime \prime} \\\left(\frac{2}{L}\right) w^{\prime}\end{array}\right]=\left[\begin{array}{@{}cc@{}}\left(\frac{2}{L}\right)^{2} \frac{d^{2}}{d \zeta^{2}} & 0 \\0 & \left(\frac{2}{L}\right) \frac{d}{d \zeta}\end{array}\right]\left[\begin{array}{@{}c@{}}v \\w\end{array}\right]=[D N(\zeta)]\{\bar{u}\},\end{aligned}\] where the 2X6 matrix relating the derivatives of the displacements to the displacement vector is given by \[\begin{aligned} \label{eq18.183} [D N(\zeta)]= \left[\begin{array}{@{}cc@{}}\left(\frac{2}{L}\right)^{2} \frac{d^{2}}{d \zeta^{2}} & 0 \\0 & \left(\frac{2}{L}\right) \frac{d}{d \zeta}\end{array}\right][N(\zeta)]= \left[\begin{array}{@{}cccccc@{}} 0 & \left(\frac{2}{L}\right)^2\phi_{1}^{\prime\prime} & \left(\frac{2}{L}\right)^2\phi_{2}^{\prime\prime} & 0 & \left(\frac{2}{L}\right)^2\phi_{3}^{\prime\prime} & \left(\frac{2}{L}\right)^2\phi_{4}^{\prime\prime}\\[4pt] \frac{2}{L} \eta_{1}^{\prime} & 0 & 0 & \frac{2}{L} \eta_{2}^{\prime} & 0 & 0 \end{array}\right]\!.\end{aligned}\] Substitute eqs. ([eq18.176]), ([eq18.178]), and ([eq18.182]) into ([eq18.180]) to get \[\begin{aligned} \label{eq18.184} \{b\}^{T} \int_{-1}^{1}\left([N]^{T}\left[\begin{array}{@{}cc@{}}m & 0 \\0 & m\end{array}\right][N] \frac{d^{2}}{d t^{2}}\{\bar{u}\}+[D N]^{T}\left[\begin{array}{@{}cc@{}}E I_{x x} & 0 \\0 & E A\end{array}\right][D N]\{\bar{u}\}\right) \frac{L}{2} d \zeta=\{b\}^{T}(\{Q\}+\{F\}),\end{aligned}\] where the external force vector due to the distributed loading is \[\begin{aligned} \label{eq18.185} \{F(t)\}=\int_{-1}^{1}[N]^{T}\left[\begin{array}{@{}ll@{}}f_{y}\left(\eta_{1}(L, t)\right) \\[4pt] f_{x}\left(\eta_{1}(L, t)\right)\end{array}\right] \frac{L}{2} d \zeta=\left[\begin{array}{@{}llllll@{}}F_{3} & F_{2} & F_{4} & F_{9} & F_{8} & F_{10}\end{array}\right]^{T}.\end{aligned}\] We identify the mass matrix for the reference element as \[\begin{aligned} \label{eq18.186} [\overline{M}]=\int_{-1}^{1}[N]^{T}\left[\begin{array}{@{}ll@{}}m & 0 \\0 & m\end{array}\right][N] \frac{L}{2} d \zeta=\frac{m L}{420}\left[\begin{array}{@{}cccccc@{}}140 & 0 & 0 & 70 & 0 & 0 \\0 & 156 & -44 & 0 & 54 & 26 \\0 & -44 & 16 & 0 & -26 & -12 \\70 & 0 & 0 & 140 & 0 & 0 \\0 & 54 & -26 & 0 & 156 & 44 \\0 & 26 & -12 & 0 & 44 & 16\end{array}\right].\end{aligned}\] and the stiffness matrix for the element as \[\begin{aligned} \label{eq18.187} [\overline{K}]=\int_{-1}^{1}[D N]^{T}\left[\begin{array}{@{}cc@{}}E I_{x x} & 0 \\0 & E A\end{array}\right][D N] \frac{L}{2} d \zeta=\frac{E I_{x x}}{L^{3}}\left[\begin{array}{@{}cccccc@{}}\left(L/r_{x}\right)^{2} & 0 & 0 & -\left(L/ r_{x}\right)^{2} & 0 & 0 \\0 & 12 & -6 L & 0 & -12 & -6 L \\0 & -6 L & 4 L^{2} & 0 & 6 L & 2 L^{2} \\-\left(L/r_{x}\right)^{2} & 0 & 0 & \left(L/r_{x}\right)^{2} & 0 & 0 \\0 & -12 & 6 L & 0 & 12 & 6 L \\0 & -6 L & 2 L^{2} & 0 & 6 L & 4 L^{2}\end{array}\right].\end{aligned}\] In the stiffness matrix ([eq18.187]) \(r_{x}^{2}=I_{x x}/A\). Equation ([eq18.184]) is valid for every choice of \(\{b\} \neq 0_{6 X 1}\), which leads to the equation of motion ([eq18.188]) for the reference element below. \[\begin{aligned} \label{eq18.188} [\overline{M}] \frac{d^{2}}{d t^{2}}\{\bar{u}(t)\}+[\overline{K}]\{\bar{u}(t)\}=(\{Q\}+\{F\}).\end{aligned}\] The matrix equation of motion for the rotated frame element in figure [fig18.20](b) is obtained by using the transformation matrix in eq. ([eq16.85]) on page . This transformation equation is repeated below as eq. ([eq18.189]): \[\begin{aligned} \label{eq18.189} [T]=\left[\begin{array}{@{}cccccc@{}}c & s & 0 & 0 & 0 & 0 \\-s & c & 0 & 0 & 0 & 0 \\0 & 0 & 1 & 0 & 0 & 0 \\0 & 0 & 0 & c & s & 0 \\0 & 0 & 0 & -s & c & 0 \\0 & 0 & 0 & 0 & 0 & 1\end{array}\right],\end{aligned}\] where \(c=\cos \theta\) and \(s=\sin \theta\). The generalized displacement vector for the rotated frame element is given by \(\left\{u^{(k)}\right\}=[T]\left\{\bar{u}^{(k)}\right\}\). The equation of motion for the rotated element is determined from eq. ([eq18.188]) by pre-multiplying it by \([T]^{T}\): \[\begin{aligned} \label{eq18.190} [T]^{T}[\overline{M}][T] \frac{d^{2}}{d t^{2}}\{u\}+[T]^{T}[\overline{K}][T]\{u\}=[T]^{T}(\{Q\}+\{F\}).\end{aligned}\] Equation Eq. ([eq18.190]) is written as \[\begin{aligned} \label{eq18.191} [M] \frac{d^{2}}{d t^{2}}\{u\}+[K]\{u\}=\left\{Q_{i-j}\right\}+\left\{F_{i-j}\right\}.\end{aligned}\] The 6X1 transformed generalized force vectors are \[\begin{aligned} \label{eq18.192} \begin{aligned}\left\{Q_{i-j}\right\} &=\left[\left(c Q_{3}+s Q_{2}\right)\left(c Q_{2}-s Q_{3}\right) Q_{9}\left(c Q_{9}+s Q_{8}\right)\left(c Q_{8}-s Q_{9}\right) Q_{10}\right]^{T} \\&=\left[Q_{3 i-2} Q_{3 i-1} Q_{3 i} Q_{3 j-2} Q_{3 j-1} Q_{3 j}\right]^{T}\end{aligned},\ \text{and}\\ \label{eq18.193} \left\{F_{i-j}\right\}=\left[\left(c F_{3}+s F_{2}\right)\left(c F_{2}-s F_{3}\right) F_{9}\left(c F_{9}+s F_{8}\right)\left(c F_{8}-s F_{9}\right) F_{10}\right]^{T}.\end{aligned}\] The mass and stiffness matrices for the rotated frame element are \[\begin{aligned} \label{eq18.194} \hspace*{-4pt}[M]&=\frac{m L}{420}\left[\begin{array}{@{}cccccc@{}}140 c^{2}+156 s^{2} & -16 c s & 44 s & 70 c^{2}+54 s^{2} & 16 c s & -26 c s \\[2pt] -16 c s & 156 c^{2}+140 s^{2} & -44 c & 16 c s & 54 c^{2}+70 s^{2} & 26 c \\[2pt] 44 s & -44 c & 16 & 26 s & -26 c & -12 \\[2pt] 70 c^{2}+54 s^{2} & 16 c s & 26 s & 140 c^{2}+156 s^{2} & -16 c s & -44 s \\[2pt] 16 c s & 54 c^{2}+70 s^{2} & -26 c & -16 c s & 156 c^{2}+140 s^{2} & 44 c \\[2pt] -26 c s & 26 c & -12 & -44 s & 44 c & 16\end{array}\right],\ \text{and}\\ \hspace*{-5pt}{[}K]&=\frac{EI}{L^3}\left[\begin{array}{@{}c@{\hspace*{5pt}}c@{\hspace*{5pt}}c@{\hspace*{5pt}}c@{\hspace*{5pt}}c@{\hspace*{5pt}}c@{}} \displaystyle\left(\frac{L}{r_x}\right)^2c^2+12s^2 & \displaystyle \left[\left(\frac{L}{r_x}\right)^2-12\right]cs & \displaystyle -6Ls & \displaystyle -\left[\left(\frac{L}{r_x}\right)^2c^2+12s^2\right] & \displaystyle \left[12-\left(\frac{L}{r_x}\right)^2\right]cs & \displaystyle 6Ls\\[14pt] \displaystyle\left[\left(\frac{L}{r_x}\right)^2-12\right]cs & \displaystyle 12c^2+\left(\frac{L}{r_x}\right)^2s^2 & \displaystyle -6Lc & \displaystyle \left[12-\left(\frac{L}{r_x}\right)^2\right]cs & \displaystyle -\left[12c^2+\left(\frac{L}{r_x}\right)^2s^2\right] & \displaystyle -6Lc\\[6pt] -6Ls & -6Lc & 4L^2 & -6Ls & 6Lc & 2L^2\\[6pt] \displaystyle-\left[\left(\frac{L}{r_x}\right)^2c^2+12s^2\right] & \displaystyle \left[12-\left(\frac{L}{r_x}\right)^2\right]cs & \displaystyle -6Ls & \displaystyle \left(\frac{L}{r_x}\right)^2c^2+12s^2 & \displaystyle \left[\left(\frac{L}{r_x}\right)^2-12\right]cs & \displaystyle -6Ls\\[14pt] \displaystyle\left[12-\left(\frac{L}{r_x}\right)^2\right]cs & \displaystyle -\left[12c^2+\left(\frac{L}{r_x}\right)^2s^2\right] & \displaystyle 6Lc & \displaystyle \left[\left(\frac{L}{r_x}\right)^2-12\right]cs & \displaystyle 12c^2+\left(\frac{L}{r_x}\right)^2s^2 & \displaystyle 6Lc\\[6pt] 6Ls & -6Lc & 2L^2 & -6Ls & 6Lc & 4L^2 \end{array}\right]. \label{eq18.195}\end{aligned}\]

[ex18.5]The frame shown in figure [fig18.21](a) consists of two identical bars of length L. It is modeled with two finite elements. Element \(\Omega_1\) is vertical with \({\rm c} = 0\) and \({\rm s} =1\), and element \(\Omega_2\) is horizontal with \({\rm c} = 1\) and \({\rm s} = 0\). The base of the vertical bar is subject to a horizontal harmonic displacement \(a \sin (\Omega t)\), where the amplitude is denoted by \(a\) and the driving frequently by \(\Omega\). There are nine degrees of freedom as shown in figure [fig18.21](b) for the unrestrained frame. model. The active, or unknown, generalized displacement vector is

\[\begin{aligned} \label{ex18.5a} \left\{u_{\alpha}(t)\right\}=\left[u_{3}(t)\enskip u_{4}(t)\enskip u_{5}(t)\enskip u_{6}(t)\enskip u_{9}(t)\right]^{T}.\end{aligned}\] The prescribed displacement vector is \[\begin{aligned} \label{ex18.5b} \left\{u_{\beta}(t)\right\}=\left[u_{1}(t)\enskip u_{2}(t)\enskip u_{7}(t)\enskip u_{8}(t)\right]^{T}=[a \sin (\Omega t)\enskip 0\enskip 0\enskip 0]^{T}.\end{aligned}\] The mass and stiffness matrices for element \(\Omega_1\) are \[\begin{aligned} \label{ex18.5c} &\left[M^{(1)}\right]=\frac{m L}{420} \begin{array}{@{}l@{}} {\begin{array}{@{}l@{}} \hspace*{10pt}\ddot{u}_1 \hspace*{15pt}\ddot{u}_2 \hspace*{15pt}\ddot{u}_3 \hspace*{15pt}\ddot{u}_4 \hspace*{15pt}\ddot{u}_5 \hspace*{15pt}\ddot{u}_6 \end{array}}\\ \left[\begin{array}{@{}cccccc@{}}156 & 0 & 44 & 54 & 0 & -26 \\0 & 140 & 0 & 0 & 70 & 0 \\44 & 0 & 16 & 26 & 0 & -12 \\54 & 0 & 26 & 156 & 0 & -44 \\0 & 70 & 0 & 0 & 140 & 0 \\-26 & 0 & -12 & -44 & 0 & 16\end{array}\right]\end{array}\nonumber\\ &\left[K^{(1)}\right]=\frac{E I_{x x}}{L^{3}} \begin{array}{@{}l@{}} {\begin{array}{@{}l@{}} \hspace*{8pt}u_1 \hspace*{30pt}u_2 \hspace*{28pt}u_3 \hspace*{20pt}u_4 \hspace*{30pt}u_5 \hspace*{30pt}u_6 \end{array}}\\ \left[\begin{array}{@{}cccccc@{}}12 & 0 & 6 L & -12 & 0 & 6 L \\0 & \left(L/r_{x}\right)^{2} & 0 & 0 & -\left(L/ r_{x}\right)^{2} & 0 \\6 L & 0 & 4 L^{2} & -6 L & 0 & 2 L^{2} \\-12 & 0 & -6 L & 12 & 0 & -6 L \\0 & -\left(L/ r_{x}\right)^{2} & 0 & 0 & \left(L/r_{x}\right)^{2} & 0 \\6 L & 0 & 2 L^{2} & -6 L & 0 & 4 L^{2}\end{array}\right]\end{array}.\end{aligned}\] For element \(\Omega_2\) the mass matrix is given by eq. ([eq18.186]), and the stiffness matrix is given by eq. ([eq18.187]), where the order of the columns is associated with generalized displacements \(\left[\begin{array}{@{}llllll@{}}u_{4} & u_{5} & u_{6} & u_{7} & u_{8} & u_{9}\end{array}\right]\). Assemble the mass matrices to get the 9X9 unrestrained mass matrix \(\left[M_{u}\right]=\left[M^{1)}\right]+\left[M^{(2)}\right]\) with due regard in the summation of the matrix elements from the individual mass matrices to their location in \(\left[M_{u}\right]\). Similarly, the 9X9 unrestrained stiffness matrix is assembled as \(\left[K_{u}\right]=\left[K^{(1)}\right]+\left[K^{(2)}\right]\). After assembly, the unrestrained matrix equation of motion is reordered by simultaneously interchanging the rows and columns to \(\left[\begin{array}{@{}lllllllll@{}}u_{3} & u_{4} & u_{5} & u_{6} & u_{9} & u_{1} & u_{2} & u_{7} & u_{8}\end{array}\right]\). The resulting matrix equation of motion in partitioned form is \[\begin{aligned} \label{ex18.5d} \left[\begin{array}{@{}l@{}}{\left[M_{\alpha \alpha}\right]\left[M_{\alpha \beta}\right]} \\{\left[M_{\beta \alpha}\right]\left[M_{\beta \beta}\right]}\end{array}\right]\left[\begin{array}{@{}l@{}}\left\{\ddot{u}_{\alpha}\right\} \\\left.\ddot{u}_{\beta}\right\}\end{array}\right]+\left[\begin{array}{@{}l@{}}{\left[K_{\alpha \alpha}\right]\left[\begin{array}{@{}l@{}}K_{\alpha \beta}\end{array}\right]} \\{\left[K_{\beta \alpha}\right]\left[\begin{array}{@{}l@{}}K_{\beta \beta}\end{array}\right]}\end{array}\right]\left[\begin{array}{@{}l@{}}\left\{u_{\alpha}\right\} \\\left\{u_{\beta}\right\}\end{array}\right]=\left[\begin{array}{@{}l@{}}\left\{Q_{\alpha}\right\} \\\left\{Q_{\beta}\right\}\end{array}\right].\end{aligned}\] The submatrices are \[\begin{gathered} \left[M_{\alpha \alpha}\right]=\frac{m L}{420}\left[\begin{array}{@{}ccccc@{}}16 & 26 & 0 & -12 & 0 \\26 & 296 & 0 & -44 & 0 \\-12 & -44 & -44 & 32 & -12 \\0 & 0 & 26 & -12 & 16 \\0 & 0 & 26 & -12 & 16\end{array}\right] \quad\left[M_{\alpha \beta}\right]=\frac{m L}{420}\left[\begin{array}{@{}cccc@{}}44 & 0 & 0 & 0 \\54 & 0 & 70 & 0 \\0 & 70 & 0 & 54 \\-26 & 0 & 0 & -26 \\0 & 0 & 0 & 44\end{array}\right]=\left[M_{\beta \alpha}\right]^{T}\label{ex18.5e}\\ \left[M_{\beta \beta}\right]=\frac{m L}{420}\left[\begin{array}{@{}cccc@{}} 156 & 0 & 0 & 0 \\ 0 & 140 & 0 & 0 \\ 0 & 0 & 140 & 0 \\ 0 & 0 & 0 & 156 \\ \end{array}\right] \quad\left[K_{\alpha \alpha}\right] =\frac{EI_{xx}}{L^3}\left[\begin{array}{@{}ccccc@{}} 4L^2 & -6L & 0 & 2L^2 & 0\\ -6L & 12+(L/r_x)^2 & 0 & -6L & 0\\ 0 & 0 & 12+(L/r_x)^2 & -6L & -6L \\ 2L^2 & -6L & -6L & 8L^2 & 2L^2 \\ 0 & 0 & -6L & 2L^2 & 4L^2\end{array}\right]\label{ex18.5f}\\ \left[K_{\alpha \beta}\right]=\frac{E I_{x x}}{L^{3}}\left[\begin{array}{@{}ccccc@{}}6 L & 0 & 0 & 0 \\-12 & 0 & -\left(L/ r_{x}\right)^{2} & 0 \\0 & -\left(L/r_{x}\right)^{2} & 0 & 12 \\6 L & 0 & 0 & 6 L \\0 & 0 & 0 & 6 L\end{array}\right]=\left[K_{\beta \alpha}\right]^{T} \quad\left[K_{\beta \beta}\right]=\frac{E I_{x x}}{L^{3}}\left[\begin{array}{@{}cccc@{}}12 & 0 & 0 & 0 \\0 & \left(L/r_{x}\right)^{2} & 0 & 0 \\0 & 0 & \left(L/ r_{x}\right)^{2} & 0 \\0 & 0 & 0 & 12\end{array}\right]\label{ex18.5g}\\ \left\{Q_{\alpha}\right\}=0_{5 X 1} \quad\left\{Q_{\beta}\right\}=\left[Q_{1}\quad Q_{2}\quad Q_{7}\quad Q_{9}\right]^{T}.\label{ex18.5h}\end{gathered}\] There are no generalized forces prescribed in degrees of freedom 3, 4, 5, 6, and 9. The unknown generalized reaction force vector \(\left\{Q_{\beta}\right\}\) corresponds to the prescribed generalized displacement vector \(\left\{u_{\beta}\right\}\). The matrix equation of motion for the unknown generalized displacements is \[\begin{aligned} \left[M_{\alpha \alpha}\right]\left\{\ddot{u}_{\alpha}\right\}+\left[M_{\alpha \beta}\right]\left\{\ddot{u}_{\beta}\right\}+\left[K_{\alpha \alpha}\right]\left\{u_{\alpha}\right\}+\left[K_{\alpha \beta}\right]\left\{u_{\beta}\right\}=0_{5 X 1}.\end{aligned}\]

Numerical data are taken from example [ex6.5] on page : \(a=0.03812\,\mathrm{m}\), \(A=171.014 \times 10^{-6}\,\mathrm{m}^{2}\), \(I_{x x}=124.25 \times 10^{-9}\,\mathrm{m}^{4}\), \(L=0.8\,\mathrm{m}\), \(r_{x}=0.0269546\,\mathrm{m}\), \(E=68.3 \times 10^{9}\,\mathrm{N}/\mathrm{m}^{2}\), \(\rho=2{,}710\,\mathrm{Kg}/ \mathrm{m}^{3}\), and \(m=\rho A=0.463448\,\mathrm{Kg}/\mathrm{m}\). The eigenvalue problem for the free vibrations of the frame is \[\begin{aligned} \label{ex18.5i} \left[K_{\alpha \alpha}\right]\{\phi\}-\lambda\left[M_{\alpha \alpha}\right]\{\phi\}=0_{5 X 1}.\end{aligned}\] The numerical results for the eigenvalues and corresponding eigen vectors are \[\begin{gathered} \left(\lambda_{1},\left\{\phi_{1}\right\}\right)=\left(9.99471 \times 10^{7},\left[\begin{array}{@{}c@{}}0.324108 \\0.268035 \\0.268035 \\0.324108\end{array}\right]\right) \quad\left(\lambda_{2},\left\{\phi_{2}\right\}\right)=\left(6.783 \times 10^{7},\left[\begin{array}{@{}c@{}}-1 \\0.595335 \\-0.595335 \\-1.026 \times 10^{-15} \\1\end{array}\right]\right)\label{ex18.5j}\\ \left(\lambda_{3},\left\{\phi_{3}\right\}\right)=\left(1.78362 \times 10^{7},\left[\begin{array}{@{}c@{}}1 \\ -0.105789 \\ -0.105789 \\ 0.972216\\ 1\end{array}\right]\right) \quad\left(\lambda_{4},\left\{\phi_{4}\right\}\right)=\left(2.89664 \times 10^{6},\left[\begin{array}{@{}c@{}}1 \\0.0104001 \\-0.0104001 \\-2.72413 \times 10^{-15} \\ -1\end{array}\right]\right)\label{ex18.5k}\\ \left(\lambda_{5},\left\{\phi_{5}\right\}\right)=\left(854,421,\left[\begin{array}{@{}c@{}}1 \\0.00365009 \\0.00365009 \\-0.991238 \\1\end{array}\right]\right).\label{ex18.5l}\end{gathered}\] The natural frequencies are given by \(\omega_{i}=\sqrt{\lambda_{i}}\), \(i=1,2, \ldots, 5\). \[\begin{aligned} \label{ex18.5m} \left[\begin{array}{@{}lllll@{}}\omega_{1} & \omega_{2} & \omega_{3} & \omega_{4} & \omega_{5}\end{array}\right]=\left[\begin{array}{@{}lllll@{}}9,997.36 & 8,235.9 & 4,223.3 & 1,701.95 & 924.346\end{array}\right] \mathrm{rad}/\mathrm{s}.\end{aligned}\] The prescribed displacement terms in the matrix equation of motion lead to a generalized force vector defined by \[\begin{aligned} \label{ex18.5n} \left\{\bar{Q}_{\alpha}\right\}=-\left[M_{\alpha \beta}\right]\left\{\ddot{u}_{\beta}\right\}-\left[K_{\alpha \beta}\right]\left\{u_{\beta}\right\}.\end{aligned}\] Numerical evaluation of this generalized force vector is \[\begin{aligned} \label{ex18.5o} \left\{\bar{Q}_{\alpha}\right\}=\left[\begin{array}{@{}c@{}}-3,032.78+0.00148063 \Omega^{2} \\7,581.96+0.00181714 \Omega^{2} \\0 \\-3,032.78-0.000874919 \Omega^{2} \\0\end{array}\right] \sin \Omega t.\end{aligned}\] The matrix equation of motion governing the unknown displacements is rewritten as \[\begin{aligned} \label{ex18.5p} \left[M_{\alpha \alpha}\right]\left\{\ddot{u}_{\alpha}\right\}+\left[K_{\alpha \alpha}\right]\left\{u_{\alpha}\right\}=\left\{\bar{Q}_{\alpha}\right\}.\end{aligned}\] The initial conditions are \(\left\{u_{\alpha}(0)\right\}=0_{5 X 1}\) and \(\left\{\dot{u}_{\alpha}(0)\right\}=0_{5 X 1}\). Transform from the physical displacements to the modal displacements by \(\left\{u_{\alpha}(t)\right\}=[\Phi]\{q(t)\}\), where the columns of the modal matrix \([\Phi]\) are the eigenvectors. Equations of motion in modal coordinates are \[\begin{aligned} \label{ex18.5q} \left[M_{g}\right] \ddot{q}+\left[K_{g}\right]\{q\}=\left\{F_{g}(t)\right\},\end{aligned}\] where the generalized mass, stiffness, and force vector are \[\begin{gathered} \left[M_{g}\right] =\operatorname{diag}\left[\begin{array}{@{}lllll@{}}0.02135 & 0.140492 & 0.0144387 & 0.0292596 & 0.098909\end{array}\right]\label{ex18.5r}\\ \left[K_{g}\right] =\operatorname{diag}\left[\begin{array}{@{}lllll@{}}2134.77 & 9529.59 & 257.532 & 84.7545 & 84.5099\end{array}\right] 10^{3},\nonumber\\ \left\{F_{g}(t)\right\}=\{A(\Omega)\} \sin \Omega t,\ \text{and where}\label{ex18.5s}\\ \{A(\Omega)\}=\left[\begin{array}{@{}c@{}}-1,983.5+92.0225 \times 10^{-6} \Omega^{2} \\7,275.15-463.879 \times 10^{-6} \Omega^{2} \\ -6{,}061.51+610.799 \times 10^{-6} \Omega^{2} \\ -2{,}953.93+1,499.53 \times 10^{-6} \Omega^{2} \\ 1.10169+2{,}354.52 \times 10^{-6} \Omega^{2} \end{array}\right].\label{ex18.5t}\end{gathered}\] The equation of motion for each modal coordinate is in the form \[\begin{aligned} \label{ex18.5u} \ddot{q}_{i}+\omega_{i}^{2} q_{i}=\frac{1}{M_{g}(i, i)} A_{i}(\Omega) \sin (\Omega t) \quad i=1,~2, \ldots, 5,\end{aligned}\] subject to initial conditions \(q_{i}(0)=0,\text{ and }\dot{q}_{i}(0)=0\). The solution for \(\Omega \neq \omega_{i}\) is \[\begin{aligned} \label{ex18.5v} q_{i}(t)=\frac{A_{i}(\Omega)}{M_{g}(i, i)\left(\omega_{i}^{2}-\Omega^{2}\right)}\left[\frac{-\Omega}{\omega_{i}} \sin \left(\omega_{i} t\right)+\sin (\Omega t)\right].\end{aligned}\] The solution for \(\Omega=\omega_{i}\) is \[\begin{aligned} \label{ex18.5w} q_{i}(t)=\frac{-A_{i}\left(\omega_{i}\right)}{2 M_{g}(i, i) \omega_{i}^{2}}\left[\left(\omega_{i} t\right) \cos \left(\omega_{i} t\right)-\sin \left(\omega_{i} t\right)\right].\end{aligned}\] Transform from modal coordinates to physical coordinates by \(\left\{u_{\alpha}(t)\right\}=[\Phi]\{q(t)\}\). The vertical displacement at the joint connecting the frame members is \(u_{5}(t)\), and for \(\Omega=0.9 \omega_{1}\) it is \[\begin{aligned} \label{ex18.5x} u_{5}(t)=\ &256.346 \times 10^{-6} \sin (\Omega t)-333.609 \times 10^{-6} \sin \left(\omega_{5} t\right)-151.009 \times 10^{-6} \sin \left(\omega_{4} t\right)\nonumber\\ &-47.4688 \times 10^{-6} \sin \left(\omega_{3} t\right)+ 41.6693 \times 10^{-6} \sin \left(\omega_{2} t\right)+20.1981 \times 10^{-6} \sin \left(\omega_{1} t\right).\end{aligned}\] The first term on the right-hand side of \(u_{5}(t)\) is the forced motion and has the same frequency as the excitation. The remaining five terms on the right-hand side constitute the natural motion with contributions from the five natural frequencies of the frame. A graph of the displacement response is shown in figure [fig18.9] for \(0 \leq t \leq 20 \tau\), where \(\tau\) is the period of the lowest natural frequency (i.e., \(\tau=(2 \pi)/\omega_{5}=0.006797\,\mathrm{s}\)). The displacement response exhibits the phenomena of beating.

Pure undamped beating, in general, is the combination of two sinusoids that have different but closely spaced frequencies, as the two sinusoids pass into and out of phase with each other. (Hallauer, 2016).

In the case of the frame response shown in figure [fig18.22], the driving frequency \(\Omega\) is close to the lowest natural frequency \(\omega_{5}\), and it is the combination of the two sinusoids \(\sin (\Omega t)\) and \(\sin \left(\omega_{5} t\right)\) that are the source of beating.

For \(\Omega=\omega_{5}\), the displacement \(u_{5}(t)\) is \[\begin{aligned} \label{ex18.5y} u_{5}(t)=&-0.406236 t \cos (\Omega t)+319.34 \times 10^{-6} \sin \left(\omega_{5} t\right)-158.115 \times 10^{-6} \sin \left(\omega_{4} t\right)\nonumber\\ &-52.3109 \times 10^{-6} \sin \left(\omega_{3} t\right)+ 45.9088 \times 10^{-6} \sin \left(\omega_{2} t\right)+22.3042 \times 10^{-6} \sin \left(\omega_{1} t\right).\end{aligned}\] A graph of the displacement response is shown in figure [fig18.23] for \(0 \leq t \leq 10 \tau\), where \(\tau\) is the period of the lowest natural frequency. As shown in the graph, the amplitude of the vibrations for \(\Omega=\omega_{5}\) increase linearly in time, which is the condition of resonance.

[sec18.10] Bathe, K-L. Finite Element Procedures in Engineering Analysis. Englewood Cliffs, NJ: Prentice-Hall Inc., 1982.[Bathe]

Boresi, A. P. Elasticity In Engineering Mechanics. Englewood Cliffs, N.J: Prentice-Hall, Inc. 1965, p. 34.

Craig, Roy R., Jr. Structural Dynamics: An Introduction to Computer Methods. New York: John Wiley & Sons, Inc., 1981, pp. 212–217, 321–340, 455–464.[Craig]

Hallauer, W. L., Jr., Introduction to Linear, Time-Invariant, Dynamic Systems for Students of Engineering. Blacksburg, VA: Self-published, 2016. http://hdl.handle.net/10919/78864.

Isakowitz, Steven J. Space Launch Systems, 2d ed. Updated by Jeff Samella.Washington DC: American Institute of Aeronautics and Astronautics, 1995, pp. 201–218.

Langhaar, H. L. Energy Methods in Applied Mechanics. New York: John Wiley, and Sons, Inc., 1962.[Langhaar]

Qu, Zu-Qing. Model Order Reduction Techniques, With Application to Finite Element Analysis. London: Springer-Verlag, 2004.[Qu2004]

Sarafin, Thomas, P., ed. Spacecraft Structures and Mechanisms – From Concept to Launch. Torrance, CA, and Dordrecht, The Netherlands: Microcosom Press, and Kluwer, 1995, pp. 49 & 50.[Sarafin1995]

Schiesser, W. E. Computational Mathematics in Engineering and Applied Science. Boca Raton, FL: CRC Press, 1994, Chapter 2.

Szabó, B., and I. Babuska. Finite Element Analysis. New York: John Wiley & Sons, 1991, pp. 163–166.

Practice exercises

1. The Lagrangian for a three-degree-of-freedom model of Atlas I is \[\begin{aligned} \label{eq18.196} L=\frac{m_{1}}{6}\left[\left(\dot{w}_{1}\right)^{2}+\dot{w}_{1} \dot{w}_{2}+\left(\dot{w}_{2}\right)^{2}\right]+\frac{1}{2} m_{2}\left(\dot{w}_{2}\right)^{2}+\frac{1}{2} m_{3}\left(\dot{w}_{3}\right)^{2}-\frac{1}{2} k_{12}\left(w_{1}-w_{2}\right)^{2}-\frac{1}{2} k_{23}\left(w_{3}-w_{2}\right)^{2},\end{aligned}\] where \(w_{1}(t)\) is the displacement at the bottom of the booster, \(w_{2}(t)\) is the displacement at the top of Centaur, and \(w_{3}(t)\) is the displacement of the payload. These displacements are defined with respect to the equilibrium state. The combined mass of the booster and Centaur is denoted by \(m_1\), the mass of the fairing by \(m_2\), and the mass of the payload by \(m_3\). Masses are determined from the weight data given in “Description of Atlas I” on page . The spring stiffness \(k_{12}\) and \(k_{23}\) are listed in “Step 1: Equations of motion about equilibrium.” on page . Lagrange’s equations of motion are \[\begin{aligned} \label{eq18.197} \frac{d}{d t}\left(\frac{\partial L}{\partial \dot{w}_{1}}\right)-\frac{\partial L}{\partial w_{1}}=R_{1} \quad \frac{d}{d t}\left(\frac{\partial L}{\partial \dot{w}_{2}}\right)-\frac{\partial L}{\partial w_{2}}=0 \quad \frac{d}{d t}\left(\frac{\partial L}{\partial \dot{w}_{3}}\right)-\frac{\partial L}{\partial w_{3}}=0,\end{aligned}\] where \({R}_1 = 77,100\). lb is the net thrust. Determine the maximum payload load factor during the initial instants of lift off. Partial answer: the value and its associated eigenvector for the smallest elastic mode is \[\begin{aligned} \label{eq18.198} \left(\lambda_{2},\left\{\phi_{2}\right\}\right)=\left(15,710.5,\left[\begin{array}{@{}c@{}}-0.795164 \\0.760969 \\1\end{array}\right]\right).\end{aligned}\] The eigenvector is normalized such that the magnitude of its largest component is a positive one.

2. Determine the natural frequencies in Hz and the corresponding modal vectors for the five-bar, pin-jointed truss shown in figure [fig18.24]. Normalize the modal vectors such that the largest component in the vector is a positive 1. Sketch the mode shapes.

3. The three-bar truss in example [ex18.3] on page is subjected to the following initial conditions \[\begin{aligned} \label{eq18.199} \left[\begin{matrix}u_{1}(0) \\u_{2}(0)\end{matrix}\right]= \left.\left[\begin{matrix}0 \\0\end{matrix}\right] \quad\left[\begin{matrix}\displaystyle\frac{d u_{1}}{d t} \\[8pt] \displaystyle\frac{d u_{2}}{d t}\end{matrix}\right]\right|_{t=0}= \left[\begin{matrix}0 \\1\end{matrix}\right] \mathrm{m}/\mathrm{s}.\end{aligned}\]

   1. Determine the generalized mass matrix \(\left[M_{g}\right]\), the generalized stiffness matrix \(\left[K_{g}\right]\), and the initial conditions in modal coordinates (begin with eqs. (d) and (e)).

   2. Determine the solution in modal coordinates and in physical coordinates.

   3. Determine the transient bar forces \(N_{1-2}(t), N_{1-3}(t), \text{and } N_{1-4}(t)\).

   4. Plot the bar forces found in part (c) for \(0 \leq t \leq 0.015\,\mathrm{s}\).

4. Model the cantilever beam in example [ex18.4] on page with two equal length elements as shown in figure [fig18.25].

   1. Determine the natural frequencies in Hz and the corresponding modal vectors. Normalize the modal vectors such that the tip displacement \(u_5\) is equal to one in each mode. Refer to table [tab18.4].

   2. Determine the percentage error of each frequency with respect to the exact frequency from the continuous beam vibration analysis.

   3. Plot the lateral displacement of the beam, \(0 \leq z \leq 0.8\,\mathrm{m}\), for each mode using eq. ([eq17.69]) on page and matrix \(\left[G_{a q}\right]\) from the last three rows of eq. ([eq18.162]).

5. The uniform beam shown in figure [fig18.26] is simply supported at each end. The material and geometric properties are the same as those given in example [ex18.4],

   1. Determine the first two natural frequencies in Hz for the beam modeled with one element. Use the condensed mass matrix ([eq18.167]) and stiffness matrix ([eq17.106]) on page .

   2. Compute the percent discrepancy of the frequencies with respect to the continuous beam solution. The frequencies for the continuous beam vibration analysis in rad/s are listed in Graig (1981) as \(\omega_{n}=\left(\frac{n \pi}{l}\right)^{2} \sqrt{\frac{E I_{x x}}{\rho A}}\) \(n=1,2\).