Dynamics of random pressure fields over bluff bodies: a dynamic mode decomposition perspective
Xihaier Luo Ahsan Kareem






label1]Xihaier Luo cor1
xluo1@nd.edu
label1]Ahsan Kareem
kareem@nd.edu
[cor1]Corresponding author. 156 Fitzpatrick Hall, Notre Dame, IN 46556, USA.
[label1]NatHaz Modeling Laboratory, University of Notre Dame,  Notre Dame, IN 46556, United States



    The proper orthogonal decomposition (POD), a modal decomposition technique, has been widely used for analyzing the spatial-temporal coherence of wind pressure fields, due to the complex nature of wind-structure interactions. The interpretation of characteristics and information of the extracted modes, however, is influenced by the POD algorithm, which tends to mix multiple temporal frequencies in a single mode. In this study, an operator-theoretic based decomposition approach referred to as the dynamic mode decomposition (DMD) is introduced to characterize aerodynamic pressure fields to offer an additional perspective. More specifically, the proposed approach can capture not only dominant spatial patterns but also identify each pattern with a specific frequency and a corresponding temporal evolution (growth/decay). To demonstrate the unique feature of the proposed decomposition algorithm, this study is directed towards the analysis of fluctuating pressures on bluff bodies immersed in boundary layer flows. In particular, aerodynamic pressure fields over a prism were analyzed using limited wind tunnel data and a comparison with the POD was carried out to observe commonalities and unique perspective of each decomposition. Quantification and interpretation of the decomposition results are provided. It is envisaged that DMD would supplement our understanding of the dynamics of pressure fields currently available via POD.




Aerodynamics of prism Koopman operator Dynamic mode decomposition Proper orthogonal decomposition Tall building







§ INTRODUCTION


Understanding bluff body aerodynamics is essential to estimate wind effects, for not only the performance-based design but also the control, operation, and maintenance of structures such as tall buildings <cit.>. Usually, aerodynamic characteristics such as fluctuating wind pressures on a bluff body are characterized as spatial-temporally varying random fields. Extensive research has been conducted to better extract dominant features embedded in these random pressure fields <cit.>.

Known under different names such as the Karhunen-Loève expansion, principal component analysis, etc, the proper orthogonal decomposition (POD) is one of the most widely used methods for analyzing random pressure fields. Specifically, POD analysis involves a procedure that transforms a multi-dimensional isotropic random field into a set of uncorrelated single-dimensional spatial modes, which are organized by a sequence of a continuous-valued random process <cit.>. Because the identified spatial modes are orthogonal to each other, the POD method optimally evaluates the random pressure field in a L_2 norm sense. Consequently, the first few modes ranked on the basis of proportional energies are generally sufficient to recover the dynamics. Despite its wide range of applications, POD analysis is restricted to estimating second-order statistics, and its dimension-reduced representation is limited to a linear transformation range <cit.>. More importantly, the POD method by definition is unable to directly identify single-frequency dynamic coherent structures <cit.>. Instead, it mixes a range of frequencies in a temporal behavior (See <ref>).

Note that fact the fluctuating wind pressure data, in the abstract sense, contains two elements: a set of states that describe the evolution of the pressure distribution in time, and an operator that defines the rule for that evolution. Therefore, integrating the evolution map into the decomposition process can potentially overcome the mixed multiple temporal frequencies issue encountered in POD analysis. This is the central idea of operator-theoretic modeling, especially the Koopman operator which lifts the dynamics from the state space to the space of observables that are functions of the state variables <cit.>. Following this principle, the nonlinear, finite-dimensional pressure dynamics can be mapped into a linear but infinite-dimensional Hilbert space.

However, solving such an infinite-dimensional representation can be very challenging. Different methods have been continually developed to obtain a numerical approximation of the Koopman operator. Dynamic mode decomposition (DMD), which is an equation-free method and was initially introduced by Schmid out of a need to identify spatio-temporal patterns buried in fluid flows, was found closely connected to the spectral analysis of the Koopman operator <cit.>. In particular, DMD decomposes a linearly dependent operator that is determined by a pair of time-shifted matrices to approximate the spectral properties associated with the Koopman operator. This results in a group of periodically oscillating eigenmodes. The resolution scheme is designed in an Arnoldi-algorithm like architecture, where the explicit formulation of evolution map can be saved, greatly improving the computational efficiency <cit.>.

In this study, we attempt to learn spatial-temporal pressure patterns utilizing operator-theoretic modeling method, i.e. the Koopman operator. The corresponding dynamic coherent structures will be resolved by means of an augmented dynamic mode decomposition algorithm. The augmentation centers on enriching a limited wind tunnel dataset, which is often the case in experimental research. Specifically, the manifold representing the state space is embedded into a higher dimensional Euclidean space with topological properties staying unchanged. This topologically equivalent representation is achieved by the implementation of Takens' Embedding Theorem <cit.>. In particular, delay coordinates containing informative content of the dynamical evolution are exploited. As a result, the augmented dataset is able to provide sufficient information for the operator-theoretic modeling of fluctuating wind pressures. We demonstrate the enhanced capability by applying the proposed operator-theoretic approach to learn the dynamics of random pressure fields over surfaces of a scaled model of a finite height prism. Compared to classic methods such as POD, learning results indicate the augmented DMD is able to isolate identified coherent structures with a specific frequency as well as a corresponding temporal growth/decay (See <ref>). Moreover, we connected the identified pressure patterns to the wind environmental conditions around the prism, explaining the pressure dynamics from a turbulence perspective <cit.>.

The remainder of this paper is organized as follows. <ref> presents the problem we are dealing with and offers a detailed discussion on the solution method. <ref> gives two mathematical examples to demonstrate the principles of different decomposition algorithms. <ref> describes the experimental setup of the comparative case study. Subsequently, <ref>, <ref>, and <ref> summarize the POD/DMD results including algorithm convergence, data sufficiency, and modal analysis. Finally, concluding remarks and an outlook on future work are provided in <ref>.





§ METHODOLOGY


In this work, we specially focus on analyzing the fluctuating wind pressures over a prism. Without loss of generality, the aerodynamic pressure field is modeled as a stationary and ergodic multivariate random process <cit.>:


    d/dt( x) = f ( x( t ), θ)


where x is the system state denoting the pressure distribution of a prism surface, and it is contained in the state space S ⊂ℝ^d. f: S →ℝ^d is a vector field that depends on the state x, time t, and a set of parameters θ controlling the wind-structure interaction. To apply this evolution representation to the discretized wind tunnel data, <ref> is reformulated as:


    x( t_k ) = F( x_k-1) = x( t_k-1) + ∫_t_k-1^t_k f ( τ, x( τ) ) d τ


where F: S → S is the flow map, k is the discrete time index, and τ is a sufficiently small number.

The objective is to provide physical insights concerning the aerodynamic characteristics of fluctuating wind pressures by analyzing limited measurements. However, the pressure data is highly nonlinear and exhibits multi-scale pattern in both space and time. To address this nonlinearity and complexity issue, an operator-theoretic approach is considered.



 §.§ The Koopman operator


Consider the dynamical system given in <ref>, let g: S →ℝ to be a real-valued observable that is an element of an infinite-dimensional Hilbert space. The collection of a complete set of observables constitutes a linear vector space. The Koopman operator, denoted by 𝒦, is an infinite-dimensional linear operator that acts on this vector space, given by <cit.>:


    𝒦 g ( x) = g ∘F( x) ⟹𝒦 g ( x_t) = g ( x_t+1)


where ∘ denotes the function composition operator. The key property of the Koopman operator exploited in this study is its linearity. It is associated with the linearity of the addition operation and follows the definition in <ref>:


    𝒦 [g_1 + g_2] ( x) = g_1 ∘F( x) + g_2 ∘F( x) = 𝒦 g_1 ( x) + 𝒦 g_2 ( x)


Note the Koopman operator provides a linear rule of evolution that is appealing to dynamical systems analysis. However, the emerging challenge is that the space of observables is infinite-dimensional. An schematic illustration of the Koopman operator for nonlinear dynamical systems is given in <ref>.





 §.§ Koopman spectral analysis


In essence, the Koopman operator a linear operator that is unitary on a transitive invariant set. It is therefore important to study its spectral properties as this will give valuable insights into the pressure dynamics <cit.>. Similar to the case of linear systems, one can view the Koopman operator as an infinite-dimensional matrix and perform eigen-decomposition of this matrix, where the computed eigenvalues and eigenmodes are physical objects in terms of explaining how the matrix acts on a vector space. Usually, the Koopman eigen-equation writes as:


    𝒦^tϕ_j ( x) = e^λ_jtϕ_j ( x)   where   j = 1, 2, …


where λ_j ∈ℂ is the j^th Koopman eigenvalue and ϕ_j: S →ℂ  denotes the paired Koopman eigenfunction, which is a complex-valued observable of the dynamical system stated in <ref>. For the discrete-time system, the definition is slightly different:


    𝒦^tϕ_j ( x) = λ_j ϕ_j ( x)   where   j = 1, 2, …




 §.§ Koopman mode decomposition


The eigenvalues and eigenfunctions obtained by solving <ref> contain spectral information about the underlying dynamical system. In particular, for systems with the Koopman spectrum consisting of only eigenvalues, the eigenfunctions of the Koopman operator form a complete basis set for the space of observables <cit.>. Thus, with the aim of representing nonlinear dynamics using a set of eigen-measurements that evolve linearly in time, one has to span an arbitrary observable like g in terms of these Koopman eigenfunctions:


    g ( x) = ∑_i^∞ϑ_i ϕ_i ( x)


with scalars ϑ_i representing coefficients of expansion. Dynamical systems involving multiple observables, which is more common in practice, can be similarly expanded:


    g( x) = [ g^1 ( x);        ⋮; g^m ( x) ] = ∑_i^∞ϑ_i ϕ_i ( x)
 

where g: S →ℝ^m denotes a vector-valued observable and the eigenfunctions ϕ_i provide a basis for Hilbert space. Given the linear Koopman expansion stated in <ref>, one can reformulate the dynamics (<ref>) of the measurements g as follows:


    𝒦^tg( x) = ∑_i^∞ϑ_i λ_i ϕ_i ( x)


This process is known as the Koopman mode decomposition (KMD), and was introduced by Mezic in 2005 <cit.>. It implies that when the initial state is fixed, say x = x_0, the dynamics recorded by measuring g over a vector-valued trajectory is an infinite summation (e.g. sum of sinusoids and exponentials in the context of a continuous-time system). The KMD method is particularly advantageous because the dynamics associated with each eigenfunction ϕ_i are determined by its corresponding eigenvalue.


For conservative dynamical systems, ϑ_i in <ref> is a unitary Koopman mode associated with the eigen-pair (λ_i, ϕ_i). It can be computed via the direct projection of the observable g onto the eigenfunction ϕ_i. These modes correspond to spatial coherent structures of direct measurements (e.g. the fluctuating wind pressures) and provide unique temporal patterns, where each pattern is characterized by linear growth or decay.




 §.§ Data-driven Koopman analysis


To proceed with the linearization of nonlinear dynamics using KMD, one has to obtain the Koopman spectral properties (eigenvalues, eigenfunctions, and modes). A variety of data-driven methods have been developed to address this issue. Despite specific algorithm differences, the first step of any type of data-driven Koopman analysis is always to assemble the experimental/simulation dataset that describes the evolution of observables <cit.>. In the context of wind pressure analysis, let observables be the direct spatial measurements of the pressure distribution over a surface of the prism. The goal is to identify a set of Koopman tuples (λ_i, ϕ_i, ϑ_i) that can link the state-space parameterization F to the evolution of observables 𝒦 in terms of analyzing aerodynamic characteristics (See <ref>).

<ref> graphically summarizes such data collection process. First, the wind pressure is measured at a number of discretized locations. Take the windward for instance, snapshots of the pressure distribution describe the aerodynamic evolution. Second, each snapshot is vectorized into a column vector. All of the column vectors are subsequently assembled in a matrix as the input dataset 𝒟 for decomposition. Third, similar to POD, the proposed decomposition method computes eigenmodes by solving an eigenvalue problem for a linear operator constructed from this input dataset. However, POD deals with the covariance matrix of 𝒟, whereas the proposed method partitions 𝒟 into a pair of time-shifted matrices (X and X^') and carries out the spectral analysis of the cross-correlation matrix that links the snapshots at time t_i+1 to the snapshots at time t_i. <ref> gives the details about the POD method and the following section discusses the proposed method.





 §.§ Dynamic mode decomposition


The dynamic mode decomposition, or DMD in short, is a data-driven regression method developed by Schmid in the fluid mechanics' community to identify spatio-temporal patterns from high-dimensional experimental data <cit.>. In later work, Rowley et al. demonstrated that linear expansion in DMD provides a finite dimensional approximation of the KMD, where the eigenvalues and eigenfunctions obtained by the DMD algorithm converge to the Koopman eigenvalues and eigenfunctions <cit.>. Several DMD algorithms have been proposed recently <cit.>. And the methodology adopted in this paper is the so-called exact dynamic mode decomposition (EDMD) algorithm <cit.>. In particular, we apply EDMD to an augmented dataset, where pressure data is arranged in a Hankel-matrix type following the Takens's embedding theorem. By combining delay-embedding of fluctuating wind pressure measurements with EDMD, the dominant dynamical patterns buried in the original time series dataset have been strengthened.

The computational procedures are given in <ref> and several properties should be addressed in terms of the numerical implementation. First, each element x_t_i (i=1, …, n) in the matrix 𝒟 is generated from a vector space S ⊂ℝ^n_s, where n_s denotes the spatial dimension, i.e. the number of pressure taps used in a wind tunnel test. Second, we form the Hankel matrix by stacking the elements of time-delayed coordinates, where the index i indicates the embedding number. Third, the goal is to find the best-fit linear operator that connects the time-shifted matrices X^' = 𝒜 X, where 𝒜 asymptotically approach the Koopman operator 𝒦. Fourth, the singular value decomposition (SVD) is a robust algorithm that provides converged POD modes on observables of ergodic dynamical systems. In <ref>, U∈ℝ^(n_s × i) × (n_s × i), V∈ℝ^(j-1) × (j-1), Σ∈ℝ^(n_s × i) × (j-1), and * denotes the conjugate transpose. U, whose columns are numerical realizations of POD modes, follows U^*U=I, and V, whose columns are temporal behaviors, holds V^*V=I (See <ref> for details). Fifth, an appropriate choice of the truncation number r can be decided by means of hard-threshold computing. Sixth, the lower rank representation 𝒜̃∈ℝ^r × r integrates the dynamical information by exploiting X^'. Seventh, the eigenvalue-eigenvector pairs for 𝒜̃ can be obtained by solving the eigen-equation stated in <ref>, where columns of Φ are eigenvectors. Λ is a diagonal matrix with non-zero entries representing eigenvalues. The growth/decay rates and corresponding frequencies can be further determined from the eigenvalue sequence λ_1, λ_2, … via σ_i = log( ℛe ( λ_i ) )/2 πΔ t and ω_i = log( ℐm ( λ_i ) )/Δ t, respectively. Eighth, it should be noted that the DMD modes computed by <ref> are commonly referred to as the exact DMD modes. An alternative approach is to multiply left singular vectors by eigenvectors ϑ^' = U_r Φ, where each column of ϑ^' is so-called the projected DMD mode. The former definition is closely related to X while the emphasis of the second definition is placed on X^'. By orthogonally projecting ϑ^' to the range of X, the projected DMD modes can be exactly represented via ϑ.





§ MATHEMATICAL EXAMPLES


To clearly illustrate the unique feature of using augmented dynamic mode decomposition (<ref>) to distill dynamic patterns directly from available datasets, two synthetic nonlinear systems are considered herein. Moreover, the DMD results are compared against the results obtained using the POD and independent component analysis (ICA) method <cit.>.



 §.§ Example 1: two mixed spatio-temporal signals


In the first example <cit.>, the nonlinear synthesis dataset of interest can be expressed by summing the following two functions:


    [          f_1(x, t)=sech(x+3) exp (i 7 t); f_2(x, t)=2 sech(x) tanh (x) exp (i 3 t) ]


The evolving system is simulated by f(x, t) = f_1(x, t) + f_2(x, t), where the two frequencies involved in the simulation are ω_1 = 7 and ω_2 = 3. In <ref>, f_1(x, t) and f_2(x, t) are simulated in the interval of [0, 4π]. For the purpose of comparison, the same amount of data have been fed into the DMD, POD and ICA algorithm. It should be noted that the imaginary values have been removed in the case of ICA as the whitening process included in ICA decomposition requires <cit.>. The decomposed results are detailed in <ref>. For DMD, the truncation number specified in <ref> is initialized to 2. The results reveal that the two dynamical patterns obtained via DMD coincide with the true evolutions and the DMD modes are almost identical to the true solution. On the other hand, POD provides a better approximation of the spatial modes than ICA while ICA excels POD in terms of learning the dynamic evolution of identified patterns.







 §.§ Example 2: mixed 2D oscillations


In this example, a synthetic dataset with two hyperbolic oscillations is created for testing the performance of different decomposition methods. The spatial modes take the following formula:


    [ u_1 = exp(-((x+10)^2+(y+10)^2)/100); u_2 = exp(-((x-10)^2+(y-10)^2)/100) ]




Specifically, the spatial domain Ω⊂ [-30, 30] × [-30, 30] with 70 elements per axis, resulting in a resolution of 70 × 70 = 4900 per snapshot. The first mode u_1 oscillates at 7 Hz and the second mode u_2 oscillates at 3 Hz. Note both modes are stationarily oscillate for the entire during of the simulation. The data is recorded with δ t = 0.01 and a total of 1000 snapshots have been collected for analysis. <ref> shows the states of the first 10 equispaced time instances. Apparently, two modes have overlapped spatial destitutions and they are evolving at different frequencies. To accelerate the decomposition process, the method of snapshot is considered for the implementation of DMD and POD <cit.>, that is, the input matrix is reshaped into the size of 𝒟∈ℝ^4900 × 1000. Again, DMD accurately captures the distribution of true modes compared to the POD and ICA. Beyond that, the evolution provided by the DMD algorithm matches well with the true dynamics. In terms of POD, it focuses on finding the eigenmodes that are energetically optimal regarding the reconstruction of the input data. Hence, the POD modes are very similar to the spatial average results. ICA modes are more isolated and results become ambiguous around the overlapped area.







§ CASE STUDY





 §.§ Wind tunnel data


In the following case study, time-series data of aerodynamic pressure fields over a finite prism (0.5 m × 0.1 m × 0.1 m) is analyzed. The experimental dataset was collected from an open-circuit wind tunnel of dimensions 14 m × 1.2 m × 1.0 m <cit.>. Specifically, the approach flow mean speed profile along the height was defined by a power law exponent of 1/4, and tests were conducted at a wind speed of 11.1438   m/s. To measure the pressure, 500 smart sensors have been installed on the scaled model. Specifically, there were 125 sensors uniformly distributed on each surface (i.e. a  5 × 25 distribution pattern in terms of the horizontal and vertical axis, respectively). All taps were synchronously recorded at a sampling frequency of 1000   Hz for a sample period of 32.768   s. The collected pressure data was transformed into a dimensionless number, that is, the pressure coefficient , and matrices X, X^' were accordingly assembled using <ref>. For more information about the experimental data, please refer to http://wind.arch.t-kougei.ac.jp/system/eng/contents/code/tpuTPU database.



 §.§ Case study outline


To demonstrate the unique features of using the Koopman operator to describe a random pressure field, we carried out a comparative study. Three types of analyses have been performed.

Convergence Study The objective of the convergence study is to determine how many snapshots are required to ensure the decomposition results are not affected by changing the number of snapshots <cit.>. In the case study, the input dataset 𝒟 consists of 32768 pressure states. Thus, a reasonable snapshot number N_snap is crucial since a small N_snap (e.g. N_snap = 30) cannot achieve converged results while a large N_snap (e.g. N_snap = 30000) would call for an additional computation. Two evaluation metrics are defined to measure the convergence behavior (See <ref>).

Sufficiency Analysis The sufficiency analysis of POD centers on constructing a low-rank representation of the fluctuating wind pressures. For the proposed augmented EDMD (See <ref>), one has to select an embedding number in such a way that the input dataset contains sufficient information for the decomposition analysis <cit.>.

Modal Analysis The goal of the modal analysis in the context of analyzing aerodynamic pressure fields is to determine the relationship between the aerodynamic loading and the identified eigenmodes. We explained DMD results from both spatial and temporal perspectives, and connected DMD modes with the turbulent flow environment around the prism.

In this study, pressure data of the wind direction α = 0^∘ was used. The complete dataset includes wind-induced pressure coefficients on four surface of the prism. For the sake of brevity, decomposition results of the windward surface have been summarized in this paper. We report the complete results in the supplementary material, which is available at <https://github.com/Xihaier/DMD-POD-aerodynamic-pressure-analysis>.



§ PART 1: CONVERGENCE STUDY RESULTS




 §.§ POD Results


<ref> shows the convergence behavior of POD analysis. The first index based on the computed eigenmodes points out the first three fundamental POD modes ϕ_1^POD, ϕ_2^POD, ϕ_3^POD exhibit consistent convergence properties. The evaluation metric κ decreases dramatically within the time span [ 0, 1000 Δ t ]. On the contrary, higher POD modes ϕ_j^POD( j = 4, …, 7 ) not only converge at a slightly slower pace but also show more random fluctuations during the convergence process. This reveals the fact that incoherent noise tend to dominate and even likely appear as higher order POD modes <cit.>. Thus, the POD method can be employed as a denoising technique for preprocessing received signals that are strongly contaminated with noise <cit.>. 



On the other hand, the second index that is defined as a function of the computed eigenvalues confirms the presence of such noise effects. The first few normalized POD eigenvalues converge at a faster speed with less local perturbation than those held by higher POD modes. Meanwhile, it is known that the normalized eigenvalues λ_j^† represent the mean energy distribution of the aerodynamic pressure data. Shown in <ref>, some POD modes have similar energy distribution as their corresponding converged eigenvalues are almost identical. From an algorithmic perspective, it is due to the data projection procedure involved in the POD analysis. This fact implies the existence of second-order statistics, i.e., correlations between projected modes <cit.>. 



 §.§ DMD Results


<ref> shows the computed κ_j^DMD and λ_j^†. According to the results, neither DMD modes nor the normalized DMD eigenvalues converge during the analysis. This phenomenon is caused by the DMD algorithm, where dynamics are integrated via a pair of time-shifted matrices X and X^' <cit.>. When the size of X is relatively small, DMD tends to produce a vastly expanded frequency band, where the amplitude spectra covers more isolated frequencies. As the size of X increases, the observed oscillation frequencies of identified DMD modes cluster themselves in a narrowband. Thus, the observed frequencies may overlap with each other, thereby obscuring the convergent process.



Convergence of the frequency band To demonstrate the aforementioned shrinkage property, we computed the DMD energy spectrum using a different number of snapshots. In <ref>, it can be seen that the frequency band shrinks rapidly when N_snap: 300 → 500 → 700. DMD modes are sparsely scattered throughout the interval [0, 500] when N_snap = 300 while the majority of resolved DMD modes are narrowed within a smaller interval [0, 100] when N_snap = 700. Specifically, some eigenmodes identified by the DMD method may highly resemble each other when N_snap is large enough. These modes share a global pattern in terms of explaining the aerodynamic pressure, but more importantly, they provide different interpretations at a local scale, which are closely related to the small scale turbulence effects. Therefore, the DMD method have achieved the convergence of the frequency band rather than the converged DMD eigenvalues and eigenmodes. 



High/low frequency dynamics To investigate the dynamics in the converged narrowband, we examined the DMD energy spectrum with N_snap = 800, 900, 1000. For the illustration purpose, frequencies have been normalized by the geometric parameter B and hourly average wind speed U. The computed results confirm the convergence of the DMD frequency band. In <ref>, the low-frequency region (blue dots) exhibits more distribution variation while the high-frequency region (yellow dots) is more unified and stable in against the N_snap change. This observation agrees with the wind environment around the prism, where large eddies of lower frequencies are more energetically unstable <cit.>. In fact, each dot has a corresponding spatial pattern, which is affected by a multitude of eddies of different sizes and strengths. Due to the richness of turbulence dynamics, the decomposition results vary from case to case. This again explains the dependence of the DMD method on the N_snap value. 





§ PART 2: SUFFICIENCY ANALYSIS RESULTS




 §.§ Sufficiency of the POD eigenvalues


According to the convergence study results, the use of 1000 successive snapshots is able to produce well-converged POD modes (See <ref>). Thus, we performed the POD analysis using 1000 snapshots. <ref> shows the sufficiency analysis results. The ratio of the cumulative sum of the POD eigenvalues is presented, where each eigenvalue is normalized by the total energy summation ∑_i=1^125λ_i. When the original dataset 𝒟 is used, the first eigenvalue ranked by the magnitude captures more than 90 percent of the energy buried in 𝒟. In many other studies, it was found that the first resolved POD mode is highly closed to the mean pressure distribution at different locations <cit.>. To investigate the aerodynamic fluctuation effects, we further extracted the mean pressure and applied POD to the remaining data 𝒟 - 𝒟̅. In the windward scenario, 7 POD modes are sufficient to reconstruct a fluctuation dataset that contains 90 percent of the energy of 𝒟. Following a similar exponential decay pattern, a need for more POD modes is required in the other three scenarios (See the supplementary material for more details). The difference reveals the fact that the intrinsic complexity of the aerodynamic pressure field on the prism surface is strongly related to the incoming turbulence, which comes from a wide range of scales. 





 §.§ Sufficiency of the delay-embedding and DMD eigenvalues


To check whether the input dataset contains sufficient information for the DMD analysis,
we compared the decomposition results of different embedding number. <ref> shows the distribution of the computed DMD eigenvalues λ_1, …, λ_125. In practice, eigenvalues are found to be unstable if their complex modulus is greater than 1 (i.e., are located outside the unit circle) and eigenvalues are deemed as stable or neutrally stable if √(ℛe ( λ_i )^2 +  ℐm ( λ_i )^2)⩽ 1 (i.e., are located inside or on the unit circle) <cit.>. According to the results, the DMD eigenvalues are scattered inside the circle when the delaying coordinates are not considered i=1. Such distribution indicates the decomposition algorithm needs more information in terms of separating different types of dynamics because ideally the eigenvalues are expected to be located at the circular perimeter. To address this issue, the input dataset is augmented by the Takens embedding theorem <cit.>. When the embedding number i: 1 → 2 → 4 → 8, it can be seen that the eigenvalues λ_1, …, λ_125 are gradually pushed to the circle boundary.



Similar to the sufficiency analysis of the POD eigenvalues, we ranked the complex-valued sequence λ_1^DMD, …, λ_125^DMD by the absolute magnitude of the real part in descending order and plotted the cumulative sum of DMD eigenvalues. The results are summarized in <ref> and several properties should be noted. First, the red line representing the complex modulus summation pattern gradually converges to diagonal reference line when the embedding number increases. The reference line is a simple linear function with zero intercepts and the slope is determined by the truncation number i, i.e. y = 1/125 x. Unlike POD results, rapid growth of the cumulative eigenvalues is undetected in ordered DMD results ℛe ( λ_1^DMD) > ℛe ( λ_2^DMD) > … (See the blue lines in <ref>). This is due to the fact that the nonlinear aerodynamics is spanned by an infinite Koopman modes. It is, therefore, energies of the computed Koopman modes asymptotically converge to a uniform distribution and the cumulative summation becomes linear <cit.>. Second, both real and imaginary lines appear to be smoother when i becomes large, affirming the enrichment effects brought by the delay coordinates. Third, the green line that covers the imaginary information of DMD eigenvalues is a downward convex function. Subsequently, the gradients of this function correspond to the phase information. It can be seen that the identified DMD modes with larger energy content tend to evolve at a slower pace than the low energy modes. This is similar to the turbulence where large eddies break into small eddies, which grow with a relatively higher frequency <cit.>.





§ PART 3: MODAL ANALYSIS RESULTS




 §.§ POD modes


Connections with statistical moments POD analysis by definition relies on the second-order statistical information <cit.>. To show the connection between POD modes and statistical moments, we summarized the mean, variance, skewness, and kurtosis of the fluctuating wind pressures in <ref>. (A). It can be seen that the the first POD mode is very similar to the first-order statistical moment and the second-order statistical moment may be generally accessed via the combination of a set of POD modes. However, variability properties included in the skewness and kurtosis cannot be explained using POD modes. From a statistical learning perspective, the POD algorithm can capture the location and variability of the data but have difficulty in describing the symmetry and flatness properties.

Spatial interpretation of POD modes Singular value decomposition (SVD) is adopted to find an orthogonal set of eigenfunctions and the computed POD modes are shown in <ref>. (B). In common practice, POD modes are assumed to be related to the energy distribution of random pressure field <cit.>. For instance, ϕ_1, ϕ_2, and ϕ_3 represent the pressure distributions that result in the along-wind overturning moment while ϕ_4 exhibits a pressure distribution pattern that is related to torsional loads. Moreover, higher-order POD modes are more associated with small scale features since eigenmodes are arranged in decreasing order of eigenvalues. For example, peak values of ϕ_4 are found around the prism edges. Similar results have been observed in the leeward and two sidewards where min-max values are around corners, suggesting the locates where large pressure fluctuations take place (See the supplementary material for more details). 

Temporal interpretation of POD modes The time coefficients associated with the computed POD modes are presented in <ref>. (C). At first glance, the evolution process looks like stochastic in nature rather than a well-ordered process. For instance, a_2 ( t ), a_3 ( t ), and a_4 ( t ) are changing abruptly in time, precluding identification of qualitative description of aerodynamic characteristics. This is caused by the zero-time-lag covariance matrix used in the POD analysis. Because such covariance matrix does not contain any evolution information, the decomposition results in assigning an eigenmode with multiple frequencies in such a way that the L_2 norm error in the POD formulation is minimized <cit.>. Consequently, the computation process precludes the prediction of the next state in time.





 §.§ DMD modes


DMD setup To compute DMD modes, we defined parameters involved in the augmented EDMD algorithm utilizing the convergence and sufficiency analysis results presented earlier. First, 730 snapshots have been used, that is, the input dataset 𝒟∈ℝ^125 × 730. Next, the input dataset was embedded into a higher dimension to ensure that most eigenvalues are driven to a neutrally stable state. In particular, the embedding number i = 30 and the embedded dataset 𝒟∈ℝ^3750 × 701. According to <ref>, the augmented decomposition algorithm would divide 𝒟 into a pair of time-shifted matrices X, X^'∈ℝ^3750 × 700. We fixed the truncation number r = 125. As a result, a low-rank representation consisting of 125 DMD modes has been extracted.

Clustering results Because the decomposition process generates 125 eigenmodes, an unsupervised machine learning method, i.e., clustering is adopted to group DMD modes of similar properties. Specifically, the k-means clustering scheme has been implemented on account of the fact that it can control the number of clusters <cit.>. We initialized the cluster number N_cluster = 3, normalized the 125 DMD eigenmodes, and specified the distance metric involved in the cluster analysis to the correlation type. As a consequence, there are three types of DMD eigenmodes categorized by the scale of spatial patterns, i.e. macro-scale, meso- scale, and micro-scale. On the other hand, we generally ranked the DMD eigenmodes by the real part of the eigenvalue in descending order. And 125 ranked DMD eigenmodes were then classified into three groups, where group 1 is {ϕ_1, …, ϕ_42}, group 2 is {ϕ_43, …, ϕ_84}, and group 3 is {ϕ_85, …, ϕ_125}. We approximated the distribution of DMD eigenvalues and frequencies associated with each cluster/group utilizing a normal kernel smoothing function. In <ref>, it can be observed that the clustering and physical grouping results are highly similar. In particular, macro-scale DMD modes tend to have larger eigenvalues while the normalized eigenvalues of micro-scale DMD modes spread over the entire interval. Such finding indicates macro-scale DMD modes contain more energy. In light of the oscillation frequency, macro-scale DMD modes have a relatively lower frequency range compared to the other two types. This observation agrees with the experimental finding on the turbulence of different length scales, where frequencies of large eddies are much smaller than small eddies <cit.>.



Energy cascade The clustering results indicate there exists a link between identified DMD modes and turbulent eddies formed in the surrounding wind environment. In this study, interest is in inferring the physical connections between clustered DMD modes and turbulent dynamics generated by the incoming winds. For example, it is known that the aerodynamic fluctuating loads acting on the windward face of the prism model are dominated by the approach flow turbulence. And aerodynamic loads acting on the two side and leeward faces are closely related to the wake dynamics including vortex shedding, which in turn is influenced by the turbulence in the approach flow <cit.>. These experimental findings show the similarity of turbulence properties with the aerodynamic loads represented by combining the DMD modes and their time coefficients. To understand how turbulence related to the aerodynamics pressure field around prisms, we focus on explaining the energy transfer by investigating the dynamical evolution of each DMD mode. In particular, the dynamical evolution process is similar to the energy transferring mechanism in turbulence, which includes energy cascade, neutrally steady, and inverse cascade <cit.>. Furthermore, it should be noted that the spatial patterns in the pressure field are similar to the distribution of length scales (eddy sizes) in the incoming turbulent flow. <ref> schematically highlights these observations.



Macro-scale DMD modes: linkage to POD modes <ref> presents the first four DMD modes in windward scenario. A general view of the resulting DMD modes shows that coherent structures at macro-scale provide a global sketch of the aerodynamic loads. Specifically, DMD modes in the windward scenario are in good agreement with the force pattern reasoning the overturning moment. Locations of local minima and maxima serve as scaled features of lateral forces. The distinct separation of minimum and maximum values indicate the aerodynamic moment on the windward surface computed using these horizontal forces directly contribute to the overall overturning moment. Very similar results have been found in the POD modes <cit.>.



Temporal evolution of macro-scale DMD modes: linkage to mean flow field Because macro-scale DMD modes are more energetic than other types of modes, they are usually assumed to be connected with, to some extent, a collection of flow mechanisms such as the mean flow field around the object <cit.>. However, one should be aware that a particular DMD mode may be simultaneously related to a few flow mechanisms. In this paper, we are mainly concerned with the mean flow field when analyzing the macro-scale DMD modes. Overall, there are two types of evolution pattern (See <ref> and the supplementary material for more information). First, the absolute amplitude of the first DMD mode appears to grow in all scenarios. This observation conforms to the physical principles of the energy injection concerning stationary boundary layer flows. Specifically, the continually incoming wind constitutes the mean energy injection that balances the mean energy dissipation <cit.>. Second, other DMD modes mostly exhibit decay trend in time. This is a strong resemblance to the turbulent energy cascade, where turbulent eddies of a wide range of scales pass energy successively from larger eddies to the smaller ones <cit.>. When wind flows over a prism, macro-scale spatial patterns reflect the large-scale fluid motions involved in the approaching wind. Observing at time t, the energy associated with these fluid motions is gradually dissipated by the viscosity, resulting in a decaying temporal behavior at t+Δ. And macro-scale coherent structures decompose to a smaller scale (meso or micro). Hence, analyzing the time coefficients of macro-scale DMD modes offer evolution dynamics regarding the break down of fluid flow.



Micro-scale DMD modes: linkage to small eddies Aerodynamic forces are usually separated into mean and fluctuating components <cit.>. And the fluctuating component related to the magnitude of pressure perturbations is directly affected by the nearby small eddies. According to the Kolmogorov length scale theory <cit.>, small eddies tend to have high frequencies, causing turbulence to be isotropic and homogeneous at a local scale. In a generic sense, small eddies represent instabilities of local interactions between different fluid motions. In <ref>, four micro-scale DMD modes from the windward scenario were selected to illustrate such pressure perturbation. Two properties of the resulting DMD modes should be noted. First, pressure perturbations are randomly distributed on each surface. Second, the active areas are relatively small and sparsely scattered. These two properties conform to the characteristics of small eddies (randomness and size) <cit.>.



Temporal evolution of micro-scale DMD modes: linkage to energy dissipation From a physical standpoint, micro-scale DMD modes do not necessarily represent any specific flow patterns regarding the surrounding flow motions. In <cit.>, it was pointed that the least energetic modes may contain information on the interaction between different mechanisms instead of actually representing one of them. The conclusion was drawn based on the POD analysis results, where temporal behavior of identified eigenmodes are hard to interpret. In <ref>, we summarize the time coefficients of selected micro-scale DMD modes. It can be clearly seen that all the selected micro-scale DMD modes exhibit decay form. In particular, the amplitude of some modes decays exponentially, for instance, mode 111 of the windward surface. This reflects the turbulent energy cascade process where large eddies generated by the mean flow cascades down to small eddies, and the kinetic energy carried by small eddies ultimately dissipated to heat in the pressure field <cit.>.



Stripe pattern: linkage to instantaneous pressure increment In light of micro-scale DMD modes, there is a specific pattern that is detected on all surfaces of the prism. We call it the stripe pattern because the computed spatial distribution indicates the surface is divided into a few bands. Pressure values of two adjacent bands come from different ranges, thus resembling the stripe pattern observed in nature. To understand these stripe patterns, we examined the evolution of the pressure distribution, finding such micro-scale pattern is closely related to the pressure increment between instantaneous spatial pressure distribution. <ref>, the first row shows the pressure distribution at two consecutive time instances t_1 and t_2 as well as the instantaneous pressure increment between p(t_1) and p(t_2). The second row displays the identified stripe pattern.



Temporal evolution of the stripe pattern: linkage to the beating phenomenon Eddies of different sizes have their corresponding locations in a frequency spectrum. In <ref>, sizes of the stripe patterns are structurally similar, indicating the size of the surrounding eddies are very close. Consequently, the natural frequencies of these turbulent fluid motions interacting with the prism are close. These factors naturally fulfill the precondition of the beating phenomenon, which in principle is a constructive interference pattern between two signal having almost identical frequencies <cit.>. The decomposition results describing the dynamical evolution of each stripe pattern mode are given in <ref>. (A). The identified DMD modes oscillate up and down and the amplitudes exhibit a decaying pattern in a general sense. This is due to the fact that small eddies are gradually damped out by viscosity. Moreover, the beat phenomenon that are highlighted by the red dots only happens in the high-frequency region (See <ref>. (B)). An intuitive explanation is based on the definition of the beating phenomenon, which not only requires the existence of two fluid motions sharing very close frequencies but also expects these two eddies are closely located in space. However, lower frequencies are associated with larger fluid motions that are energetically unstable. This precludes the existence of the beating phenomenon in the low-frequency region.



Description in discrete Fourier domain Chen, K.K et al. pointed the DMD method provides identical results to the temporal discrete Fourier transform (DFT) if the input data is linearly independent and zero centered <cit.>. Despite the fact that our aerodynamic pressure data is highly nonlinear and the mean subtraction is not carried out, DMD spectrum was found to provide useful information regarding the pressure dynamics from the frequency perspective <cit.>. We studied the mode amplitude as a function of frequency, which can be computed from the imaginary part of the DMD eigenvalues. Such spectrum captures the magnitude and phase of pressure dynamics through a linear combination of eigenvectors, where every eigenvector grows or decays at a specified frequency according to its associated eigenvalue. For comparison, we computed the amplitude of oscillations at a wide range of frequencies by means of the Fourier transform, where the temporal DFT model is essentially a superposition of a series of harmonic oscillators. <ref> shows that the DMD spectrum closely resembles the FFT spectrum and the proposed augmented EDMD algorithm can capture the leading amplitudes of the power spectrum computed by the FFT. The difference between DMD and FFT spectrum comes from the algorithm. In particular, FFT spectrum is computed independently using the collected time series data recorded via each pressure tap while every single point in the DMD spectrum denotes a certain coherent structure covering 125 pressure taps over each surface.



Review of the temporal evolution As wind passes by a building, it results in the formation of Karman vortex shedding. As these vortices are advected by the mean flow, they result in interactive periodic pressure signatures on the building side faces <cit.>. Few attempts have been made to connect identified eigenmodes to the convection of the vortices because of the spatiotemporal nature of these recurrent patterns. In <cit.>, the effects of the advection of vortices on the pressure patterns have been identified by repeated operations of digital filter and Hilbert transform of the convolutive data. In this context, the essence of the DMD lies at the use of a paired time-shifted matrices in the computation of eigenmodes, which naturally integrates the evolution dynamics into the decomposition process. In addition, the proposed augmented DMD algorithm (<ref>) exploits the delaying coordinates, serving as an alternative to the convolution operation adopted in <cit.>.



<ref> shows the temporal evolution of the pressure field at different scales. The time between snapshots equals the sampling period of the measured data 10 Δ t, 5 Δ t, and Δ t for the macro-, meso-, and micro-scale modes respectively. A total length of 7 sequences are presented and a better animated illustration of the dynamical pressure patterns is available at <https://xihaier.github.io/>. The macro-scale sequences reported in <ref> clearly shows the development process of the horseshoe vortex and the meso-scale captures the recurrent Karman vortex shedding. Meanwhile, the results confirm the spatial-temporal evolution of selected macro-scale DMD modes can capture the distribution of overall pressure dynamics. Developed at lower frequencies, macro-scale DMD modes serve as the base of aerodynamic pressure fields from an energy perspective. In the case of meso/micro modes, we can observe more local fluctuations with relatively higher frequencies and lower spatial correlation.



§ CONCLUDING REMARKS


This study investigates the dynamics of random pressure fields over bluff bodies with an application to pressure field over a prism immersed in a turbulent boundary layer. An operator-theoretic approach is introduced for the aerodynamic characterization using limited wind tunnel test data. First, the highly nonlinear pressure dynamics is linearized in the entire basin by an infinite-dimensional linear operator, i.e. the Koopman operator. Subsequently, spectral analysis of the Koopman operator provides a spatio-temporal characterization of the pressure data. To compute the Koopman eigen-elements, an augmented dynamic mode decomposition method is proposed. The augmentation is achieved by the use of the Takens’s embedding theorem, where time delay coordinates are considered. As a result, the computed DMD representation is spatio-temporally orthogonal where each identified coherent structure is assigned to a specific frequency and a corresponding temporal growth/decay.

To demonstrate the unique feature of the proposed decomposition method, a comparative study between classical POD and proposed augmented DMD is carried out. Unlike POD, which tends to mix the temporal frequencies in the identified modes, DMD provides single-frequency eigenmodes. Furthermore, it is found that the computed DMD modes are closely connected to the wind field approaching the prism. Specifically, similar to the turbulence cascade, macro-scale DMD modes break up and lead to the energy cascade, and micro-scale DMD modes are related to the instantaneous pressure increments and small eddy dynamics.

Future directions regarding applying the Koopman operator to learn similar nonlinear dynamical systems include: (1) Extension of the current work to the extreme learning condition where available data is merely a sequence of a scalar, and (2) Moving from the current Laplacian determinism to learning the intrinsic orderliness governing chaotic system.



§ ACKNOWLEDGMENTS

This work was supported in part by the National Science Foundation (NSF) under Grant No. 1562244. This support is gratefully acknowledged. The authors also acknowledge gratefully the use of Tokyo Polytechnic University aerodynamic database. The codes and data used in this work will be made available at <https://xihaier.github.io/> upon publication of this manuscript.




elsarticle-num




§ APPENDIX. PROPER ORTHOGONAL DECOMPOSITION


This appendix provides a quick summary of the proper orthogonal decomposition (POD) method. Let 𝒟 be the input dataset. Under finite energy assumption, 𝒟 belongs to the Hilbert space ℋ = L^2 ( Ω_T×Ω_x) of Lebesque measurable. Because the tensor space is dense in L^2 ( Ω_T×Ω_x), the input dataset takes a tensor product form:


    𝒟 = 𝒜( Ω_T) ⊗Φ( Ω_x)


where 𝒜( Ω_T) ≐ L^2 ( Ω_T), Φ( Ω_x) ≐ L^2 ( Ω_x), and Ω_x represents the spatial domain. As a result, the spatio-temporal representation of 𝒟 is explicitly converted to a bilinear composition form, with the first vector space 𝒜 producing the time-dependent coefficients and the second vector space Φ giving spatial modes. Mathematically speaking, the preceding equation can be equivalently expressed as a linear superposition using a complete set of basis functions according to the Schmidt decomposition theorem:


    𝒟 = ∑_j=1^N a_j ( t ) ϕ_j ( x )


Unfortunately, N is a considerably large number or infinity in most cases. The essence of POD is to seek an orthogonal projector υ( x, t ) that minimizes the approximation error:


    𝒥( υ) = || 𝒟 - υ( x, t ) ||_n_F


with || · ||_n_F denoting the Frobenius norm. In the context of POD, n_F = 2, which guarantees the kinetic energy interpretation of its optimized modes ϕ_j, j=1, 2, …, N, and the projector υ( x, t ) can be directly constructed through utilizing a finite set of orthogonal basis functions (M ≪ N) following the principles stated in <ref>:


    𝒟≈υ( x, t ) = ∑_j=1^M a_j ( t ) ϕ_j ( x )


Then, the optimal basis functions can be obtained by maximzing the average projection of 𝒟 onto Φ:


    < 𝒟, Φ( x) >_Ω_x/< Φ( x), Φ( x) >_Ω_x


By applying the calculus of variations, the optimization problem stated in <ref> is transformed to a Fredholm eigenvalue setting, where the optimal POD basis is exactly the eigenfunctions in the following integral equation:


    ∫_Ω_x R ( x, x^') ϕ_j ( x^') d x^' = λ_j ϕ_j ( x)


where R ( x, x^') is the symmetric spatial correlation tensor that is computed by taking the inner product of the pressure data in Ω_x:


    R ( x, x^') = < υ( x, t ) ⊗υ( x^', t ) >_Ω_x


Two important properties associated with <ref> should be observed. First, Hilbert-Schmidt theory assures the eigenvalue problem has an infinitely countable pair of eigenfunctions and eigenvalues {ϕ_j, λ_j }, out of which there exists a sequence of positive eigenvalues λ_1 ⩾λ_2 ⩾…⩾ 0. POD modes are ranked according to the eigenvalues of R ( x, x^') ϕ_j ( x^') by a descending order. Second, eigenfunctions are orthogonal ∫_Ω_xϕ_i ( x) ϕ_j ( x) d x = δ_ij, that is, the POD algorithm seeks the least number of modes for describing 𝒟.



§ APPENDIX. EVALUATION METRICS FOR CONVERGENCE STUDY


We checked the convergence of the decomposition results using both eigenmodes and eigenvalues. The first evaluation metric is a function of identified eigenmodes. It is used to analyze the convergence behavior of the spatial pattern. Consider two datasets 𝒟^1 = {x( t_1 ), …, x( t_m) } and 𝒟^2 = {x( t_1 ), …, x( t_n) }, which contain pressure signals of the same prism surface but with a different sampling interval [t_1, t_m] and [t_1, t_n]. Two sets of eigenmodes Φ^m and Φ^n can be accordingly computed by means of either POD or DMD method. To measure their relative difference, a L_2 norm based metric is defined:  


    κ_j = < ϕ^m_j - ϕ^n_j, ϕ^m_j - ϕ^n_j >_Ω_x


In this study, n = 30000, serving as the reference point in investigating convergence properties, and m is sampled from the interval [Δ t, n Δ t].

Meanwhile, the second evaluation metric is a function defined based on the normalized eigenvalues:


    λ_j^† = λ_j/∑_i=1^nλ_i ( m )


where λ_j^† primarily examines the temporal convergence properties. It should be further noted that the DMD eigenvalues λ_j^DMD and eigenmodes ϕ_j^DMD are complex-valued. For this reason, complex modulus that measures the magnitude of a tensor from the coordination center to its location in the complex plane is utilized:


    λ_j^DMD    = √(( ℛe ( λ_j^DMD) )^2 + ( ℐm ( λ_j^DMD) )^2)
    ϕ_j^DMD    = √(( ℛe ( ϕ_j^DMD) )^2 + ( ℐm ( ϕ_j^DMD) )^2)









