# Numerical Investigation of Aeroelastic Mode Distribution for Aircraft Wing Model in Subsonic Air Flow

Authors: Marianna A. Shubov; Stephen Wineberg; Robert Holt
Journal: Mathematical Problems in Engineering
 (2010)

---

## Abstract

In this paper, the numerical results on two problems originated in aircraft wing modeling have been presented.The first problem is concerned with the approximation to the set of the aeroelastic modes, which are the eigenvalues of a certain boundary-value problem. The affirmative answer is given to the following question: can the leading asymptotical terms in the analytical formulas be used as reasonably accurate description of the aeroelastic modes? The positive answer means that these leading terms can be used by engineers for practical calculations. The second problem is concerned with the flutter phenomena in aircraft wings in a subsonic, incompressible, inviscid air flow. It has been shown numerically that there exists a pair of the aeroelastic modes whose behavior depends on a speed of an air flow. Namely, when the speed increases, the distance between the modes tends to zero, and at some speed that can be treated as the flutter speed these two modes merge into one double mode.


---

## Body

## 1. Introduction and Formulation of Problem

The main goal of this paper is to present the numerical results for two problems arising in aircraft wing modeling.The first one is related to the numerical approximation of the discrete spectrum of a certain matrix differential operator that represents the structural part of an aircraft wing model in a subsonic incompressible inviscid air flow. More precisely, it is shown that the asymptotic approximation for the set of aeroelastic modes has a very regular structure, that is, there are two distinct branches which are called the β-branch and the δ-branch (see formulas (2.2), (2.3)). Each branch is described by leading asymptotical terms and remainder terms. An affirmative answer is given to the following question: can the leading asymptotical terms be used as a reasonably accurate description of the aeroelastic modes? The question can be rephrased as: if one discretizes the main matrix differential operator with spectral methods, can the spectrum of the discretized finite-dimensional operator be described by the leading asymptotical terms of the continuous operator? The answer is affirmative, and one can see almost perfect spectral branches for the finite-dimensional approximation (see Figure 2).The second problem is related to the entire model involving both the structural and aerodynamical parts. In the second problem, the authors provide a numerical justification of a certain fact well-known in the aircraft engineering community. Namely, it is shown that there exists a pair of the aeroelastic modes whose behavior depends on a speed of an air flow in the following way: the distance between these aeroelastic modes tends to zero with increasing speed. At some speed that can be treated as “the flutter speed”, the two modes merge in one double mode. If a speed continues to increase, the double aeroelastic mode splits up into two simple modes that are moving apart. Such a dynamics clearly indicates that there exists a specific speed (“the flutter speed”) that causes an appearance of a double aeroelastic mode leading to a flutter instability. In this study, the authors deal with a linear wing model, which is obviously a significant simplification of a realistic situation. In a more precise setting, when the speed is near what is called “the flutter speed”, one should use a nonlinear wing model, which results in limit cycle oscillations rather than flutter. One of the important conclusions of this investigation is that even though the model studied in the paper is linear, it is capable of capturing the flutter instability, which is definitely supported by numerical simulations. Before the numerical analysis results, a mathematical statement of the problem and a description of those analytical results that are relevant to the numerical analysis will be given.An ultimate goal of an aircraft wing modeling is to find an approach to flutter control. Flutter is a structural dynamical instability, which consists of violent vibrations of the solid structure with rapidly increasing amplitude [1, 2]. The physical reason for this phenomenon is that under special conditions, the energy of air flow is rapidly absorbed by the structure and transformed into the energy of mechanical vibrations. In engineering practice, flutter must be avoided either by design of the structure or by introducing control mechanism capable to suppress harmful vibrations. Flutter is an inherent feature of fluid-structure interaction and, thus, it cannot be eliminated completely. However, the critical conditions for the flutter onset can be shifted to the safe range of the operating parameters. This is exactly the goal for designing flutter control mechanisms.Flutter is an in-flight event, which happens beyond some speed-altitude combinations. High-speed aircrafts are most susceptible for flutter although no speed regime is truly immune from flutter. Flutter instabilities occur in a variety of different engineering and even biomedical situations. Namely in aeronautic engineering, flutter of helicopter, propeller, and turbine blades is a serious problem. It also affects electric transmission lines, high-speed hard disk drives, and long-spanned suspension bridges. Flutter of cardiac tissue and blood vessel walls is of a special concern to medical practitioners.At the present moment, there exist only a few models of fluid-structure interaction involving flutter development, for which precise mathematical formulations are available. The authors' main objects of interest are analytical and numerical results on such models, which can be used, in particular, for flutter explanation. It is certainly important for designing flutter control mechanisms.Ideally, a complete picture of a fluid-structure interaction should be described by a system of partial differential equations, a system which contains both the equations governing the vibrations of an elastic structure and the hydrodynamic equations governing the motion of gas or fluid flow. The system of equations of motion should be supplied with appropriate boundary and initial conditions. The structural and hydrodynamic parts of the system must be coupled in the following sense. The hydrodynamic equations define a pressure distribution on the elastic structure. This pressure distribution in turn defines the so-called aerodynamic loads, which appear as forcing terms in structural equations. On the other hand, the parameters of the elastic structure enter the boundary conditions for the hydrodynamic equations.In the present paper it is assumed that the model describes a wing of high-aspect ratio (i.e., the length of a wing is much greater than its width, though both quantities are finite) in asubsonic, inviscid, incompressible air flow. The hydrodynamic equations have been solved explicitly and aerodynamic loads are represented as forcing terms in the structural equations as time convolution-type integrals with complicated kernels. Thus, the model is described by a system of integro-differential equations. Analytical results on this model include the following. The system of equations of motion is treated as a single evolution-convolution equation in the Hilbert state space of the model. (The integral convolution part of this equation vanishes if a speed of an air stream is zero, and the equations of motion describe the so-called ground vibrations.) After applying a Laplace transformation with respect to the time variable to both sides of the evolution equation, one obtains a matrix differential equation involving the complex parameter λ. For this new equation, the following results have been shown.(a) The representation of the solution of the original initial boundary-value problem in the frequency domain has been given in terms ofthe generalized resolvent, which is an analytic operator-valued function of the spectral parameter λ. The aeroelastic modes are defined as the poles of the generalized resolvent; the corresponding mode shapes are defined in terms of the residues at these poles [3–5].(b) Explicit asymptotic formulas for the aeroelastic modes have been derived [4, 6]. (To the best of the authors' knowledge, these are the first such formulas in the literature on aeroelasticity.) The entire set of aeroelastic modes splits into two branches, which are asymptotically close to the eigenvalues of the structural part of the system [3, 5, 7].(c) It has been shown that the set of the mode shapes forms a nonorthogonal basis (Riesz basis) of the state space of the system. (The set of the generalized eigenvectors of the structural part of the system has a similar property [4, 8].)


### 1.1. Statement of Problem. Operator Setting in Energy Space

The so-called “Goland model” [1, 6, 9] is considered, that is, the simplest structural model—a uniform rectangular beam (Figure 1) with two types of motion, plunge and pitch. To formulate a mathematical model, the following dynamical variables  are introduced [3, 4, 9, 10]:(1.1)X(x,t)=(h(x,t)α(x,t)),-L≤x≤0,t≥0,
where h(x,t) is the bending, α(x,t) is the torsion angle, and x is the span variable. The model can be described by the following linear system:(1.2)(Ms-Ma)Ẍ(x,t)+(Ds-uDa)Ẋ(x,t)+(Ks-u2Ka)X(x,t)=[f1(x,t)f2(x,t)],
where the “overdot” denotes the differentiation with respect to time t. The subscripts “s” and “a” are use to distinguish the structural and aerodynamical parameters, respectively. All 2×2 matrices in (1.2) are given by the following formulas:(1.3)Ms=[mSSI],Ma=(-πρ)[1-a-a(a2+18)],Ds=[0000],
where m is the density of the flexible structure, S is the mass moment, I is the moment of inertia, ρ is the density of air, a is the linear parameter of the structure, and -1≤a≤1 (a is a relative distance between the elastic axis of a model wing and its line of center of gravity). Even though in the current paper the matrix Ds has only trivial entries, it is preferable to keep it in (1.2) since the problem with a nontrivial structural damping Ds will be studied in the authors' forthcoming paper:(1.4)Da=(-πρ)[01-10],Ks=[E∂4∂x400-G∂2∂x2],Ka=(-πρ)[000-1],
where E is the bending stiffness, and G is the torsion stiffness. The parameter u in (1.2) denotes the stream speed. The right hand side of system (1.2) can be represented as the following system of two convolution-type integral operations:(1.5)f1(x,t)=-2πρ∫0t[uC2(t-σ)-Ċ3(t-σ)]g(x,σ)dσ≡∫0tC̃1(t-σ)g(x,σ)dσ,(1.6)f2(x,t)=-2πρ∫0t[12C1(t-σ)-auC2(t-σ)+aĊ3(t-σ)+uC4(t-σ)+12Ċ5(t-σ)]g(x,σ)dσ≡∫0tC̃2(t-σ)g(x,σ)dσ,(1.7)g(x,t)=uα̇(x,t)+ḧ(x,t)+(12-a)α̈(x,t).
Explicit formulas for the aerodynamical functions Ci, i=1,…,5, can be found in [11]. It is known that the self-straining control actuator action can be modeled by the following boundary conditions [5, 9, 12–14]:(1.8)Eh′′(0,t)+βh′̇(0,t)=0,h′′′(0,t)=0,Gα′(0,t)+δα̇(0,t)=0,β,δ∈ℂ+∪{∞},
where ℂ+ is the closed right half-plane. The boundary conditions at x=-L are (1.9)h(-L,t)=h′(-L,t)=α(-L,t)=0.
In (1.8) and (1.9) and below, the prime designates derivative with respect to x. The initial state of the system is given as follows:(1.10)h(x,0)=h0(x),ḣ(x,0)=h1(x),α(x,0)=α0(x),α̇(x,0)=α1(x).
The solution of the problem given by (1.2) and conditions (1.8)–(1.10) are given in the energy space ℋ. It is assumed that the parameters satisfy the following two conditions:(1.11)det[mSSI]>0,0<u≤2GLπρ.
The second condition in (1.11) means that the flow speed must be below the “divergence” or static aeroelastic instability speed. The state space ℋ, which is a Hilbert space, is described as the set of 4-component vector-valued functions Ψ=(h,ḣ,α,α̇)T≡(ψ0,ψ1,ψ2,ψ3)T (“T” stands for transposition) obtained as a closure of smooth functions satisfying the boundary conditions  (1.12)ψ0(-L)=ψ0′(-L)=ψ2(-L)=0
in the following energy norm, which is well-defined under conditions (1.11):(1.13)∥Ψ∥ℋ2=12∫-L0[(ψ3(x)ψ̅1(x)+ψ̅3(x)ψ1(x))E|ψ0′′(x)|2+G|ψ2′(x)|2+m̃|ψ1(x)|2+Ĩ|ψ3(x)|2+S̃(ψ3(x)ψ̅1(x)+ψ̅3(x)ψ1(x))-πρu2|ψ2(x)|2]dx,
where the following notations have been used:(1.14)m̃=m+πρ,S̃=S-aπρ,Ĩ=I+πρ(a2+18),Δ=m̃Ĩ-S̃2.
The initial-boundary value problem defined by (1.2) and conditions (1.8)–(1.10) can be represented in the following form  (1.15)Ψ̇=iℒβδΨ+ℱ̃Ψ̇,Ψ=(ψ0,ψ1,ψ2,ψ3)T,Ψ|t=0=Ψ0.ℒβδ is a matrix differential operator in ℋ given by the following expression:(1.16)ℒβδ=-i[0100-EĨΔd4dx4-πρuS̃Δ-S̃Δ(Gd2dx2+πρu2)-πρuĨΔ0001ES̃Δd4dx4πρum̃Δm̃Δ(Gd2dx2+πρu2)πρuS̃Δ]
defined on the domain  (1.17)𝒟(ℒβδ)={Ψ∈ℋ:ψ0∈H4(-L,0),ψ1∈H2(-L,0),ψ2∈H2(-L,0),ψ3∈H1(-L,0);ψ1(-L)=ψ1′(-L)=ψ3(-L)=0;ψ0′′′(0)=0;Eψ0′′(0)+βψ1′(0)=0,Gψ2′(0)+δψ3(0)=0H4(-L,0)}.ℱ̃ is a linear integral operator in ℋ given by the formula (1.18)ℱ̃=1Δ[10000[Ĩ(C̃1*)-S̃(C̃2*)]000010000[-S̃(C̃1*)+m̃(C̃2*)]][000001u(12-a)000001u(12-a)].
In (1.18), the star “*” stands for the convolution operation and the kernels C̃1 and C̃2 are defined in (1.5), and (1.6).Figure 1
Wing structure beam model.Figure 2
Distinctβ-and δ-branches.Remark 1.1.
It is important to emphasize that (1.15) is not an evolution equation. It does not have a dynamics generator and does not define any semigroup in the standard sense. By applying the Laplace transformation to both sides of (1.15), one obtains the following Laplace transform representation of the solution:
(1.19)Ψ̂(λ)=-(λI-iℒβδ-λ𝔉̂(λ))-1(I-𝔉̂(λ))Ψ0.
To find the solution in the space-time domain, one has to “calculate” the inverse Laplace transform of Ψ̂ by contour integration in the complex λ-plane. In this connection, the properties of the “generalized resolvent operator”
(1.20)ℛ(λ)=(λI-iℒβδ-λ𝔉̂(λ))-1
are of crucial importance. It has been proved that ℛ(λ) has a countable set of poles which we called the eigenvalues or the aeroelastic modes. The residues of ℛ(λ) at these poles are precisely the projectors on the corresponding generalized eigenspaces. The differential part and the role of the convolution part of the system have been analyzed in [3–8, 15]. In particular, it has been shown that the convolution part does not “destroy” the main characteristics of the discrete spectrum, which is produced by the differential part of the problem. Namely, it has been proved that the aeroelastic modes are asymptotically close to the discrete spectrum of the operator iℒβδ, and the rate at which the aeroelastic modes approach that spectrum has been calculated.


## 1.1. Statement of Problem. Operator Setting in Energy Space

The so-called “Goland model” [1, 6, 9] is considered, that is, the simplest structural model—a uniform rectangular beam (Figure 1) with two types of motion, plunge and pitch. To formulate a mathematical model, the following dynamical variables  are introduced [3, 4, 9, 10]:(1.1)X(x,t)=(h(x,t)α(x,t)),-L≤x≤0,t≥0,
where h(x,t) is the bending, α(x,t) is the torsion angle, and x is the span variable. The model can be described by the following linear system:(1.2)(Ms-Ma)Ẍ(x,t)+(Ds-uDa)Ẋ(x,t)+(Ks-u2Ka)X(x,t)=[f1(x,t)f2(x,t)],
where the “overdot” denotes the differentiation with respect to time t. The subscripts “s” and “a” are use to distinguish the structural and aerodynamical parameters, respectively. All 2×2 matrices in (1.2) are given by the following formulas:(1.3)Ms=[mSSI],Ma=(-πρ)[1-a-a(a2+18)],Ds=[0000],
where m is the density of the flexible structure, S is the mass moment, I is the moment of inertia, ρ is the density of air, a is the linear parameter of the structure, and -1≤a≤1 (a is a relative distance between the elastic axis of a model wing and its line of center of gravity). Even though in the current paper the matrix Ds has only trivial entries, it is preferable to keep it in (1.2) since the problem with a nontrivial structural damping Ds will be studied in the authors' forthcoming paper:(1.4)Da=(-πρ)[01-10],Ks=[E∂4∂x400-G∂2∂x2],Ka=(-πρ)[000-1],
where E is the bending stiffness, and G is the torsion stiffness. The parameter u in (1.2) denotes the stream speed. The right hand side of system (1.2) can be represented as the following system of two convolution-type integral operations:(1.5)f1(x,t)=-2πρ∫0t[uC2(t-σ)-Ċ3(t-σ)]g(x,σ)dσ≡∫0tC̃1(t-σ)g(x,σ)dσ,(1.6)f2(x,t)=-2πρ∫0t[12C1(t-σ)-auC2(t-σ)+aĊ3(t-σ)+uC4(t-σ)+12Ċ5(t-σ)]g(x,σ)dσ≡∫0tC̃2(t-σ)g(x,σ)dσ,(1.7)g(x,t)=uα̇(x,t)+ḧ(x,t)+(12-a)α̈(x,t).
Explicit formulas for the aerodynamical functions Ci, i=1,…,5, can be found in [11]. It is known that the self-straining control actuator action can be modeled by the following boundary conditions [5, 9, 12–14]:(1.8)Eh′′(0,t)+βh′̇(0,t)=0,h′′′(0,t)=0,Gα′(0,t)+δα̇(0,t)=0,β,δ∈ℂ+∪{∞},
where ℂ+ is the closed right half-plane. The boundary conditions at x=-L are (1.9)h(-L,t)=h′(-L,t)=α(-L,t)=0.
In (1.8) and (1.9) and below, the prime designates derivative with respect to x. The initial state of the system is given as follows:(1.10)h(x,0)=h0(x),ḣ(x,0)=h1(x),α(x,0)=α0(x),α̇(x,0)=α1(x).
The solution of the problem given by (1.2) and conditions (1.8)–(1.10) are given in the energy space ℋ. It is assumed that the parameters satisfy the following two conditions:(1.11)det[mSSI]>0,0<u≤2GLπρ.
The second condition in (1.11) means that the flow speed must be below the “divergence” or static aeroelastic instability speed. The state space ℋ, which is a Hilbert space, is described as the set of 4-component vector-valued functions Ψ=(h,ḣ,α,α̇)T≡(ψ0,ψ1,ψ2,ψ3)T (“T” stands for transposition) obtained as a closure of smooth functions satisfying the boundary conditions  (1.12)ψ0(-L)=ψ0′(-L)=ψ2(-L)=0
in the following energy norm, which is well-defined under conditions (1.11):(1.13)∥Ψ∥ℋ2=12∫-L0[(ψ3(x)ψ̅1(x)+ψ̅3(x)ψ1(x))E|ψ0′′(x)|2+G|ψ2′(x)|2+m̃|ψ1(x)|2+Ĩ|ψ3(x)|2+S̃(ψ3(x)ψ̅1(x)+ψ̅3(x)ψ1(x))-πρu2|ψ2(x)|2]dx,
where the following notations have been used:(1.14)m̃=m+πρ,S̃=S-aπρ,Ĩ=I+πρ(a2+18),Δ=m̃Ĩ-S̃2.
The initial-boundary value problem defined by (1.2) and conditions (1.8)–(1.10) can be represented in the following form  (1.15)Ψ̇=iℒβδΨ+ℱ̃Ψ̇,Ψ=(ψ0,ψ1,ψ2,ψ3)T,Ψ|t=0=Ψ0.ℒβδ is a matrix differential operator in ℋ given by the following expression:(1.16)ℒβδ=-i[0100-EĨΔd4dx4-πρuS̃Δ-S̃Δ(Gd2dx2+πρu2)-πρuĨΔ0001ES̃Δd4dx4πρum̃Δm̃Δ(Gd2dx2+πρu2)πρuS̃Δ]
defined on the domain  (1.17)𝒟(ℒβδ)={Ψ∈ℋ:ψ0∈H4(-L,0),ψ1∈H2(-L,0),ψ2∈H2(-L,0),ψ3∈H1(-L,0);ψ1(-L)=ψ1′(-L)=ψ3(-L)=0;ψ0′′′(0)=0;Eψ0′′(0)+βψ1′(0)=0,Gψ2′(0)+δψ3(0)=0H4(-L,0)}.ℱ̃ is a linear integral operator in ℋ given by the formula (1.18)ℱ̃=1Δ[10000[Ĩ(C̃1*)-S̃(C̃2*)]000010000[-S̃(C̃1*)+m̃(C̃2*)]][000001u(12-a)000001u(12-a)].
In (1.18), the star “*” stands for the convolution operation and the kernels C̃1 and C̃2 are defined in (1.5), and (1.6).Figure 1
Wing structure beam model.Figure 2
Distinctβ-and δ-branches.Remark 1.1.
It is important to emphasize that (1.15) is not an evolution equation. It does not have a dynamics generator and does not define any semigroup in the standard sense. By applying the Laplace transformation to both sides of (1.15), one obtains the following Laplace transform representation of the solution:
(1.19)Ψ̂(λ)=-(λI-iℒβδ-λ𝔉̂(λ))-1(I-𝔉̂(λ))Ψ0.
To find the solution in the space-time domain, one has to “calculate” the inverse Laplace transform of Ψ̂ by contour integration in the complex λ-plane. In this connection, the properties of the “generalized resolvent operator”
(1.20)ℛ(λ)=(λI-iℒβδ-λ𝔉̂(λ))-1
are of crucial importance. It has been proved that ℛ(λ) has a countable set of poles which we called the eigenvalues or the aeroelastic modes. The residues of ℛ(λ) at these poles are precisely the projectors on the corresponding generalized eigenspaces. The differential part and the role of the convolution part of the system have been analyzed in [3–8, 15]. In particular, it has been shown that the convolution part does not “destroy” the main characteristics of the discrete spectrum, which is produced by the differential part of the problem. Namely, it has been proved that the aeroelastic modes are asymptotically close to the discrete spectrum of the operator iℒβδ, and the rate at which the aeroelastic modes approach that spectrum has been calculated.


## 2. Asymptotic and Spectral Results on Matrix Differential and Integral Operators

### 2.1. Matrix Differential Operatorℒβδ

Theorem 2.1.
(a)ℒβδ is a closed linear operator in ℋ, whose resolvent is compact, and therefore, the spectrum is discrete [3, 4, 16].
(b) Operatorℒβδ is nonselfadjoint unless ℜβ=ℜδ=0. If ℜβ≥0 and ℜδ≥0, then this operator is dissipative, that is, ℑ(ℒβδΨ,Ψ)≥0 for Ψ∈𝒟(ℒβδ). The adjoint operator ℒβδ* is given by the matrix differential expression (1.16) on the domain obtained from (1.17) by replacing the parameters β and δ with (-β̅) and (-δ̅), respectively.Theorem 2.2.
(a) The operatorℒβδ has a countable set of complex eigenvalues. If
(2.1)δ≠GĨ,
then the set of eigenvalues is located in a strip parallel to the real axis.
(b) The entire set of eigenvalues asymptotically splits into two different subsets. We call them theβ-branch and the δ-branch and denote these branches by {νnβ}n∈ℤ and {νnδ}n∈ℤ, respectively. If ℜβ≥0 and ℜδ>0, then each branch is asymptotically close to its own horizontal line in the closed upper half-plane. If ℜβ>0 and ℜδ=0, then both horizontal lines coincide with the real axis. If ℜβ=ℜδ=0, then the operator ℒβδ is selfadjoint and, thus, its spectrum is real.
(c) The following asymptotics are valid for theβ-branch of the spectrum as |n|→∞:(2.2)νnβ=(sgnn)(π2L2)EĨΔ(|n|-14)2+κn(ω),ω=|δ|-1+|β|-1,
with Δ being defined in (1.14). A complex-valued sequence {κn} is bounded above in the following sense: supn∈ℤ{|κn(ω)|}=C(ω),C(ω)→0asω→0.
(d) The following asymptotics are valid for theδ-branch of the spectrum:(2.3)νnδ=πnLĨ/G+i2LĨ/Glnδ+GĨδ-GĨ+O(|n|-1/2),|n|→∞.
In (2.3), “ln” means the principal value of the logarithm. If β and δ stay away from zero, that is, |β|≥β0>0 and |δ|≥δ0>0, then the estimate O(|n|-1/2) in (2.3) is uniform with respect to both parameters.
(e) There may be only a finite number of multiple eigenvalues of a finite multiplicity each. Therefore, only a finite number of the associate vectors may exist.The theorem below deals with the Riesz basis property of the generalized eigenvectors of the structural operatorℒβδ. The Riesz basis is a mild modification of an orthonormal basis, namely a linear isomorphism of an orthonormal basis.Theorem 2.3.
The set of generalized eigenvectors of the operatorℒβδ forms a Riesz basis in the energy space ℋ.The next statement deals with the asymptotic distribution of the aeroelastic modes [4, 7, 15].Theorem 2.4.
(a) The set of the aeroelastic modes (which are the poles of the generalized resolvent operator) is countable and does not have accumulation points on the complex planeℂ. There might be only a finite number of multiple poles of a finite multiplicity each. There exists a sufficiently large R>0 such that all aeroelastic modes, whose distance from the origin is greater than R, are simple poles of the generalized resolvent. The value of R depends on the speed u of an airstream, that is, R=R(u).
(b) The set of the aeroelastic modes splits asymptotically into two series, which we call theβ-branch and the δ-branch. Asymptotical distribution of the β-and the δ-branches of the aeroelastic modes can be obtained from asymptotical distribution of the spectrum of the operator ℒβδ. Namely if {λnβ}n∈ℤ is the β-branch of the aeroelastic modes, then λnβ=iλ̂nβ and the asymptotics of the set {λ̂nβ}n∈ℤ is given by the right-hand side of formula (2.2). Similarly, if {λnδ=iλ̂nδ}n∈ℤ is the δ-branch of the aeroelastic modes, then the asymptotical distribution of the set {λ̂nδ}n∈ℤ is given by the right-hand side of formula (2.3).


### 2.2. Structure and Properties of the Matrix Integral Operator

Important information on the Laplace transform of the convolution-type matrix integral operator (1.18) is collected below.Lemma 2.5.
Let𝔉̂ be the Laplace transform of the kernel of matrix integral operator (1.18). The following formula is valid for 𝔉̂:(2.4)𝔉̂(λ)=[00000ℒ(λ)uℒ(12-a)ℒ(λ)00000𝒩(λ)u𝒩(λ)(12-a)𝒩(λ)],
where
(2.5)ℒ(λ)=-2πρΔuλ{-S̃2+[Ĩ+(a+12)S̃]T(λu)},𝒩(λ)=-2πρΔuλ{m̃2-[S̃+(a+12)m̃]T(λu)},
and T is the Theodorsen function defined by the formula
(2.6)T(z)=K1(z)K0(z)+K1(z);K0 and K1 are the modified Bessel functions [17].


### 2.3. Properties of the Theodorsen Function

The Theodorsen function is a bounded analytic function on the complex plane with the branch-cut along the negative real semi-axis. As|z|→∞, the following asymptotic representation holds:(2.7)T(z)=12+O(11+|z|),as|z|→∞.
Taking into account that z=λ/u, one can write λ𝔉̂(λ) as the following sum: λ𝔉̂(λ)=𝔐+𝔑(λ), where the matrix 𝔐 is defined by the formula  (2.8)𝔐=[00000AuA(12-a)A00000BuB(12-a)B]
with A and B being given by (2.9)A=-πρuΔ-1[Ĩ+(a-12)S̃],B=πρuΔ-1[S̃+(a-12)m̃].
The matrix-valued function 𝔑(λ) is defined by the formula  (2.10)𝔑(λ)=[00000A1(λ)uA1(λ)(12-a)A1(λ)00000B1(λ)uB1(λ)(12-a)B1(λ)],
where A1(λ)=-2πρuΔ-1(T(λ/u)-1/2)[Ĩ+(a+1/2)S̃], and B1(λ)=2πρuΔ-1(T(λ/u)-1/2)[S̃+(a+1/2)m̃]. For each λ, 𝔑(λ) is a bounded operator in ℋ with the following estimate for its norm:(2.11)∥𝔑∥ℋ≤C(1+|λ|)-1,
where C is an absolute constant the precise value of which is immaterial for us. Therefore, the generalized resolvent (1.20) can be written in the form  (2.12)ℛ(λ)=𝔖-1(λ),where𝔖(λ)=λI-iℒβδ-𝔐-𝔑(λ).Theorem 2.6.
𝔐 is a bounded linear operator in  ℋ. The operator  𝒦βδ defined by
(2.13)𝒦βδ=ℒβδ-i𝔐
is an unbounded nonselfadjoint operator in ℋ with compact resolvent. The spectral asymptotics of 𝒦βδ coincide with the spectral asymptotics of ℒβδ (see Theorem 2.2). In contrast to ℒβδ, the operator 𝒦βδ is not dissipative for any boundary control gains. However, 𝒦βδ is also a Riesz spectral operator, that is, the set of its generalized eigenvectors forms a Riesz basis of ℋ.


## 2.1. Matrix Differential Operatorℒβδ

Theorem 2.1.
(a)ℒβδ is a closed linear operator in ℋ, whose resolvent is compact, and therefore, the spectrum is discrete [3, 4, 16].
(b) Operatorℒβδ is nonselfadjoint unless ℜβ=ℜδ=0. If ℜβ≥0 and ℜδ≥0, then this operator is dissipative, that is, ℑ(ℒβδΨ,Ψ)≥0 for Ψ∈𝒟(ℒβδ). The adjoint operator ℒβδ* is given by the matrix differential expression (1.16) on the domain obtained from (1.17) by replacing the parameters β and δ with (-β̅) and (-δ̅), respectively.Theorem 2.2.
(a) The operatorℒβδ has a countable set of complex eigenvalues. If
(2.1)δ≠GĨ,
then the set of eigenvalues is located in a strip parallel to the real axis.
(b) The entire set of eigenvalues asymptotically splits into two different subsets. We call them theβ-branch and the δ-branch and denote these branches by {νnβ}n∈ℤ and {νnδ}n∈ℤ, respectively. If ℜβ≥0 and ℜδ>0, then each branch is asymptotically close to its own horizontal line in the closed upper half-plane. If ℜβ>0 and ℜδ=0, then both horizontal lines coincide with the real axis. If ℜβ=ℜδ=0, then the operator ℒβδ is selfadjoint and, thus, its spectrum is real.
(c) The following asymptotics are valid for theβ-branch of the spectrum as |n|→∞:(2.2)νnβ=(sgnn)(π2L2)EĨΔ(|n|-14)2+κn(ω),ω=|δ|-1+|β|-1,
with Δ being defined in (1.14). A complex-valued sequence {κn} is bounded above in the following sense: supn∈ℤ{|κn(ω)|}=C(ω),C(ω)→0asω→0.
(d) The following asymptotics are valid for theδ-branch of the spectrum:(2.3)νnδ=πnLĨ/G+i2LĨ/Glnδ+GĨδ-GĨ+O(|n|-1/2),|n|→∞.
In (2.3), “ln” means the principal value of the logarithm. If β and δ stay away from zero, that is, |β|≥β0>0 and |δ|≥δ0>0, then the estimate O(|n|-1/2) in (2.3) is uniform with respect to both parameters.
(e) There may be only a finite number of multiple eigenvalues of a finite multiplicity each. Therefore, only a finite number of the associate vectors may exist.The theorem below deals with the Riesz basis property of the generalized eigenvectors of the structural operatorℒβδ. The Riesz basis is a mild modification of an orthonormal basis, namely a linear isomorphism of an orthonormal basis.Theorem 2.3.
The set of generalized eigenvectors of the operatorℒβδ forms a Riesz basis in the energy space ℋ.The next statement deals with the asymptotic distribution of the aeroelastic modes [4, 7, 15].Theorem 2.4.
(a) The set of the aeroelastic modes (which are the poles of the generalized resolvent operator) is countable and does not have accumulation points on the complex planeℂ. There might be only a finite number of multiple poles of a finite multiplicity each. There exists a sufficiently large R>0 such that all aeroelastic modes, whose distance from the origin is greater than R, are simple poles of the generalized resolvent. The value of R depends on the speed u of an airstream, that is, R=R(u).
(b) The set of the aeroelastic modes splits asymptotically into two series, which we call theβ-branch and the δ-branch. Asymptotical distribution of the β-and the δ-branches of the aeroelastic modes can be obtained from asymptotical distribution of the spectrum of the operator ℒβδ. Namely if {λnβ}n∈ℤ is the β-branch of the aeroelastic modes, then λnβ=iλ̂nβ and the asymptotics of the set {λ̂nβ}n∈ℤ is given by the right-hand side of formula (2.2). Similarly, if {λnδ=iλ̂nδ}n∈ℤ is the δ-branch of the aeroelastic modes, then the asymptotical distribution of the set {λ̂nδ}n∈ℤ is given by the right-hand side of formula (2.3).


## 2.2. Structure and Properties of the Matrix Integral Operator

Important information on the Laplace transform of the convolution-type matrix integral operator (1.18) is collected below.Lemma 2.5.
Let𝔉̂ be the Laplace transform of the kernel of matrix integral operator (1.18). The following formula is valid for 𝔉̂:(2.4)𝔉̂(λ)=[00000ℒ(λ)uℒ(12-a)ℒ(λ)00000𝒩(λ)u𝒩(λ)(12-a)𝒩(λ)],
where
(2.5)ℒ(λ)=-2πρΔuλ{-S̃2+[Ĩ+(a+12)S̃]T(λu)},𝒩(λ)=-2πρΔuλ{m̃2-[S̃+(a+12)m̃]T(λu)},
and T is the Theodorsen function defined by the formula
(2.6)T(z)=K1(z)K0(z)+K1(z);K0 and K1 are the modified Bessel functions [17].


## 2.3. Properties of the Theodorsen Function

The Theodorsen function is a bounded analytic function on the complex plane with the branch-cut along the negative real semi-axis. As|z|→∞, the following asymptotic representation holds:(2.7)T(z)=12+O(11+|z|),as|z|→∞.
Taking into account that z=λ/u, one can write λ𝔉̂(λ) as the following sum: λ𝔉̂(λ)=𝔐+𝔑(λ), where the matrix 𝔐 is defined by the formula  (2.8)𝔐=[00000AuA(12-a)A00000BuB(12-a)B]
with A and B being given by (2.9)A=-πρuΔ-1[Ĩ+(a-12)S̃],B=πρuΔ-1[S̃+(a-12)m̃].
The matrix-valued function 𝔑(λ) is defined by the formula  (2.10)𝔑(λ)=[00000A1(λ)uA1(λ)(12-a)A1(λ)00000B1(λ)uB1(λ)(12-a)B1(λ)],
where A1(λ)=-2πρuΔ-1(T(λ/u)-1/2)[Ĩ+(a+1/2)S̃], and B1(λ)=2πρuΔ-1(T(λ/u)-1/2)[S̃+(a+1/2)m̃]. For each λ, 𝔑(λ) is a bounded operator in ℋ with the following estimate for its norm:(2.11)∥𝔑∥ℋ≤C(1+|λ|)-1,
where C is an absolute constant the precise value of which is immaterial for us. Therefore, the generalized resolvent (1.20) can be written in the form  (2.12)ℛ(λ)=𝔖-1(λ),where𝔖(λ)=λI-iℒβδ-𝔐-𝔑(λ).Theorem 2.6.
𝔐 is a bounded linear operator in  ℋ. The operator  𝒦βδ defined by
(2.13)𝒦βδ=ℒβδ-i𝔐
is an unbounded nonselfadjoint operator in ℋ with compact resolvent. The spectral asymptotics of 𝒦βδ coincide with the spectral asymptotics of ℒβδ (see Theorem 2.2). In contrast to ℒβδ, the operator 𝒦βδ is not dissipative for any boundary control gains. However, 𝒦βδ is also a Riesz spectral operator, that is, the set of its generalized eigenvectors forms a Riesz basis of ℋ.


## 3. Numerical Results for Two-Branch Discrete  Spectrum of Operatorℒβδ

The first question is related to the accuracy of the asymptotic approximations of the eigenvalues of the operatorℒβδ. By their nature, asymptotic formulas (2.2) and (2.3) should be understood in the following way.Formula (2.2) for the β-branch eigenvalues means that there exist a positive number N1 and a small constant 0<ɛ≪1 such that for all |n|≥N1, the β-branch eigenvalues (2.2) satisfy the estimate (3.1)|νnβ-(sgnn)(π2L2)EĨΔ(|n|-14)2|≤ε.In other words, for|n|≥N1, all eigenvalues {νnβ} are located in the ɛ-vicinity of the points {ν°nβ}, given by the leading asymtotical term of (2.2). Obviously, ɛ can be chosen as small as desired by manipulating the control parameters β and δ.Formula (2.3) for the δ-branch means that there exists N2>0 and ɛ̃>0 such that for all |n|≥N2 the δ-branch eigenvalues satisfy the estimate  (3.2)|νnδ-πnLĨ/G-ı2LĨ/Glnδ+GĨδ-GĨ|≤ɛ̃|n|.
This formula means that for |n|≥N2 each eigenvalue νnδ is located in a small circle of radius that tends to zero at the rate |n|-1/2 and is centered at the points {ν°nβ}; each center coincides with the leading asymptotical term from (2.3).Regarding this description, the following important question holds: from which numbersN1 and N2 can the eigenvalues be approximated by the leading asymptotical terms with acceptable accuracy? In other words, can one claim that asymptotical formulas (2.2) and (2.3) are valuable for practitioners or are they just important analytical results of the spectral analysis?As is well known, from practical applications only the first dozen of the lowest eigen-frequencies are important for engineers. The results of numerical simulations below show that the asymtotical formulas are indeed quite accurate, that is, if one places on the complex plane the numerically producedβ-branch and δ-branch eigenvalues, then the theoretically predicted branches can be seen practically immediately (see Figure 2).Figure2 means that the leading terms from asymptotics (2.2) and (2.3) can be used by practitioners as good approximations for the eigenvalues.The numerical procedure which is quite nontrivial, is briefly outlined below. The numerical approximation is based on Chebyshev polynomials. First, a finite-dimensional approximation for the main differential operatorℒβδ, (denoted by Lβδ) is described. Let ℋN be an N-dimensional subspace of polynomials of degree (N-1); obviously ℋN⊂ℋ, where ℋ is the main Hilbert space with norm (1.13). Each polynomial is uniquely determined by its values at N points. These N points are the N extremal points of the degree (N-1) Chebyshev polynomial TN-1(x) on [-1,1], defined by (3.3)TN-1(y)=cos(N-1)θ,y=cos(θ),0≤θ≤π,-1≤y≤1.
These N extrema are yk=cos(kπ/(N-1)), k=0,…,N-1. Transforming them to [-L,0] with xk+1=(L/2)(yk-1) yields  (3.4)-L=xN<xN-1<⋯<xj<⋯<x2<x1=0.
Next, the action of Lβδ is considered on ℋN, defined to be the subspace of ℋ consisting of 4-component vectors, each of whose components is a polynomial of degree N-1:(3.5)Ψ⃗(x)≡[ψ0(x)ψ1(x)ψ2(x)ψ3(x)],
where ψ0, ψ1, ψ2, and ψ3 are polynomials of degree N-1 on -L≤x≤0.The boundary conditions are imposed on the finite dimensional system by restricting the polynomialsψi(x), i=0,1,2,3, to a subspace 𝒳N⊂ℋN of the functions satisfying the boundary conditions as follows: (3.6)Ψ⃗(x)≡[ψ0(x)ψ1(x)ψ2(x)ψ3(x)]suchthat{ψ1(-L)=ψ1′(-L)=ψ3(-L)=0,ψ0′′′(0)=0,Eψ0′′(0)+βψ1′(0)=0,Gψ2′(0)+δψ3(0)=0,Aeroelasticcontrolsψ0(-L)=ψ0′(-L)=ψ2(-L)=0.
Here, derivative(s) at the boundary points means derivative(s) of the relevant polynomial component at the boundary points. In the matrix representation of Lβδ, the derivatives d2/dx2 and d4/dx4 are computed with differentiation matrices, represented by D2 and D4, respectively, [18].The finite dimensional eigenvalue problem, with polynomials represented as vectors of nodal values, then becomes the following. Find(3.7)Ψ⃗≡[ψ⃗0ψ⃗1ψ⃗2ψ⃗3]∈𝒳N
such  that ıLβδΨ⃗=λΨ⃗, that is,  (3.8)ıLβδΨ⃗=[0I00-EĨΔD4-πρuS̃Δ-S̃Δ(GD2+πρu2)-πρuĨΔ000IES̃ΔD4πρum̃Δ-m̃Δ(GD2+πρu2)πρuS̃Δ]·[ψ0ψ1ψ2ψ3]=λ[ψ0ψ1ψ2ψ3].
This system's eigenvalues are so sensitive to roundoff errors that conventional approaches to solving the eigenvalue problem are numerically intractable, even in double precision. Fortunately, there is a better way to calculate the eigenvalues of this system, which we now briefly describe. (For a more detailed explanation of these details, and the numerical difficulties, see [19].)Conventional collocation methods would impose boundary values by replacing some rows in the matrix representation ofıLβδ with rows that impose the appropriate restrictions on the vector components. To overcome the numerical difficulties, an alternative method of imposing boundary values is used—namely, representing𝒳N as the kernel of a linear operator defined as follows.The nine boundary conditions can be described by a mapping ofΛ:ℋN→ℝ9 given by  (3.9)ΛΨ⃗≡Λ[ψ0ψ1ψ2ψ3]=[ψ0(-L)ψ0′(-L)ψ2(-L)Gψ2′(0)+δψ3(0)Eψ0′′(0)+βψ1′(0)ψ0′′′(0)ψ1(-L)ψ1′(-L)ψ3(-L)],soΛΨ⃗=0means{ψ0(-L)=0ψ0′(-L)=0ψ2(-L)=0Gψ2′(0)+δψ3(0)=0Eψ0′′(0)+βψ'1(0)=0ψ0′′′(0)=0ψ1(-L)=0ψ1′(-L)=0ψ3(-L)=0}
which are exactly the nine boundary conditions. In other words, 𝒳N is the kernel of the operator Λ. It can easily be seen that  (3.10)k=dim(Ker(Λ))=dim(ℋN)-rank(Λ)=4N-9,ℝk≅𝒳N.
Let B represent a 4N×k matrix whose columns form an orthonormal basis for Ker(Λ), so  (3.11)B:ℝk→𝒳N=Ran(B),BTB=identityonℝk,BBT=identityon𝒳N.
Finally, transform Lβδ with the basis change B to an operator with the domain ℝk rather than 𝒳N, according to the diagram on Figure 3: the relationship between the eigenvalues of BTLβδB and Lβδ is as follows:
(3.12)IfΨ⃗∈𝒳N,LβδΨ⃗=λΨ⃗,thenΨ⃗=Bv⃗,forsomev⃗∈ℝk,BTLβδBv⃗=λBTBv⃗=λv⃗whichmeans{eigenvaluesofLβδrestrictedto𝒳N}⊂{eigenvaluesofBTLβδB}.
However, the reverse inclusion is not necessarily true. This means that after solving the standard eigenvalue problem BTLβδBv⃗=λv⃗, one has to use some criterion to select those λ such that LβδBv⃗=λBv⃗. The chosen criterion is to calculate both eigenvalues and eigenvectors v⃗ and choose those λ such that ∥LβδBv⃗-λBv⃗∥<ɛ. This is done in two steps.(1)
Use a QZ algorithm to calculate eigenvalue-eigenvector pairs, and calculate individual eigenvalue condition numbers (e.g., with MATLAB's package “condeig”).(2)
Use an inverse power iteration on each eigenvalue-eigenvector pair from QZ  algorithm, which improves the accuracy of the eigenvalues. The transformation of the eigenvalue problemLβδψ⃗=λψ⃗ to BTLβδBv⃗=λv⃗ produces the spectrum with the asymptotic distribution predicted by Theorem 2.2.Figure 3
Graphical representation of basis change.


## 4. Numerical Results on Flutter Modes

The second problem is related to the effect of the integral operator. As it is shown in the first round of calculations, the spectrum of the structural part of the entire aeroelastic dynamics generator splits into two branches. It is shown that the leading asymptotical terms of theβ-and δ-branches can be used as practically acceptable approximations for the eigenvalues. However, analysis of “structural eigenvalues” was restricted to the ground vibration model. It is certainly important to clarify how the integral part of the problem (i.e., the aerodynamics) affects the distribution of the eigenvalues.Obviously, one needs to select an acceptable approximation for the matrix-valued function𝔉̂(λ), or to find out which approximation would be the best for the Theodorsen function T(λ). It is well-known that the problem of finding a good approximation for this function is an important component of computational aeroelasticity. Typically, a rational function is the standard approximation for T(λ). In this work a different approximation, which appears naturally, is suggested. Namely, the asymptotic approximation for the Theodorsen function as its argument λ tends to infinity is used. Splitting the Theodorsen function yields corresponding splitting for the matrix-valued function 𝔉̂(λ).So,  the main assumption is that the influence of the aerodynamics takes place through the matrix𝔐 whose entries depend on the wind speed u. The goal of the numerical simulations below is to control the aeroelastic modes' response to increasing u, and/or changing the mass moment, S. In particular, one hopes to reach a specific speed, the flutter speed. As is known, flutter instability can be caused by one of the following two reasons. The first reason is the existence of an “unstable aeroelastic mode,” that is, the mode whose real part is positive. The second reason is that for some speed u, the two aeroelastic modes initiated by two different branches merge together creating one double mode. For speeds that are close to the critical speed (or the flutter speed), one should use the modified model problem (most likely, a nonlinear model).Below, a summary of numerical findings is presented.The graphs below show the distribution of the eigenvalues of a nonselfadjoint operatorℒβδ defined in (1.16) and (1.17). However, to find approximations for the aeroelastic modes, one should rotate each picture by 90 degrees in a counterclockwise direction.As explained earlier, numerical calculations are based on spectral methods using Chebyshev polynomials. For a typical run from 60 to 80 Chebyshev nodes were used. As many as 120 nodes were used to see whether more nodes yielded significantly better accuracy, but there actually was higher numerical error with 120 nodes than with 80 or 100 nodes. This was apparently due to the extremely large fourth-order derivatives in the matrix differential operator at nodes near the boundaries.Figure4 shows a typical display of the lowest eigenvalues of the system. The diagram shows the clear presence of the projected β-and δ-branches, when the speed u is relatively low.Figure 4
Smaller eigenvalues.Changing the parameters of the system causes the eigenvalues to shift positions. It was put forth in [1, 2] that the presence of double eigenvalues would create the flutter condition. In the search for these double eigenvalues, the only parameters of the system that have been changed were the wind speed, u, and the mass moment of the wing, S.The first significant change in eigenvalue position was observed while examining the system withu=2.1 and S on the order of 0.2. The eighth and ninth eigenvalues with positive real parts were seen to jump branches as S changed. En route in this transition, the eigenvalues are clearly near-double. In addition, the effects of fixing S and changing u were examined. It was noticed that the same eigenvalues were near-double. It should be noted that this near-double condition is also true of the mirrored eighth and ninth eigenvalues, but only the set with positive real components was examined numerically. The near-double nature of these eigenvalues is subsequently examined.The following five plots (Figures5, 6, 7) show the progression of two pairs of near-double eigenvalues with changing mass moment, S.Progression of near-double eigenvalues with changing mass moment.
(a)(b)Progression of near-double eigenvalues with changing mass moment.
(a)(b)Figure 7
Progression of near-double eigenvalues with changing mass moment.The next five plots (Figures8, 9, 10) show the progression of two pairs of near-double eigenvalues with changing wind speed, u.Progression of near-double eigenvalues with changing wind speed.
(a)(b)Progression of near-double eigenvalues with changing wind speed.
(a)(b)Figure 10
Progression of near-double eigenvalues with changing wind speed.To examine the conditions under which these eigenvalues are the closest, one can fixS and calculate the absolute distance between these near-double eigenvalues as a function of u. Figure 11 shows a typical example of such a function. Note the presence of a distinct minimum distance.Figure 11
Absolute distance between first set of close eigenvalues.Many production runs with fixedS have been performed, and it was found that both u and the absolute distance at the minimum were not unique. That is, the separation of these near-double eigenvalues is not only a function of u, but also of S. It becomes clear that in order to get a full handle on this problem, it is most useful to develop a fully three-dimensional plot. Figure 12 shows the absolute distance between these eigenvalues of note as a function of both S and u. It was interesting to see that there is a distinct low point in this figure, which occurs in the area of u=950 and S=0.19. Because eigenvalues didn’t cross over into the negative imaginary half plane, and because there is such a well defined low point in this figure, one can assume that it is in this region that flutter will be observed.Figure 12
3D Plot of absolute distance between near-double eigenvalues.It is important to note that these near-double eigenvalues were never observed to become exact double eigenvalues. Because the eigenvalues are numerical approximations, one cannot say that two eigenvalues have the same exact value, and a close examination of the path these eigenvalues take when passing close to each other reveals that they are indeed distinct. Figure13 shows characteristic paths that always have clear separation. This confirms Figure 11, which illustrates that the absolute distance between these near-double eigenvalues is never zero.Figure 13
Typical paths of close eigenvalues.


## 5. Concluding Remarks

The present paper is concerned with numerical investigation of two problems arising in the area of theoretical aeroelasticity. Namely, it has been shown that analytical formulas representing the asymptotical distribution of aeroelastic modes for a specific aircraft wing model can be used by practitioners. Specifically, the leading terms in the spectral asymptotics represent vibrational frequencies of a wing quite accurately. The second problem is to clarify the nature of such dynamical instabilities as flutter. It has been shown that the model can capture the flutter phenomenon. In particular, it was observed that two different aeroelastic modes moved towards each other to create a double mode when the speed of the aircraft approached a certain value called the “flutter speed”. This direction of research could be extended as follows. It is worthwhile to do a full study on the path of the near-double eigenvalues. The paths that they take while moving past each other may indicate the conditions in which flutter is most rapidly approached. This information may prove useful when an aircraft changes the angle of attack to circumvent flutter. A more important future research will involve the examination of other near-double eigenvalues of this system. It is shown that this near-double condition is not unique to the eigenvalues examined in this paper. However, we do not yet know whether these additional eigenvalues can be trusted numerically to be near-double. Also, it is not yet known what their minimum absolute separation is, and if it is small enough to produce flutter.Another future step in this project will be to include the full integral operator. This would be a more complete physical model of the system, and we expect that it will yield more accurate results. However, it will no doubt be more demanding computationally.



---

*Source: 879519-2010-03-07.xml*