HQQA (DUI-h . THESlS f) X,‘ f\ é—ULL/ ”BRA? Michigan Stat; University l This is to certify that the thesis entitled FICTITIOUS DOMAIN SOLUTION OF FREQUENCY RESPONSE PROBLEMS presented by Sudarsanam C hellappa has been accepted towards fulfillment of the requirements for MS Mechanical Engineering degree in I V ' Major pilifessor Date 6‘30““ loop 0-7639 MS U is an Affirmatiw Action/Equal Opportunity Institution PLACE IN RETURN BOX to remove this checkout from your record. TO AVOID FINES return on or before date due. MAY BE RECALLED with earlier due date if requested. DATE DUE DATE DUE DATE DUE woo WM 14 FICTITIOUS DOMAIN SOLUTION OF FREQUENCY RESPONSE PROBLEMS By Sudarsanam Chellappa A THESIS Submitted to Michigan State University In partial fulfillment of the requirements For the degree of MASTER OF SCIENCE Department of Mechanical Engineering 2000 ABSTRACT FIC'I'I’I‘IOUS DOMAIN SOLUTION OF FREQUENCY RESPONSE PROBLEMS. by Sudarsanam Chellappa An alternative solution strategy for frequency response problems in elasticity is proposed. This strategy uses a mesh-less finite element method, derived in a periodic domain setting, extended to allow the analysis of arbitrary geometries with the use of a fictitious domain. This solution technique is designed to solve problems with a large number of degrees of freedom and it requires the use of iterative solvers where preconditioning is necessary for reasonable convergence rates. A preconditioned conjugate gradient-type solver is used with a preconditioner specially designed for this problem so that the convergence of the iterative algorithm is quick and insensitive to the size of the problem. Some examples illustrating this approach are provided. To my parents ACKNOWLEDGEMENTS I would like to thank my advisor Dr. Alejandro Diaz. Your guidance and support has been invaluable. I would also like to thank the members of my committee: Dr. Andre Benard and Dr. Alan Haddow. Thank you for the time spent in reviewing this document and providing valuable suggestions. iv TABLE OF CONTENTS LIST OF FIGURES".-- - - - - . _ _ _ _ - - - - vii Chapter 1 Introduction - - - - - - .............. - - . - . - - --1 Chapter 2 Fictitious Domain Approach - - - - - — --7 2.1 Introduction .............................................................................................. 7 2.2 Field Equations of Elastostatics .............................................................. 8 2.3 Problem Statement .................................................................................. 9 2.4 Summary ................................................................................................ 12 Chapter 3 Applying Fictitious Domain Methods To Solve Frequency Response Problems in Two Dimensions - - - _ _ - _ -- _ -- _ - 14 3.1 Introduction ............................................................................................ 14 3.2 Field Equations in Elastodynamics ....................................................... 15 3.3 Problem Statement ................................................................................ 16 3.4 Problem Discretization .......................................................................... 18 3.5 Solution Strategy ................................................................................... 20 3.6 Preconditioning ...................................................................................... 21 3.7 Derivation of the Preconditioner ........................................................... 22 3.8 Inversion of the Preconditioner ............................................................. 24 3.9 Examples ................................................................................................. 26 3.9.1 Example 1 ..................................................................................... 26 3.9.2 Example 2 ..................................................................................... 33 3.9.3 Example 3 ..................................................................................... 36 3.9.4 Example 4 ..................................................................................... 39 Chapter 4 Solving Frequency Response Problems in Three Dimensions ............ 42 4.1 Introduction ............................................................................................ 42 4.2 Problem Statement ................................................................................ 42 4.3 Problem Discretization in Three Dimensions ....................................... 44 4.4 Solution Strategy ................................................................................... 46 4.5 Derivation of the Preconditioner ........................................................... 47 4.6 Examples ................................................................................................ 51 4.6.1 Example 1 ..................................................................................... 51 4.6.2 Example 2 ..................................................................................... 54 Chapter 5 Conclusions ................................................................................................. 57 APPENDIX ................................................................................................. 59 BIBLIOGRAPHY ............................................................................................ 61 List of Figures 2.1 Original Problem on Q, ............................................................................... 9 2.2 Fictitious Domain Problem ........................................................................ 10 2.30 -Periodic domain ...................................................................................... 1 1 2.4 A circular domain at various resolutions .................................................. 13 3.1 Frequency response problem on an arbitrary domain .............................. 16 3.2 Frequency response problem on Q ........................................................... 17 3.3 Structure of block-circulant matrices ........................................................ 23 3.4 Cantilever beam with a tip load ................................................................ 26 3.5 Comparison of fictitious domain solution with exact result ..................... 27 3.6 Performance of the solver for various preconditioners ............................. 28 3.7 Performance of the solver for difl'erent frequencies and resolutions ........ 30 3.8 Frequency response plots ........................................................................... 32 3.9 One quarter of a plate with a center hole .................................................. 33 3.10 Performance of the solver for difl'erent frequencies and resolutions ....... 34 3.11 Frequency response plot ............................................................................ 35 3.12 A curved beam, fixed at one end and a tip load at the other .................. 36 3.13 Performance of the solver for difi'erent frequencies and resolutions ....... 37 3.14 Frequency response plot ............................................................................ 38 3.15 An arbitrary geometry fixed at two ends and a force in the middle ....... 39 3.16 Performance of the solver for different frequencies and resolutions ....... 40 3.17 Frequency response plot ............................................................................ 41 4.1 Example of a voxel discretization. ............................................................. 44 4.2 Constrained and fictitious degrees of freedom .......................................... 47 4.3 Structure of two-block circulant matrices ................................................. 48 4.4 A rectangular bar with a tip load .............................................................. 51 4.5 Performance of the solver for difl'erent frequencies and resolutions ........ 52 4.6 Frequency response plots ........................................................................... 53 4.7 One half of a crankshaft, constrained along its length ............................. 54 4.8 Performance of the solver for difl‘erent frequencies and resolutions ........ 55 4.9 Frequency response plot ............................................................................. 56 vii Chapter 1 Introduction The goal of solving frequency response problems is to determine the response of a structure subject to a periodic forcing. This is of interest in design problems in which the mechanical and structural systems are subjected to harmonically varying external loads caused by reciprocating power train or other rotating machinery such as motors, fans, compressors and forging hammers. When a machine or any structure oscillates in some form of periodic or random motion, the motion generates alternating pressure waves that propagate from the moving surface at the velocity of sound. These motions with frequencies in the audible range can cause discomfort to the people nearby. For instance, in an automobile the input forces transmitted from the road and the power train excite the vehicle compartment boundary panels. Airplane body and wing structures are also subjected to a harmonic load transmitted from the propulsion system causing noise problem, cracks, fatigue failure of the tail shaft and discomfort to the crew. This provides a motivation for us to study the frequency response of a structure. The equations associated with this problem are the partial differential equations associated with linear elasticity. Typically, these problems are solved using the finite element method. In this approach, the problem domain is discretized into small elements and the local response is approximated using polynomial shape functions associated with a particular choice of a finite element. This results in discretized set of linear equations of the form Ax =b, where A is the coeficient matrix resulting from the finite element discretization, b is the force vector and x is the response. Typically, a solution to the above set of linear equations is obtained by inverting the coemcient matrix (A) using some standard techniques like LU-decomposition. These techniques however require that the matrix be stored in care. There are difliculties, however, that arise when the size of the problem is very large as in the case when a very fine resolution of the geometry or the material distribution is necessary. As the problem size increases, the memory requirements and CPU time necessary to obtain finite element solutions dominate the overall process. For large-scale problems, high computational demands often make the problems impractical to solve using reasonable computational resources. In very large problems one must resort to iterative procedures that do not require storage of the coeflicient matrix. Iterative solvers need more CPU time than direct methods but use substantially less in-core memory. Now, the rate of convergence becomes an important issue. It is well known that when the standard finite element approach is used, the rate of convergence of iterative schemes such as preconditioned conjugate gradient methods decays as the resolution increases, eventually rendering the approach useless [I]. This motivates the search for a solution scheme that is better suited for these kinds of large-scale problems. This dissertation presents a new method to solve large-scale frequency response problems by modifying the standard finite element approach using mesh-less techniques that are specially tailored to solve these problems. Our goal is to develop a new analysis technique to determine approximate solutions to the elasticity equations associated with the frequency response problem. We seek an analysis technique that performs analysis at rates that are independent of resolution. To meet this goal we turn to the fictitious domain technique. Our approach is based on mesh-less methods for analysis that accurately predicts the elastic response of highly detailed mechanical components. The discretization is developed from an image of the component embedded into a simple fictitious domain: a square in two dimensions and a cube in three dimensions. The fictitious domain is used to simplify the analysis by embedding the originally complex geometry of the structure in a simple regular domain. In two dimensions, this process corresponds to a pixel level discretization of the component, while in three dimensions the discretization is similar to a voxel discretization (the three-dimensional equivalent to pixels). Due to potentially high resolution of the image discretization for reasons of accuracy, an iterative preconditioned conjugate gradient type solver is incorporated to solve the resulting linear equations. A special preconditioner that results in convergence rates that are insensitive to the resolution of the image discretization is used. Each pixel or voxel in the image discretization is associated with a single finite element, thus avoiding the time consuming meshing process often required for standard finite element techniques. One important property of this numerical scheme is that the coding for the solution of the boundary value problem is independent of the geometry of the boundary. Normally in finite element approaches to such a boundary value problem, one has a finite element grid, which is adapted to the boundary. This is not necessary for this scheme, which gives it great flexibility in applications. These mesh-less techniques have been used efi‘ectively by DeRose and Diaz for solving topology optimization problems in elastostatics using a wavelet- Galerkin approach [2,3,4]. In this approach, bases constructed from refinable shift-invariant wavelets were combined with traditional Galerkin techniques to perform the required analysis in the generalized topology optimization problems. The boundary conditions were implemented in the form of Lagrange multipliers. This resulted in a discretized set of linear equations whose coefficient matrix is real, symmetric and positive semi-definite. They solved this with a Preconditioned Conjugate Gradient solver by incorporating a preconditioner based on the work of Bramble and Pasciak [5]. Earlier work on the solution of partial difi'erential equations using fictitious domain techniques has been mainly due to R.O. Wells, X. Zhou, R. Glowinski, T.W. Pan, et a1, some of which are discussed briefly here. In 1993, Wells and Zhou [6] published a modified form of the classical penalty method for solving a Dirichlet boundary value problem using wavelets. This was based on the fact that one could expand the boundary measure under the chosen basis, which led to a fast, approximate calculation of boundary integrals. Glowinski et al in 1995 [7] presented a new fictitious domain formulation for the strongly elliptic boundary value problem with Neumann boundary conditions for a bounded domain in a finite-dimensional Euclidean space with a smooth boundary and extended that to a larger rectangular domain with periodic boundary conditions. They have presented the application of this method using both wavelet—Galerkin and finite element approximation schemes. They have also shown that the convergence rates of both schemes are comparable. In 1996, Glowinski et al [8] presented a wavelet multigrid preconditioner for the conjugate gradient method that acted as an efiicient solver for the linear system arising fi'om a wavelet-Galerkin discretization of a Dirichlet boundary value problem via a penalty/fictitious domain formulation. They have also described some numerical experiments to confirm the emciency of the iterative solver. The remainder of the dissertation is organized as follows. Chapter 2 discusses the fundamentals of the fictitious domain approach using a simple elastostatic case for better understanding. Chapter 3 discusses the application of the fictitious domain approach to solve frequency response problems in two dimensions. Chapter 4 talks about solving three-dimensional frequency response problems and describes a slightly modified version of the solution technique used in the case of two-dimensional problems. Conclusions about this method are provided in Chapter 5. Chapter 2 Fictitious Domain Approach 2.1 Introduction This chapter discusses the basics of the fictitious domain approach and its application to solve problems in linear elasticity. The goal of the fictitious domain approach is to solve problems with arbitrary domains by embedding them into a simple domain, a square in two dimensions or a cube in three dimensions, and solving instead an equivalent problem in the simplified geometry. Hopefully, recasting the problem in the simplified domain facilitates the use of efficient solution techniques without compromising the accuracy in the domain of interest. The following sections discuss this approach. We have shown its application to solve the simple static elasticity problems so as to facilitate better understanding of the approach. For a more complete discussion of the fictitious domain approach see the work of Wells, Glowinski et a1 [6,7,8]. 2.2 Field Equations of Elastostatics - Equilibrium Equations: UN +f, =0 (2.1) o Strain-Displacement relationship: 61 = all” + u“) (2.2) l o Stress-Strain relationship: 0'” = D be], (2-3) ., Here, a is the elastic stress tensor, 8 is the elastic strain tensor, 1’ is the body force per unit volume, u is the displacement field and D is the elastic material tensor. Combining the above three equations we get D blah/k + f. = 0 (2.4) I} Multiplying this equation with a weight function (w) and integrating over the volume we get it into the familiar weak form, which is [De(u)g(w)dV = [1”de for all w, e H'(V) (2.5) V V Here H‘ (V) denotes the Sobolev space of 1.2 functions in V whose It" derivative is also in L2(V). L2 is the set of functions f : R —+ IR such that the integral of f2 in V is finite, i.e., f e L’(V) ¢::> {/de < 00 (2.6) 2.3 Problem Statement We consider a two-dimensional elasticity problem on an arbitrary domain (1,, with displacement boundary conditions specified on F, a subset of 69,, as shown in Figure 2.1. F Figure 2.] Original problem on Q, The above problem can be expressed as Find v, 6 V0, suchthat IDm£(v)£(w)dQ = L0,”) {0, all "I 6 Val (2.7) a. where Ln. (w) = [1"de (2.8) 02 and w,eVo,E{w,eH‘(92):w,=0inF;602} (2.9) Note that v=(v,,vy), f=(f,,fy) and D02 =Dn’(x,y). Now we embed 9, into a regular (fictitious) domain 0, to get a new heterogeneous domain Q(=[0,d]x[0,d])=9, US), with (2, n0, =¢ as shown in Figure 2.2 and introduce another constraint, such that the solution is Q- periodic. Figure 2.2. Fictitious dormin problem The problem statement can be rewritten as Find u, e Kn such that [Dne(u)e(w)dfl = If’wdfl for all w, 6 KO (2'10) 0 Q where Ka£{u,eVn: u,=0inl‘;602} (2.11) and Va 5 {u, e H'(O): u, is O-periodic} (2.12) Now, this problem looks exactly like a problem produced by the homogenization procedure for computing the homogenized coefficients of composites of essentially different components with a periodic structure. The Q -periodicity is introduced so as to facilitate a particular solution technique (to be discussed at a later stage) specially tailored to such kind of problems. The meaning of Q-periodicity can be mathematically stated as shown in (2.13) and illustrated in Figure 2.3. ".25. .3 A function v :Q -) IR, where Q is [0, d]x [0,d], is Q-periodic , if v(x,y)=v(x+d,y+d) forall(x,y)eQ (2-13) III ‘ Ci (3* (j C} C (j d c. land‘— Figure 2.3. 0 -Periodic domain Convergence proofs and error estimates relevant to this method can be found in the work of Bakhvalov and Knyazev [9]. ll Our goal is to obtain solutions It to the Q-periodic problem that are sumciently close to v in 02. To accomplish this we would like the material in the fictitious domain to have little effect, if any, on the solution in Q, . In the elasticity problem this is guaranteed by requiring that the fictitious domain be composed of a very weak material. The material tensor in the heterogeneous domain 0 is given as DH x,yeQ, (2.14) Do, (x. y) x. y e 92 Dfl(xay) ={ Similarly, the mass-density distribution is given as pm xv)’ e QI (x. )= (2.15) pa y {p0, (x9 y) x9y E 02 DWIIII and pm are properties of the material in the fictitious domain. Typically, we assume O< D,“ < DO, and 0 < p“ « pa: (2.16) and that I)“.It is isotropic. 2.4 Summary We have now successfully formulated an equivalent problem in a regular domain in place of the original problem. This has been done at the price of having significantly increased the size of the problem this fact is illustrated in Figure 2.4. Thus we need some special solution techniques so that the convergence of the solution process is independent of the size of the problem. mum". Dani-32:32 Figure 2.4 A circular domain at various resolutions Chapter 3 Applying Fictitious Domain Methods To Solve Frequency Response Problems in Two- Dimensions. 3.1 Introduction This chapter discusses the solution of two dimensional frequency response problems using fictitious domain techniques. In fi'equency response analysis, the goal is to determine the behavior of a structure subjected to a periodic forcing function with specified boundary conditions or constraints. The approach discussed in the following sections uses a fictitious domain to recast a given frequency response problem in elasticity with an arbitrary domain into a periodic problem with a simple, regular domain. Some examples are shown which illustrate the approach. 3.2 Field Equations in Elastodynamics In this case the equilibrium equations are given as UN +f, = pi}, (3.1) The strain-displacement relations and the stress-strain relationships are the same as those given in equations 2.2 and 2.3. Here, p is the mass-density. Combining equation 3.1 with equations 2.2 and 2.3 we get Duhuw, + f, = pit, (3.2) Multiplying this equation with a weight function (w) and integrating over the volume we get the weak form, which is jpwiidV+jDe(u)s(w)dr/= jr’de for all w, eH'(V) (3.3) We account for damping forces by treating them as a special type of non- conservative body force in the above equation [10]. We assume that the damping force at any point is proportional to the velocity and opposite in direction to the velocity, or f, = —Cu, (3.4) where C is a constant. Now, substituting this in equation 3.3, we get IpwiidV+ICwudV+[De(u)s(w)dV=[1”de forallwleH'W) (3.5) V V V V Since the above equation holds for all w, e H l(V), the coemcients of w, must be equal to zero, thus, in matrix form, mii+cu+ku=f (3.6) Where k, m and c are the stiffness, mass and damping matrices respectively. We are concerned only with the steady state solutions to the above equations. We assume solutions of the form u, = v 1e" (3.7) Using the above approximation the discretized equations (3.6) can be written as [k-w2m+iwc]v=f (3.8) 3.3 Problem Statement Figure 3.1. Frequency response problem on an arbitrary domain Figure 3.1 illustrates a typical frequency response problem in elasticity. In the figure F is the section of the boundary where the displacement constraints are applied, f(a)) is the periodic forcing function with to being the forcing frequency. The discretized form of the above problem can be mathematically expressed as, [k-w’mh’we] u=f (3.9) such that u, = 0 on F Here It is the stifl'ness matrix, In is the mass matrix, c is the damping matrix, f is the force vector, a) is the forcing frequency and u is the complex response vector, which contains the magnitude and phase of the displacement. To solve problem (3.9) using the fictitious domain approach, this arbitrary domain is embedded inside a regular (square) fictitious domain as shown in Figure 3.2. Figure 3.2. Frequency response problem on 9 l7 Here 0, e 912 is the original arbitrary domain that is embedded in a square fictitious domain (2, e ‘Rz to get a new domain (1:0, U02, with Q, on, =¢. As before the new domain acts as the domain for a new Q-periodic problem stated as, Az=F (3.10) such that zi = 0 on F and z is periodic where, A = [K -w2 M + iwC], K, M and C are the stiffness, mass and damping matrices of the new heterogeneous domain (1, respectively, and z is the new complex response vector. The constraints in the original problem are enforced in the new fictitious domain problem also. 3.4 Problem Discretization The domain 0 is treated as a square of dimensions NxN where N=2) for a fixed integer j that determines the level of refinement in the approximation. We then designate each unit square X (It, I)=[k, k+l) x [1, 1+1) in Q as a pixel. This efl'ectively yields a level j, NxN raster description of the domain. The material tensor in two-dimensional elasticity has the form d”(x,y) d”(x,y) d”(x,y) D(x.y)= d2'(x.y) d”(x,y) d”(x.y) (3.11) d3'(x,y) d’2(x,y) d”(x,y) We simplify the material model by assuming that the material tensor at any location is represented by a scaled version of the fixed isotropic material tensor D°, where 0" has the form F1 v o ‘ 0 = __E__ " 1 0 (3.12) l — v2 _ ( ) 0 0 l V _ 2 i and E and v are the Youngs modulus and Poisson’s ratio respectively. The material distribution in Q is now expressed as D(x.y) = Way) D° (3-13) where , 6 Way) = { W“ x y 0‘ andw e (0.11 (3.14) W(x9y) x,ye 02 As discussed earlier, with reasonable choices of the fictitious material tensor, the body force and the constraints set, the solution of the fictitious domain problem (3.10) approaches the solution of the original problem (3.9) in the domain of interest. In addition to the suggested choice of material properties, for the material in the fictitious domain to have little or no effect on the solution in Q, (see equations 2.13-2.15), we also make the ratio of the stifi'ness tensors to the mass-densities in the fictitious domain much larger compared to those in the actual problem domain. This is done so that the eigenvalues associated with the modes of deformation in the fictitious domain are very large when compared to those in the domain of interest, that is, Low»—l—Dnz (3.15) peak Po, A (diagonal) lumped mass matrix is assumed. The damping is assumed to be proportional or Raleigh damping, that is C = 0M + flK (3.16) where a and ,B are constants. 3.5 Solution Strategy In the large-scale problems of interest the discretized system of equations is usually very large and may not be solved by direct methods. We consider a conjugate gradient type method for the solution of the sparse linear system in equation (3.10) with complex symmetric coefficient matrices. The boundary conditions are applied as in the standard finite element approach, that is, by zeroing out the rows and columns corresponding to the constrained degree of freedom. Mathematically this can be expressed as, A A, = F (3.17) where, A = [L +(l-L)A(l -L)], I is the identity matrix and L is such that, L". =1 for d.o.f i in F (fixed) (3.18) Ly = 0 otherwise It is often very enticing to avoid solving complex linear systems by converting them into real systems instead. There are some theoretical and numerical results presented by Freund [11] which show that this is a ‘fatal’ approach, at 20 least for Krylov subspace methods. In most cases the resulting real systems are harder to solve by the conjugate gradient-type algorithms than the original complex ones. Also the convergence behavior in this case is very erratic. In this work, one of the simplified quasi-minimal residue (QMR) algorithms for linear systems with J-symmetric or J-Hermitian coefficient matrices has been implemented. This algorithm is one of the simplest possible implementations of the QMR method without look-ahead [12]. This is the quasi-minimal residue from bi-conjugate gradient (QMR-from-BCG) algorithm proposed by Freund and Szeto [13]. By adding updates for the QMR iterates and the QMR search directions and some minor scalar operations to these BCG recursions, they have obtained an algorithm that is only marginally more expensive than BCG and requires on additional matrix- vector product per iteration. 3.6 Preconditioning The conjugate gradient-type algorithms require an effective preconditioner for rapid convergence. The preconditioner (P) is a matrix that in some sense approximates the original matrix (A) and is easy to invert. This is because the preconditioned conjugate gradient-type algorithms require the inversion of the preconditioner at each step. 21 Then, for example, the following preconditioned system P"Az = P"F (3.19) could be solved instead of equation (3.17). In general, this system is no longer symmetric. However, we observe that P“'A for the P-inner product (x, y), has the property, (Ly), -=- (Pr. y) = (x. Py) (320) since (P"Axiy)i- = (Ax. y) = (X. Ay) = (xiP(P"A)y) = (x. F'Py). (3.21) Thus, a simple way of preserving symmetry is to replace the usual Euclidean inner product in the algorithm by the P-inner product [14]. 3.7 Derivation of the Preconditioner Consider the case where the material distribution in Q is homogeneous. In this situation, the coefiicient matrix ( K ,, —a)2M H +iwCH), denoted A“ , possess some special properties. Each of its sub-matrices in its partitioned structure is an Ni’xN2 block-circulant matrix A - A“ A"" (3 22) H An, Ari». . A circulant matrix is a square matrix that can be generated from a single row of data. The (i+l)"' row of a circulant matrix is generated by performing a circulant shift of the 1'" row once to the right. An N 2xN 2 block-circulant 22 matrix is made up of N2 blocks of N xN circulant matrices arranged in a circulant fashion, i.e., the (i+l)”' block row is generated by shifting the in block row once to the right and wrapping around the right end. This is illustrated in Figure 3.3. I II III IV ‘4 \1234 4l23 III IV I II 34l2 234] III IV I ll II III IV I Figure 3.3. Structure of block-circulant matrices Block-circulant matrices can be stored very efficiently (only one row is needed) and inverted very efficiently using Fourier Transforms [15]. This is due to the fact that block-circulant matrices are diagonal in the Fourier basis. Recalling that each sub-matrix of A" is block circulant, A" is tri-diagonal when represented in the Fourier space. For our particular case a good candidate for the preconditioning matrix (P) would be p = A" +xL (3.23) 23 where, L is as described in equation (3.18) and x is a large penalty parameter. Here, the boundary conditions are applied in the form of a penalty; this is so as to make the inversion of the preconditioner possible as otherwise the matrix would be singular. 3.8 Inversion of the Preconditioner The inversion of the preconditioning matrix (P) can be performed very efficiently using a two-step procedure that exploits the fact that A" is a symmetric, block-circulant matrix. First P" is expressed as P" = We +pp'>+xB'B} -pp‘r' (3.24) where BTB=L 1...1 o...o ' (3.25) I”[0...0 1m] Here p is the matrix of rigid body modes of the structure. The matrix pp'l' is added to A11 to remove the rigid body mode from the Q-periodic solution. Matrix (Au +ppT) is also symmetric and block-circulant and thus it can be inverted very efficiently using Fourier 'h'ansforms. The subtraction of pp” and the addition of 878 are accomplished using the Sherman-Morrison- Woodbury formulae for the matrix inversion of rank-r updates [16] 24 Given the inverse of a square matrix A, the Sherman-Morrison formula gives us a way of computing the inverse for changes in A of the form A —> (A+ n ® v) (3.26) for some vector u and v. (A+u®v)" = A'I _(A"u)®(vA") (3.27) 1+ 2. where A 5 vA"u (3.28) In our case, we want to add more than a single correction term, then the Sherman-Morrison formula cannot be directly applied for the simple reason that the storage of the whole inverse matrix A“' is not feasible. Instead, we use the Woodbury formula, which is the block matrix form of the Sherman- Morrison formula. (A + IN)" = A" —[A"'U(l + v’A"U)"v“A"] (3.29) Here A is an nxnmatrix, while U and V are nx p matrices with p