Version 1.0

AeroElastic analysis using Dynamic Mode Interpolation

Methodology

aedmi computes the solution of the aeroelastic equation u/l * M * p^2 + K - 1/2*r*u^2 * Q(k) = 0, where M and K are the modal mass and stiffness matrices of the structure, Q is the modal aerodynamic loads matrix, u and r are the freestream fluid velocity and density, and l is a reference length. p is a non-dimensional parameter defined as p = (g+1i)*k, where g is the true damping and k is the reduced frequency.
The flutter solution is obtained using a kind of reduced order modeling technique called dynamic mode interpolation. This technique was originally developed by Hüseyin Güner during his doctoral thesis. Practically, the mode shapes of the structure are first pre-computed. Then, the aerodynamic loads at some reference Mach numbers and reduced frequencies are pre-computed by imposing the motion for each structural mode shape. Finally, aedmi will compute the modal load matrix Q for each Mach number and reduced frequency, and solve the aeroelastic system by interpolating this matrix for any reduced frequency.

Dependencies

• python 3 with numpy and scipy
• [optional] VTK python bindings to read VTK formatted data
• [optional] matplotlib to save figures

Code architecture

The code is split into several submodules:

• iops: classes to read and write from/to disk
• msh: simple mesh data structure
• flutter: classes to perform flutter analysis
• utils: various general utilities

The flutter (core) submodule contains:

• structural: read mode shapes and structural matrices from disk
• aero: read unsteady aerodynamic pressures from disk and compute reduced aerodynamic load matrix
• fluid: compute fluid density and velocity from Mach number using a matched or unmatched strategy
• solution: compute the flutter speed

Issues

Several issues are opened (#1, #2, #3), but none is critical.

Tests

Battery successfully passes on ubuntu 18.04, python 3.6 and on windows 10, python 3.8.

flutter/aero.py

+#!/usr/bin/env python3
+# -*- coding: utf8 -*-
+# test encoding: à-é-è-ô-ï-€
+
+# Copyright 2021 University of Liège
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+# TODO Add interpolation from structural grid (mode shapes) to aerodynamic grid
+
+import time
+import numpy as np
+
+class Aero:
+    """Class for computing aerodynamic loads
+    """
+    def __init__(self, machs, ks):
+        self.machs = machs # list of mach numbers
+        self.ks = ks # list of reference reduced frequencies
+        # Note: the code can read data from several time instances and produce the same number of flow modes
+        #       however, only the first mode can be currently used to compute the modal load matrix Q
+        #       as such, self.nh is constrained to 1, so that 3 time instances are read, and 3 modes (0th, 1st and its complex conjugate) are computed (with only the 1st being used)
+        self.nh = 1 # number of retained harmonics
+    def __str__(self):
+        return 'Aerodynamic Loads with\nMach numbers: {}\nReduced frequencies: {}'.format(self.machs, self.ks)
+
+    def __loads(self, grid, hcp, ndim, struct):
+        """Compute modal loads matrix
+            grid: aerodynamic surface grid
+            hcp: modal pressure coefficients of 1st mode
+            ndim: problem dimension
+            struct: structural model
+           Q(nmodes, nmodes) = W(nmodes, npts) * N(npts, nmodes)
+           where W contains the vertical component of the mode shapes, and N contains the vertical component of the normal forces (divided by the freestream dynamic pressure)
+        """
+        # Timer
+        print('Computing  modal load matrix...', end=' ', flush=True)
+        cpu = time.perf_counter()
+        # Compute averaged modal displacements at cell
+        nc = len(grid.cells) # number of cells
+        w = np.zeros((struct.ns, nc))
+        for k in range(struct.ns):
+            for l in range(nc):
+                cell = grid.cells[l]
+                for m in range(len(cell.nodes)):
+                    w[k,l] += struct.modes[k][cell.nodes[m].idx] / len(cell.nodes)
+            w[k,:] *= struct.pfs[k] / struct.amps[k] # mode shapes (w_z * participation / amplification)
+        # Compute cell area and unit normal vector
+        a = np.zeros(nc)
+        n = np.zeros((nc,3))
+        for i in range(nc):
+            n[i,:], a[i] = grid.cells[i].evalNormalArea()
+        # Compute loads at cell from averaged nodal pressure
+        qs = [[np.zeros((struct.ns, struct.ns), dtype=np.complex64) for _ in range(len(self.ks))] for _ in range(len(self.machs))]
+        for i in range(len(self.machs)):
+            for j in range(len(self.ks)):
+                hcn = np.zeros((nc, struct.ns), dtype=np.complex64)
+                for k in range(struct.ns):
+                    for l in range(nc):
+                        cell = grid.cells[l]
+                        for m in range(len(cell.nodes)):
+                            hcn[l,k] += hcp[i][j][k][cell.nodes[m].idx] / len(cell.nodes)
+                    hcn[:,k] *= -a * n[:,ndim-1] * struct.pfs[k] / struct.amps[k] # normalized load (cn_z = -A*n_z*cp * participation / amplification)
+                qs[i][j] = 1j * 2 * np.matmul(w, hcn)  # load matrices at reference frequencies and Mach numbers
+        print('done! Wall-clock time:', time.perf_counter() - cpu, 's')
+        return qs
+
+    def fromTime(self, reader, cpname, ndim, struct):
+        """Read pressure coefficient from time domain and convert it to frequency domain to compute modal load matrix
+            reader: reader used to read data from disk
+            cpname: name of the pressure coefficient variable
+            ndim: problem dimension
+            struct: structural model
+        """
+        # Compute direct Fourier transmform (time -> frequency)
+        nf = 2 * self.nh + 1 # number of flow modes
+        harmonics = np.zeros(nf) # sequence of harmonic factors (3 modes: [0, 1, -1])
+        for i in range(self.nh):
+            harmonics[2*(i+1)-1] = i+1
+            harmonics[2*(i+1)] = -(i+1)
+        instances = 2 * np.pi * np.arange(0, nf, 1) / nf # sequence of instances (3 modes: 2pi/N*[0, 1, 2])
+        dft = 1 / nf * np.exp(-1j * np.outer(harmonics, instances)) # DFT[k,n] = 1/N * exp(-j*2pi/N*k*n)
+        # Read the aerodynamic grid
+        grid = reader.readGrid('mach{0:03d}_freq{1:03d}_mode{2:03d}_sol{3:1d}'.format(0, 0, 0, 0))
+        # Read pressure coefficient at nf time instances
+        hcp = [[[np.zeros(len(grid.nodes)) for _ in range(struct.ns)] for _ in range(len(self.ks))] for _ in range(len(self.machs))]
+        for i in range(len(self.machs)):
+            for j in range(len(self.ks)):
+                for k in range(struct.ns):
+                    cp = np.zeros((len(grid.nodes),nf)) # time data
+                    for l in range(nf):
+                        cp[:,l] = np.squeeze(reader.readValues('mach{0:03d}_freq{1:03d}_mode{2:03d}_sol{3:1d}'.format(i, j, k, l), [cpname])[cpname])
+                    # Convert from time to frequency domain and retain first mode
+                    hcp[i][j][k] = np.matmul(cp, np.transpose(dft))[:,1]
+        # Compute modal load matrix
+        self.Q = self.__loads(grid, hcp, ndim, struct)
+
+    def fromFreq(self, reader):
+        """Read data from frequency domain to compute modal load matrix
+            reader: reader used to read data from disk
+        """
+        # Read the aerodynamic grid
+        #grid = reader.readGrid('mach{0:03d}_freq{1:03d}_mode{2:03d}_sol0'.format(0, 0, 0))
+        # Read modal pressure coefficient from first mode
+        #hcp[i][j][k] = np.squeeze(reader.readValues('mach{0:03d}_freq{1:03d}_mode{2:03d}_sol1'.format(i, j, k), [cpname])[cpname])
+        # Compute modal load matrix
+        #self.Q = self.__loads(...)
+        raise NotImplementedError()
+
+    def compute(self, m, k):
+        """Compute aerodynamic load matrix by linear interpolation
+            m: Mach number
+            k: reduced frequency
+           Q = (k2 - k)/(k2-k1)*Q(k1) + (k - k1)/(k2-k1)*Q(k2)
+        """
+        return (self.Q[m][1] - self.Q[m][0]) / (self.ks[1] - self.ks[0]) * (k - self.ks[0]) + self.Q[m][0]