An Interval Process Method for Non-Random Uncertain Aeroelastic Analysis

In this paper, we use the Monte Carlo simulation to study aeroelastic behavior caused by non-random uncertain free-stream velocity. For sampling, we use the interval process method. Each family of samples is defined by a correlation function and upper and lower bounds. By using this sampling method, there is no need for constructing precise probability distribution functions; therefore, this method is suitable for practical engineering applications. We studied the aeroelastic behavior of an airfoil and a high aspect-ratio wing.


Introduction and Background
Dynamical systems are sensitive to physical uncertainties, especially in unstable cases where any perturbation from a steady-state solution can be amplified. Aerodynamic forces and moments are often uncertain loads acting on aerospace structures. For example, uncertainty in free-stream airspeed causes uncertainty in aerodynamic loads, and consequently, creates uncertainty in the aeroelastic behavior of structures. Analysis of the uncertain dynamic responses is essential to the safety and reliability of structures. This type of uncertainty analysis is considered parametric uncertainty. There are other causes of uncertainty in the study of a dynamic system [1], mainly uncertainty in physical and mathematical models [2,3] (model-form uncertainty) and uncertainty in methods of solving a problem (predictive uncertainty). In this paper, we focused on studying aeroelastic parametric uncertainty.
The literature on uncertainty quantification is rich [4]. Many such methods have been applied to aeroelastic analysis [1,5]. These methods include intrusive and nonintrusive methods and are mostly illustrated on two-degree-of-freedom (2 DOF) systems [1,5]. Khodaparast and their co-authors compared the performance of fuzzy and probabilistic methods against Monte Carlo simulations (MCS). Polynomial Chaos Expansion method (PCE) was used on a composite wing [6]. The perturbation method was used on the Goland wing [7]. However, the Monte Carlo simulation (MCS) remains the gold-standard of uncertainty analysis. Hence, many authors compare their results against results from MCS. Singular value decomposition (SVD) [8] and reduced-order response surface methods [9] have been used to improve computational efficiency of MCS. MCS is a nonintrusive uncertainty analysis techniques; therefore, it can be used without any modification of the governing equations.
We use the Monte Carlo simulation with the interval process method to study the uncertainty in aeroelastic behavior of structures as a result of uncertainty in free-stream velocity. Jiang and their co-authors have initially introduced the interval process method for nonrandom vibrational analysis [19][20][21][22]. The interval process model is used to describe the uncertain time-variant loads instead of stochastic process models [19]. This method allows us to obtain the upper and lower bounds of the response instead of the probability distribution function of the response. The nonrandom vibration analysis decreases the dependency on large experimental samples because it does not require precise probability distribution functions. The bounds of dynamic responses are easy to understand conceptually and convenient to use in practical applications. In designing structures, we need to ensure that aeroelastic instabilities do not occur. Therefore it is enough to calculate the limits of aeroelastic instabilities under uncertain conditions. Jiang and their co-authors focused on structural dynamic response [19][20][21][22]. Khodaparast and their co-authors used the interval process method along with an optimization method (response surface method) to study the effect of uncertain structural parameters on the flutter speed of a 2 DOF system and the Goland wing [16].
In this paper, using the interval process method and MCS for uncertain aeroelastic analysis, we will explore the limit cycle oscillation of an airfoil and a high aspect-ratio wing (the Goland wing) under uncertain free-stream velocity.

Theory
We use the Monte Carlo simulation with non-random uncertain variables. Figure 1 shows our Monte Carlo simulation. Monte Carlo simulation has been used for decades in the statistical study of structural dynamics; however, using the interval process method [19] in aeroelastic analysis is novel.

Sampling
Each sample family is defined by the upper and lower bounds and a correlation function. One can also specify a sample family using an average value and change from the average value instead of the upper and lower limits. The following equations show the relationships between the upper (X U ) and lower (X L ) bounds, the average value ( middle point, X ave ), and the half-width value (change from the average value, δX).
One way to represent δX is δX = (∆)(X ave ), where ∆ is a fraction of X ave .
Figure 2 [21] shows a typical interval process with variable upper and lower limits in time. In general, the value of X(t i ) at a given time t i has influence on the value of X(t j ) at another time t j , and vice versa. In other words, the values of X(t) at different times are mutually correlated. A correlation function is a function that quantifies this correlation, and it can be derived from the experimental sample functions of X(t). This method has five steps [21].
Step 1: Calculate the average and half-width functions from the given upper and lower limits. Then calculate the correlation matrix, R xx , from the given correlation function, r xx (t i , t i ).
Step 2: Use the modified Cholesky decomposition [23] to create the lower triangular matrix, L x , from the correlation matrix, R xx .
Step 4: Using the lower triangular transformation of the correlation matrix, L x , transform y to x , x = L x y.
x is a new sample in the interval of [−1, 1] with correlation matrix R xx .
Step 5: Finally, use linear translation and stretching transformation on x to get the sample in [X U , X L ] with the R xx correlation matrix.
This method is particularly useful for time-dependent variables with a time scale similar to the flight time. For example, free-stream velocity can have uncertain variation during flight, while the structural stiffness changes during the lifetime of the aircraft but not during the flight time. We use samples of free-stream velocity and will study its effect on structure deformation. The free-stream velocity samples are distributed according to different probability distributions, thereby making the process generated not strictly stationary. It is only weakly stationary because of the imposed correlation and keeping the middle point and half-width functions constant.
This sampling method can be expanded to sample multi-variables, with either independent correlation functions for each sample or with a mutual correlation function. For example, researchers can use this method to study a test case with variable free-stream velocity and air density. However, in this paper, we focus on a variable free-stream velocity during flight.

Deterministic Solver for an Airfoil Cross-Section
Figure 3 [24] shows a rigid airfoil cross-section attached to two linear springs. The equations for flutter analysis of this linear aeroelastic system are derived using Lagrange's equations [24].
where m is mass, I P is the mass moment of inertia about point P, k h and k θ are the spring stiffnesses, M ac is the aerodynamic pitching moment at aerodynamic center and L is lift. We calculate L and M ac using known airfoil parameters and add unsteady effects by adding Peters' dynamic inflow and apparent mass [24]. For the results presented in this paper, we use six dynamic inflow states. The free-stream velocity is the uncertain variable, and plunge, h, and pitch angle, θ are output of this study.

Deterministic Solver for a High-Aspect Ratio Wing
A high aspect-ratio wing can be modeled with a beam theory. We have implemented the fully intrinsic beam theory [25] in our in-house research program, the Nonlinear Aeroelastic Trim Furthermore, Stability for HALE Aircraft (NATASHA) [26]. The fully intrinsic beam theory [25] is a set of geometrically exact dynamic formulation that includes initial curvature and twist, shear deformation, rotary inertia, and general anisotropy. NATASHA includes an unsteady aerodynamic solver based on strip theory. It calculates the lift, drag, and aerodynamic pitching moment for each cross-section using the given airfoil aerodynamic coefficients. The unsteady aerodynamic effects are added by adding the apparent mass and Peters' dynamic inflow [27]. NATASHA is coupled with Variational Asymptotic Beam Section (VABS) [28] for cross-sectional analysis. VABS solves a 2D finite element problem over a beam's cross-section and calculates stiffness and mass properties of the cross-section, which are NATASHA's inputs. Figure 4 shows the flowchart for this deterministic solver. The free-stream velocity is the uncertain input variable. We present the tip deformation of the wing as output variable.  Figure 4. The high-fidelity structure model, coupling VABS [28] and NATASHA [26], for aeroelastic analysis of high-aspect ratio wings.

Numerical Results
First, we show results of sampling using different correlation functions. Then we present the Monte Carlo simulations (MCS) of aeroelastic behavior of a typical airfoil section and a high-aspect-ratio wing, the Goland Wing [29].

Sampling
In this section, we show samples created using trigonometric and exponential correlation functions. Trigonometric and exponential functions are building blocks of other smooth functions, therefore we use them in this paper. For non-academic applications, a suitable correlation function is chosen based on measured data of free-stream velocity. Figure 5 shows five samples with a constant correlation function, r xx (t i , t j ) = r xx = const, when the lower and upper limits are 42 and 44, respectively. The samples are constant as expected. Figure 6 shows five samples using correlation function, r xx (t i , t j ) = r xx = Cos(ω 0 τ) for different values of ω 0 , where τ = |t j − t i | is the time span. The lower and upper limits are 42 and 44, respectively. For ω 0 = 0, r xx = 1, and we get constant lines as expected. Figure 6a-d show the larger the ω 0 , the larger the samples' frequencies. Figure 7 shows a few samples using correlation function r xx (t i , t j ) = r xx = e −α|τ| for different values of α, where τ = |t j − t i | is the time span. Again the lower and upper limits are 42 and 44, respectively. As expected, r xx = 1 when α = 0 and it creates constant line samples. As α becomes larger and larger, the samples become closer to the upper and lower bounds but never exceed the bounds. By combining an exponential function and a trigonometric func-tion, we can control both the frequency and amplitude of the samples. Figure 8a,b show five samples for r xx = Cos(ω 0 τ)e −ατ for ω 0 = 1, α = 1, and ω 0 = 5, α = 0.5, respectively.

Airfoil Cross-Section
Equation (4) are usually solved using normalized parameters. Table 1 shows the definition of these normalized parameters and the values we used. For these parameters, flutter boundary is between V = U bω θ = 2.15 and V = 2.2, where U is free-stream velocity.  Figure 9a shows five samples of V with V ave = 2.2, ∆ = 0.05, and correlation function r xx = e −5τ . Figure 9b,c show the two degrees of freedom associated with each sample, h/b and θ. All the samples show unstable behavior, despite there being instances where the free-stream velocity is below the flutter velocity. The effect of instability is dominating; therefore, all samples are unstable. Figure 10a-c show velocity (Figure 10a), and the two degrees of freedom (Figure 10b,c), respectively. In this case V ave = 2.2, ∆ = 0.1. Comparing  Figures 9a-c and 10a-c shows that increasing ∆, the half-width of velocity samples, does not change the behavior of the system, but it increases the amplitude of solution. Figure 11a   As we increase the number of velocity samples, the average of the response to all samples of velocity, V, converges to the response of the system to the average velocity, V ave . This statement is the law of large numbers in probability theory and statistical energy analysis [30]. Figure 13a-c show 500 velocity samples, V, and the response to each sample for the two degrees of freedom, h/b and θ, for V ave = 2.15, r xx = e −ατ , α = 20, and ∆ = 0.1. Figure 13d shows the average of 100, 200, 300, 400, and 500 velocity samples and the average velocity, V ave = 2.15. Figure 13e   Other correlation functions can be used for sampling. For example, Figures 15a-c and 16a-c show results using r xx = Cos(ω 0 τ), with ω 0 = 1, ∆ = 0.05, and V ave = 2.2 (Figure 15a-c) and V ave = 2.15 (Figure 16a-c), respectively. These two cases show a stable and an unstable case. Changing ω 0 will change the frequency of samples; therefore, it will not affect the stability behavior of the response.

High-Aspect Ratio Wing
We use the Goland wing [29] as a test case. Table 2 shows the properties of the Goland wing [26]. A eigenvalue analysis shows that the Goland wing flutter boundary is between 446 ft/s and 448 ft/s. Figure 17 shows the real and imaginary eigenvalues against time for this wing.  We use three average velocities, U ave , for Monte Carlo simulations of the Goland wing's aeroelastic behavior: a velocity much smaller than the instability boundary U ave = 430 ft/s; a velocity slightly larger than the instability boundary, U ave = 450 ft/s; and a velocity much larger than the instability boundary, U ave = 460 ft/s. Figures 18-20 show the results for five free-stream velocity samples with an exponential correlation function, r xx = e −ατ , with ∆ = 0.1, and α = 5 or α = 50 for U ave = 430, 450, and 460 ft/s, respectively. Figure 18b,d show the out-of-plane wing tip deformation for U ave = 430 ft/s with α = 5 and α = 50, respectively. For both cases, all samples show stable damped behavior, since U ave is smaller than the instability boundary. Figure 19b,d show the out-of-plane wing tip deformation for U ave = 450 ft/s with α = 5 and α = 50, respectively. As expected, by increasing α, the sample velocity gets closer to the upper and lower limits but remains within the bounds. Results for α = 5 (Figure 19b) show that while the initial disturbance damps out to some extent, the structure continues to oscillate. α = 50 (Figure 19d) creates unstable oscillations. This behavior is the effect of dominant velocities larger than the instability boundary. Figure 20b-d show similar results for U ave = 460 ft/s. For both α = 5 and α = 50 all samples result in unstable but bounded (limit cycle oscillation) response, because the average velocity is larger than the instability boundary.   Figure 21 shows a test case with U ave = 430 ft/s, r xx = e −ατ , where α = 5, and ∆ = 0.1. This case is stable. Figure 21a shows the average of 100 free-stream velocity samples and the average velocity. Figure 21b shows the average of all responses and the response to the average velocity. Figure 22 shows similar results for a test case with U ave = 460 ft/s, r xx = e −ατ , where α = 5, and ∆ = 0.1, which is unstable. These results show that the average of 100 responses predicts the behavior of the response to the average velocity; however, for an unstable case, more samples are required, which means modeling uncertainty is more critical when the structure is close to the stability boundary. A trigonometric correlation function controls the phase of samples. Figure 23a,c show samples of free-stream velocity for U ave = 430 f t/s with correlation function r xx = Cos(ω 0 τ) with ω 0 = 1 and ω 0 = 2π, respectively. As ω 0 increases, the frequency of U increases. Figure 23b,d show the out-of-plane tip deformation for ω 0 = 1 and ω 0 = 2π, respectively. Since the free-stream velocity samples remain in the stable region, all results are stable. Figure 24a,c show samples of free-stream velocity for U ave = 450 f t/s with correlation function r xx = Cos(ω 0 τ) with ω 0 = 1 and ω 0 = 2π, respectively. As ω 0 increases, the frequency of U increases. Figure 24b,d show the out-of-plane tip deformation for ω 0 = 1 and ω 0 = 2π, respectively. In both cases, the average velocity is in the unstable region. However, Figure 24d shows an unstable behavior while Figure 24b shows a beating behavior. Figure 25a,c show samples of free-stream velocity for U ave = 460 f t/s with correlation function r xx = Cos(ω 0 τ) with ω 0 = 1 and ω 0 = 2π, respectively. As ω 0 increases, the frequency of U increases. Figure 25b,d show the out-of-plane tip deformation for ω 0 = 1 and ω 0 = 2π, respectively. In both cases, the average velocity is in the unstable region, and the oscillations are unstable with a beating frequency.

Conclusions
We used the Monte Carlo simulation and interval process method to study the aeroelastic behavior caused by non-random uncertain free-stream velocity. For sampling, we used the interval process method. Each family of samples was defined by a correlation function and upper and lower bounds. We studied two academic correlation functions, exponential and trigonometric. Any smooth function can be modeled as a series of exponential and trigonometric functions. For non-academic applications, the suitable correlation function should be used based on experimental data of samples.
We studied the aeroelastic behavior of a 2 DOF airfoil and the Goland wing. In both cases, The general behavior depends on the average free-stream velocity. However, if velocity, within the upper and lower bounds, becomes dominantly larger than the flutter speed, an unstable behavior occurs. If the upper bound of the velocity sample is smaller than the stability boundary, then all samples cause a stable and bounded oscillation.
We also illustrated the law of large numbers: as we increase the number of free-stream velocity samples, the average of the response to all samples converges to the response of the system to the average velocity.
Author Contributions: Conceptualization and methodology: Z.S.; software: All authors; validation: All authors; formal analysis:All authors; data curation: All authors; writing-original draft preparation, Z.S.; writing-review and editing, Z.S.; visualization: All authors; supervision: Z.S.; project administration: Z.S. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.

Data Availability Statement:
Interested readers can request the presented data by contacting the first author.

Conflicts of Interest:
The authors declare no conflict of interest.