Panel flutter analysis

The aerodynamic stiffness matrix of the linear piston theory is calculated with Shell.calc_kA(), as a function of the aerodynamic parameter beta of the Shell object. The flutter onset is found by increasing beta and solving the eigenvalue problem \((\boldsymbol{K}_C + \boldsymbol{K}_G + \boldsymbol{K}_A + \lambda^2 \boldsymbol{M})\boldsymbol{c} = \boldsymbol{0}\) with structsolve.freq(), until two natural frequencies coalesce and the eigenvalues become complex. In this convention \(\boldsymbol{K}_A = -\int \left(\beta (\boldsymbol{S}^w_{,x})^\top \boldsymbol{S}^w + \gamma (\boldsymbol{S}^w)^\top \boldsymbol{S}^w\right) dA\), with the row \(\boldsymbol{S}^w\) of the approximation functions of \(w\), for the load \(q = \beta w_{,x} + \gamma w\) of a flow along \(x\), where the curvature correction \(\gamma\) of Krumhaar softens a cylindrical shell, see Shell.calc_kA(). The term in \(\beta\) is the transpose of that of Eq. (21) of [castro2016FlutterPanel], which for \(w = 0\) at the leading and trailing edges corresponds to reversing the flow direction and does not change the flutter boundary, whereas the sign of the term in \(\gamma\) of Eqs. (21) and (29) of that paper is reversed, a mistake of the paper. The aerodynamic damping, Shell.calc_cA(), is \(i \boldsymbol{C}_A\) with \(\boldsymbol{C}_A\) of Eq. (23) of the paper, entering the eigenvalue problem as \(+\omega \boldsymbol{C}_A\). The code below is extracted from one of the panels unit tests, which verifies the critical beta with and without pre-stress:

def test_aero():
    for model in ['plate_clpt_donnell',
                  'cylshell_clpt_donnell',
                  'cylshell_clpt_sanders',
                  'plate_fsdt_donnell',
                  'plate_tsdt_donnell']:
        for pre_stress in [False, True]:
            print('Flutter Analysis Piston Theory, consider_pre_stress={0}, model={1}'.
                  format(pre_stress, model))
            s = Shell()
            s.model = model
            s.a = 1.
            s.b = 0.5
            s.r = 1.e8
            s.stack = [0, 90, -45, +45]
            s.plyt = 1e-3*0.125
            E2 = 8.7e9
            s.laminaprop = (142.5e9, E2, 0.28, 5.1e9, 5.1e9, 5.1e9)
            s.rho = 1.5e3
            s.m = 8
            s.n = 9

            if pre_stress:
                s.Nxx = -60.
                s.Nyy = -5.
            else:
                s.Nxx = 0
                s.Nyy = 0

            # testing commong methodology based on betastar
            if pre_stress:
                betasstar = np.linspace(50, 350, 40)
            else:
                betasstar = np.linspace(300, 690, 40)
            betas = betasstar/(s.a**3/E2/(len(s.stack)*s.plyt)**3)
            s.beta = betas[0]
            kC = s.calc_kC()
            kG = s.calc_kG()
            kA = s.calc_kA()
            kM = s.calc_kM()

            eigvals, eigvecs = freq(kC + kG + kA, kM, sparse_solver=True, silent=True)
            out = np.zeros((len(betasstar), eigvals.shape[0]),
                    dtype=eigvals.dtype)
            for i, beta in enumerate(betas):
                s.beta = beta
                kA = s.calc_kA()
                eigvals, eigvecs = freq(kC + kG + kA, kM, sparse_solver=True, silent=True)
                omegan = (-eigvals)**0.5
                omegan = omegan*s.a**2/(np.pi**2*sum(s.plyts))*np.sqrt(s.rho/E2)
                out[i, :] = omegan

            ind = np.where(np.any(out.imag != 0, axis=1))[0][0]
            if pre_stress:
                assert np.isclose(betas[ind], 129.663, atol=0.1, rtol=0)
            else:
                assert np.isclose(betas[ind], 358.875, atol=0.1, rtol=0)