NACA airfoil wing model

NACA airfoil wing model

Contents

Finite element model

Model from the following reference:

Jain, S., Tiso, P., Rutzmoser, J. B., & Rixen, D. J. (2017). A quadratic manifold for model order reduction of nonlinear structural dynamics. Computers & Structures, 188, 80�94. <https://doi.org/10.1016/J.COMPSTRUC.2017.04.005>

Finite element code taken from the following package:

Jain, S., Marconi, J., Tiso P. (2020). YetAnotherFEcode (Version v1.1). Zenodo. <http://doi.org/10.5281/zenodo.4011282>

clear all; close all; clc
run ../../install.m
% change to example directory
exampleDir = fileparts(mfilename('fullpath'));
cd(exampleDir)

system parameters

epsilon = 1;

generate model

[M,C,K,fnl,f_0,outdof] = build_model_semiIntrusive();
n = length(M);
disp(['Number of degrees of freedom = ' num2str(n)])
disp(['Phase space dimensionality = ' num2str(2*n)])
Reading mesh from Wing.msh
Building FE model
Assembling M,C,K matrices
Applying boundary conditions
Solving undamped eigenvalue problem

ans = 

  Patch (Deformed Mesh) with properties:

    FaceColor: 'interp'
    FaceAlpha: 1
    EdgeColor: [0 0 0]
    LineStyle: '-'
        Faces: [49968×3 double]
     Vertices: [149904×3 double]

  Use GET to show all properties

Using Rayleigh damping
Assembling external force vector
Getting nonlinearites
Number of degrees of freedom = 133920
Phase space dimensionality = 267840

Dynamical system setup

We consider the forced system

which can be written in the first-order form as

where

.

DSorder = 2;
DS = DynamicalSystem(DSorder);
set(DS,'M',M,'C',C,'K',K,'fnl_semi',fnl);
set(DS.Options,'Emax',5,'Nmax',10,'notation','multiindex')
set(DS.Options,'Intrusion','semi')

% set(DS.Options,'Emax',5,'Nmax',10,'notation','tensor')

We assume periodic forcing of the form

Fourier coefficients of Forcing

kappas = [-1; 1];
coeffs = [f_0 f_0]/2;
DS.add_forcing(coeffs, kappas,epsilon);

Linear Modal analysis and SSM setup

[V,D,W] = DS.linear_spectral_analysis();
Due to high-dimensionality, we compute only the first 5 eigenvalues with the smallest magnitude. These would also be used to compute the spectral quotients
Assuming a proportional damping hypthesis with symmetric matrices
modal damping ratio for 1 mode is 2.000000e-03
modal damping ratio for 2 mode is 2.000000e-03
modal damping ratio for 3 mode is 2.264041e-03
modal damping ratio for 4 mode is 3.612822e-03
modal damping ratio for 5 mode is 5.538531e-03
the left eigenvectors may be incorrect in case of asymmetry of matrices

 The first 10 nonzero eigenvalues are given as 
   1.0e+02 *

  -0.0006 + 0.2934i
  -0.0006 - 0.2934i
  -0.0030 + 1.5226i
  -0.0030 - 1.5226i
  -0.0041 + 1.8088i
  -0.0041 - 1.8088i
  -0.0113 + 3.1381i
  -0.0113 - 3.1381i
  -0.0274 + 4.9385i
  -0.0274 - 4.9385i

Choose Master subspace (perform resonance analysis)

S = SSM(DS);
set(S.Options, 'reltol', 0.1,'notation','multiindex')
% set(S.Options, 'reltol', 0.1,'notation','tensor')
masterModes = [1,2];
S.choose_E(masterModes);
No (near) outer resonances detected in the (truncated) spectrum
sigma_out = 46
(near) inner resonance detected for the following combination of master eigenvalues
     2     1
     3     2
     4     3
     5     4
     1     2
     2     3
     3     4
     4     5

These are in resonance with the follwing eigenvalues of the master subspace
  -0.0587 +29.3428i
  -0.0587 +29.3428i
  -0.0587 +29.3428i
  -0.0587 +29.3428i
  -0.0587 -29.3428i
  -0.0587 -29.3428i
  -0.0587 -29.3428i
  -0.0587 -29.3428i

sigma_in = 46

Forced response curves using SSMs

Obtaining forced response curve in reduced-polar coordinate

order = 3; % Approximation order

setup options

set(S.Options, 'reltol', 1,'IRtol',0.02,'notation', 'multiindex','contribNonAuto',false)
set(S.FRCOptions, 'nt', 2^7, 'nRho', 200, 'nPar', 200, 'nPsi', 100, 'rhoScale', 2 )
% set(S.FRCOptions, 'method','level set')
set(S.FRCOptions, 'method','continuation ep', 'z0', 1e-4*[1; 1])
set(S.FRCOptions, 'outdof',outdof)

choose frequency range around the first natural frequency

omega0 = imag(S.E.spectrum(1));
omegaRange = omega0*[0.9 1.1];

extract forced response curve

FRC = S.extract_FRC('freq',omegaRange,order);
figFRC = gcf;
*****************************************
Calculating FRC using SSM with master subspace: [1  2]
(near) outer resonance detected for the following combination of master eigenvalues
     5     0
     6     0
     6     1
     7     1
     7     2
     8     2
     0     5
     0     6
     1     6
     1     7
     2     7
     2     8
     6     0
     7     0
     7     1
     8     1
     8     2
     0     6
     0     7
     1     7
     1     8
     2     8
    10     0
     0    10

These are in resonance with the follwing eigenvalues of the slave subspace
   1.0e+02 *

  -0.0030 + 1.5226i
  -0.0030 + 1.5226i
  -0.0030 + 1.5226i
  -0.0030 + 1.5226i
  -0.0030 + 1.5226i
  -0.0030 + 1.5226i
  -0.0030 - 1.5226i
  -0.0030 - 1.5226i
  -0.0030 - 1.5226i
  -0.0030 - 1.5226i
  -0.0030 - 1.5226i
  -0.0030 - 1.5226i
  -0.0041 + 1.8088i
  -0.0041 + 1.8088i
  -0.0041 + 1.8088i
  -0.0041 + 1.8088i
  -0.0041 + 1.8088i
  -0.0041 - 1.8088i
  -0.0041 - 1.8088i
  -0.0041 - 1.8088i
  -0.0041 - 1.8088i
  -0.0041 - 1.8088i
  -0.0113 + 3.1381i
  -0.0113 - 3.1381i

sigma_out = 46
(near) inner resonance detected for the following combination of master eigenvalues
     2     1
     3     2
     4     3
     5     4
     1     2
     2     3
     3     4
     4     5

These are in resonance with the follwing eigenvalues of the master subspace
  -0.0587 +29.3428i
  -0.0587 +29.3428i
  -0.0587 +29.3428i
  -0.0587 +29.3428i
  -0.0587 -29.3428i
  -0.0587 -29.3428i
  -0.0587 -29.3428i
  -0.0587 -29.3428i

sigma_in = 46
Due to (near) outer resonance, the exisitence of the manifold is questionable and the underlying computation may suffer.
Attempting manifold computation
Manifold computation time at order 2 = 00:01:59
Estimated memory usage at order  2 = 2.29E+02 MB
Manifold computation time at order 3 = 00:03:14
Estimated memory usage at order  3 = 2.46E+02 MB

 Run='freqSubint1Order3.ep': Continue equilibria along primary branch.

    STEP   DAMPING               NORMS              COMPUTATION TIMES
  IT SIT     GAMMA     ||d||     ||f||     ||U||   F(x)  DF(x)  SOLVE
   0                          2.20e-02  4.17e+01    0.0    0.0    0.0
   1   1  1.00e+00  9.79e-02  1.10e-04  4.17e+01    0.0    0.0    0.0
   2   1  1.00e+00  1.10e-03  4.97e-08  4.17e+01    0.0    0.0    0.0
   3   1  1.00e+00  1.74e-07  9.50e-16  4.17e+01    0.0    0.1    0.0

 STEP      TIME        ||U||  LABEL  TYPE            om         rho1          th1          eps
    0  00:00:00   4.1699e+01      1  EP      2.9343e+01   7.5879e-01   2.7048e+00   1.0000e+00
   10  00:00:00   4.0892e+01      2          2.8745e+01   1.7381e-01   3.0448e+00   1.0000e+00
   20  00:00:01   3.7621e+01      3  EP      2.6409e+01   3.5977e-02   3.1216e+00   1.0000e+00

 STEP      TIME        ||U||  LABEL  TYPE            om         rho1          th1          eps
    0  00:00:01   4.1699e+01      4  EP      2.9343e+01   7.5879e-01   2.7048e+00   1.0000e+00
   10  00:00:01   4.2598e+01      5          3.0015e+01   1.7601e+00   1.6752e+00   1.0000e+00
   14  00:00:01   4.2609e+01      6  SN      3.0030e+01   1.7677e+00   1.5257e+00   1.0000e+00
   14  00:00:01   4.2609e+01      7  FP      3.0030e+01   1.7677e+00   1.5257e+00   1.0000e+00
   20  00:00:01   4.2526e+01      8          2.9988e+01   1.6911e+00   1.2675e+00   1.0000e+00
   30  00:00:01   4.1867e+01      9  FP      2.9587e+01   6.3662e-01   3.6246e-01   1.0000e+00
   30  00:00:01   4.1867e+01     10  SN      2.9587e+01   6.3657e-01   3.6243e-01   1.0000e+00
   30  00:00:01   4.1867e+01     11          2.9587e+01   6.3444e-01   3.6115e-01   1.0000e+00
   40  00:00:01   4.2110e+01     12          2.9766e+01   2.5560e-01   1.4258e-01   1.0000e+00
   50  00:00:01   4.3587e+01     13          3.0812e+01   7.1868e-02   3.9952e-02   1.0000e+00
   55  00:00:01   4.5658e+01     14  EP      3.2277e+01   3.5984e-02   1.9999e-02   1.0000e+00
Time spent on mapping the FRC to physical coordinates maptime 57.0748 seconds. 
Total time spent on FRC computation upto O(3) = 00:06:20