## Abstract

Development and employment of dynamic elements can result in substantial gain in solution convergence for vibration analysis, when compared with the usual finite element discretization. Efficient solution of the associated quadratic matrix eigenproblem is crucial in achieving a superior and relatively economical solution. This paper first describes a novel eigensolution of the quadratic matrix equation by a progressive simultaneous iteration method. Free vibration analysis results of a rectangular prestressed membrane are next presented in detail that provides an assessment of relative solution convergence and computing expenses of the two idealization procedures.

## 1. Introduction

An accurate and economical solution of the structural vibration problem is a vital preliminary for the subsequent dynamic response analysis of aerospace vehicles such as Hyper-X, the space shuttle and others. In this process, first the continuous structural system is discretized by a suitable idealization technique such as the finite element method (FEM), yielding a set of large sparse matrix equations. These equations are then solved by a suitable eigensolver to yield the natural frequencies and associated modes, which may then be used for an effective response calculation that certifies flight worthiness. Finite dynamic elements (FDEs) were first developed for line elements by Przemieniecki (1966). Subsequently, such elements were developed for various continuums such as prestressed membrane (Gupta 1976), plane (Gupta 1979) and shell (Martin & Lung 1991) among others.

In the standard FEM, the continuous displacement field vector within an element is derived in terms of the unknown nodal values, which is valid only for static loading. For dynamic loading, the displacement field is a function of the entire time history of the nodal displacements. However, for the case of structural free vibration involving harmonic motions, such a relationship is expressed as
1.1
where the shape function matrix **a** is also a function of the natural frequencies *ω*, resulting in frequency-dependent stiffness and inertia matrices. Solution of the resulting eigenproblem is uneconomical and rather intractable, since these matrices are functions of the frequencies that are not known, a priori. To circumvent such a difficult proposition, the shape function matrix is expanded in a series form in ascending powers of *ω* as
1.2
that results in the following expressions for the element stiffness and inertia matrices
1.3
and
1.4
as well as the corresponding system matrices
1.5
and
1.6
This results in the equation of free vibration
1.7**q** being the amplitude of **U**, and **K**_{0} and **M**_{0} are static stiffness and inertia matrices, respectively, the rest being higher-order dynamic correction terms, each of order *N*. In the dynamic element method (DEM; Gupta 1976), the series in equation (1.7) is usually truncated after the third term yielding the quadratic matrix eigenvalue problem
1.8
where **A**=**K**_{0}, **B**=**M**_{0}−**K**_{2} and **C**=**M**_{2}−**K**_{4} are symmetric and positive definite matrices of order *N*, in which *λ*=*ω*^{2}, the solution of which yields the desired roots and vectors. For most practical problems, these matrices are highly sparse in nature and a desirable solution algorithm must necessarily exploit such sparsity, enabling solution of large, complex problems in an efficient and economical fashion. An eigenvalue solution algorithm for equation (1.8) and an associated code, based on a progressive simultaneous iteration (PSI; Gupta & Lawson 1999) method, was developed that testifies to its superior convergence characteristics, which renders the FDE (Gupta 1976) method, involving the **A**, **B** and **C** matrices, highly efficient and attractive. The first FDE for a continuum was developed as a prestressed membrane element (Gupta 1976). In this procedure, the progressive process involves only the unconverged roots during the various stages of its iteration, thereby effecting considerable economy in solution time. Relevant numerical examples involving prestressed membrane elements are also presented in detail that testifies to the efficacy of the techniques presented herein.

## 2. Solution of the quadratic matrix eigenproblem

A solution of the quadratic matrix equation may be achieved by first reducing equation (1.8) into a system that involves a single matrix (Przemieniecki 1966, 1968) and solving by any standard procedure.

However, the solution of such a full matrix will be prohibitive for large, practical problems. The matrices **A**, **B** and **C** are banded and sparse for most problems, and an effective eigensolution needs to exploit such sparsity to yield an economical solution. Such a solution can be effected by first rearranging equation (1.8) in the following form:
2.1
with and which is of the form
2.2
where the matrices **E** and **F** are both symmetric and sparse, **E** being also positive definite; the associated roots and vectors are real in nature. Solution of the quadratic matrix equation has been achieved earlier employing a combined Sturm sequence and inverse iteration (Gupta & Meek 2003), Lanczos (Gupta & Meek 2003) and subspace iteration (Gmuer 1990) methods. A simultaneous iteration (SI; Jennings & Orr 1971; Bathe & Wilson 1973) method with progressive features is described herein that pertains to the current quadratic matrix equations and proves to be comparatively much more efficient and economical for solution of large-order problems. In particular, this procedure fully exploits matrix sparsity even within the matrix bandwidth.

A PSI procedure was developed earlier for solution of the usual linear eigenproblem (Gupta & Lawson 1999) and also quadratic eigenproblem for gyroscopic (Gupta & Lawson 2001) systems involving a pair of Hermitian and real symmetric sparse matrices of order 2*N*×2*N*. An associated numerical algorithm for the solution of the current quadratic matrix eigenproblem, equation (2.2), is given below.

*Step 1*. Generate an *NRT* (*NRT*>*NR*) set of random vectors
each being of order 2*N*, *NR* being the number of desired roots.

*Step 2*. Perform a Cholesky factorization, initially once only
in which **L** is a unit lower triangular matrix, **D** being diagonal.

*Step 3*. Start of iteration; at a typical *i*th step, solve
by the usual back substitution method, progressively, using only the unconverged roots.

*Step 4*. Estimate the magnitude of the *j*th unconverged root *λ*^{j} as
and check for convergence if this is not the first iteration, using the test
EPSN being a suitable root convergence accuracy parameter.

Assuming that *NR*1 roots have passed the convergence test and if *NR*1=*NR*, then go to step 7.

*Step 5*. Perform progressive vector orthonormalization with respect to the **E** matrix as below:
and
employing the *NR*1–*NR* unconverged roots.

*Step 6*. Go to step 3 and continue iteration using the computed set of vectors
*Step 7*. Form the reduced eigenvalue problem
and using the progressive strategy, compute the unconverged roots only, employing the Rayleigh–Ritz method, where being an approximation of *λ* and in which
*Step 8*. Convergence test for each typical *j*th root under consideration is performed next
*Step 9*. Recompute the eigenvectors
and set
*Step 10*. In step 8, if convergence is achieved for only *NR*1 roots, *NR*1<*NR*, then go to steps 3, 5 and 7–9 only; otherwise go to step 11.

*Step 11*. End of iteration.

Progressive saving in the solution process is effected in operations in steps 3, 5, 7 and 9 resulting in a very efficient quadratic matrix eigenproblem solution procedure when compared to other similar existing techniques. For further reference, the process of checking root norm convergence involving steps 3–6 will be called phase 1 operation. Phase 2 is designated for the process that involves steps 3, 5 and 7–10. In phase 1, an initial estimate of the roots and vectors is achieved that accelerates into an accurate estimate in phase 2. A typical value around 0.01–0.05 is set for the accuracy factor EPSN for vector norm check in phase 1, whereas EPS for the Rayleigh–Ritz operation is set to a rather small value as 0.00001 or so that yields a sufficient accurate solution in phase 2. It may be noted that the usual SI approach (Gmuer 1990) does not include the progressive features. Typically, six iterations or so are needed in phase 1, that is followed by about two iterations in phase 2 to complete the solution.

## 3. Prestressed membrane dynamic element

Figure 1 shows a uniformly prestressed thin membrane element of thickness *h*. The differential equation of motion of the membrane can be written as
3.1
in which *u* is the out-of-plane deformation at any point on the membrane, *σ**h* the constant uniform tensile force per unit length and *ρ* the membrane mass per unit area. The membrane deformation *u* can be expressed in series form in terms of its nodal values as
3.2
where **a** is the shape function vector (Gupta & Meek 2003). Substituting equation (3.2) into equation (3.1) yields the equation of motion (Gupta 1976)
3.3
Appropriate solution (Gupta & Meek 2003) of equation (3.3) yields the various terms of the shape function series
in which
3.4
3.4
and **a**_{1}=0 for the present problem. Similarly the strain–displacement relationship is expressed as
3.6
which then yields the individual stiffness and inertia matrices **k**_{0}, **k**_{2}, **k**_{4}, **m**_{0} and **m**_{2} by standard procedures resulting in the element matrix equation
3.7
which are then assembled for the entire membrane yielding equation (1.7).

## 4. Numerical examples

This example pertains to a simply supported rectangular membrane problem with the following details: side length, *a*=8; side length, *b*=12; thickness, *h*=1; density, *ρ*=1; uniform tension, *σ*=1.

The membrane is discretized with varying mesh sizes to study solution convergence characteristics. Table 1 provides details of such an analysis for the first 20 roots with mesh sizes varying from 4×4 to 100×100. Figure 2 presents the pattern of such a convergence study for the computed roots. In an effort to provide the relative efficacies of the usual FEM and the current DEM procedures, their respective computing solution times are depicted in figure 3; the MAINIB refers to the matrix generation module of the STARS (Gupta 1997) code used for all solutions. Corresponding mode shapes are depicted in figure 4.

All computations were performed with a single CPU, using Dell Precision M4300, Intel Core 2 Duo (T9500), 3.5 GB RAM and MS Windows XP Professional with SP3. The related algorithms were incorporated in the software STARS (Gupta 1997), which was used to solve the example problems, presented herein.

## 5. Concluding remarks

A dynamic element discretization procedure generally results in better solution convergence when compared with the usual FEM and hence is the justification (Voss 1987) for further development in this area; this involves generation of new type of elements and very importantly development of efficient quadratic matrix eigensolution techniques and tools. This paper describes a novel quadratic matrix eigensolver in detail, brief details of a prestressed membrane dynamic element and elaborate solution results pertaining to an example problem, affording comparison of the FEM and DEM techniques.

The numerical example relates to the free vibration analysis of a rectangular prestressed membrane; dynamic element formulation for other various elements such as shell, plates and solids are given in Gupta (1976) and Gupta & Meek (2003), among others. It is apparent from solution results that the total DEM solution is considerably faster than the corresponding FEM solution and the difference grows larger with increasing degree of discretization. It is interesting to note that although the matrix generation module (MAINIB) in DEM is more time-consuming than that in FEM, their corresponding PSI solutions are opposite in nature. This is because of much faster PSI solution convergence for the more refined, higher-order DEM idealizations. A highly efficient sparse matrix scheme using a combination of the minimum fill method along with a multifrontal decomposition procedure ensures the most economical solution in terms of both storage and CPU time.

In the element formulation, it is known (Gupta 1976) that **k**_{0} is a direct function of the prestress *σ*, the higher-order matrix **k**_{4} is a function of 1/*σ* and the rest of the matrices are not influenced by *σ*. A series of test runs were made with differing prestress values *σ*, and the frequencies were found to be modified by mostly the square root of *σ*, indicating the much more dominant role of **k**_{0}.

Alternative formulation (Dai & Liu 2007; Nguyen-Xuan *et al*. 2008) of quadrilateral elements show improved convergence rates for static and dynamic problems. It will be beneficial to develop such elements for more commonly used triangles and tetrahedrons for ready applications in aerospace engineering.

## Acknowledgements

The authors thank S. Choi for his sincere effort in the preparation of this paper.

## Footnotes

- Received August 21, 2009.
- Accepted September 23, 2009.

- © 2009 The Royal Society