Esta página es para la versión anterior. La página correspondiente en inglés ha sido eliminada en la versión actual.

Análisis modal de un avión flexible de ala voladora

Este ejemplo muestra el cálculo de los modos de flexión de un avión de ala flexible. La respuesta de vibración del ala se recoge en varios puntos a lo largo de su palmo. Los datos se utilizan para identificar un modelo dinámico del sistema. Los parámetros modales se extraen del modelo identificado. Los datos del parámetro modal se combinan con la información de posición del sensor para visualizar los distintos modos de flexión del ala. Este ejemplo requiere cuadro de herramientas de procesamiento de señal™.

El avión flexible del ala

En este ejemplo, estudiamos los datos recopilados de un pequeño avión flexible de ala voladora construido en los Laboratorios de Vehículos Aéreos Deshabitados, Universidad de Minnesota [1]. La geometría de la aeronave se muestra a continuación.

El ala del avión puede sufrir grandes deformaciones bajo carga. Las frecuencias de modo flexible son más bajas que las de los aviones comunes con alas más rígidas. Este diseño flexible reduce los costos de material, aumenta la agilidad y el rango de vuelo de la aeronave. Sin embargo, si no se controlan, los modos flexibles pueden conducir a inestabilidades aeroelásticas catastróficas (flutter). Diseñar leyes de control efectivas para suprimir estas inestabilidades requiere una determinación precisa de los diferentes modos de flexión del ala.

Configuración experimental

El objetivo del experimento es recoger la respuesta de vibración de la aeronave en varios lugares en respuesta a una excitación externa. El avión está suspendido de un marco de madera utilizando un solo resorte en su centro de gravedad. El muelle es lo suficientemente flexible para que la frecuencia natural de la oscilación de masa de resorte no interfiera con las frecuencias fundamentales de la aeronave. Una fuerza de entrada se aplica a través de un agitador electrodinámico Unholtz-Dickie Modelo 20 cerca del centro de la aeronave.

Veinte acelerómetros PCB-353B16 se colocan a lo largo del ala para recoger la respuesta de vibración como se muestra en la siguiente figura.

El comando de entrada del agitador se especifica como una entrada chirp de amplitud constante del formulario

<math display="inline">
<mrow>
<mi mathvariant="italic">A</mi>
<mtext></mtext>
<mi mathvariant="normal">sin</mi>
<mrow>
<mo>(</mo>
<mrow>
<mi>ω</mi>
<mrow>
<mo>(</mo>
<mrow>
<mi mathvariant="italic">t</mi>
</mrow>
<mo>)</mo>
</mrow>
<mi mathvariant="italic">t</mi>
</mrow>
<mo>)</mo>
</mrow>
</mrow>
</math>
. La frecuencia del chirrido varía linealmente con el tiempo, es decir,
<math display="inline">
<mrow>
<mi>ω</mi>
<mrow>
<mo>(</mo>
<mrow>
<mi mathvariant="italic">t</mi>
</mrow>
<mo>)</mo>
</mrow>
<mo stretchy="false">=</mo>
<msub>
<mrow>
<mi mathvariant="italic">c</mi>
</mrow>
<mrow>
<mn>0</mn>
</mrow>
</msub>
<mo stretchy="false">+</mo>
<msub>
<mrow>
<mi mathvariant="italic">c</mi>
</mrow>
<mrow>
<mn>1</mn>
</mrow>
</msub>
<mi mathvariant="italic">t</mi>
</mrow>
</math>
. El rango de frecuencia cubierto por la señal de entrada es de 3 a 35 Hz. Los datos son recogidos por dos acelerómetros (acelerómetros de borde delantero y final en una ubicación) a la vez.x Por lo tanto, se llevan a cabo 10 experimentos para recoger todas las 20 respuestas del acelerómetro. Las mediciones del acelerómetro y del transductor de fuerza se muestrean a 2000 Hz.

Preparación de datos

Los datos están representados por 10 conjuntos de señales de dos salidas/una entrada, cada una de las cuales contiene muestras de 600K. Cargar los datos desde el archivo MAT:FlexibleWingData.mat

load FlexibleWingData.mat MeasuredData

La variable es una estructura con campos , , ..., .MeasuredDataExp1Exp2Exp10 Cada campo es una estructura con campos y que contiene las dos respuestas y los datos de fuerza de entrada correspondientes.y Trazar los datos para el primer experimento.

fs = 2000;                % data sampling frequency Ts = 1/fs;                % sample time y = MeasuredData.Exp1.y;  % output data (2 columns, one for each accelerometer) u = MeasuredData.Exp1.u;  % input force data t = (0:length(u)-1)' * Ts; figure subplot(211) plot(t,y) ylabel('Outputs (m/s^2)') legend('Leading edge','Trailing edge') subplot(212) plot(t,u) ylabel('Input') xlabel('Time (seconds)')

Para preparar los datos para la identificación del modelo, los datos se empaquetan en objetos.iddata Los objetos son una forma estándar de empaquetar datos de dominio de tiempo en System Identification Toolbox™.iddata La señal de entrada se trata como limitada por bandas.

ExpNames = fieldnames(MeasuredData); Data = cell(1, 10); for k = 1:10    y =  MeasuredData.(ExpNames{k}).y;    u =  MeasuredData.(ExpNames{k}).u;    Data{k} = iddata(y, u, Ts, 'InterSample', 'bl'); end

Combine los objetos de conjunto de datos en un objeto de datos multiexperimento.

Data = merge(Data{:})
Data =  Time domain data set containing 10 experiments.  Experiment   Samples      Sample Time                Exp1      600001            0.0005                Exp2      600001            0.0005                Exp3      600001            0.0005                Exp4      600001            0.0005                Exp5      600001            0.0005                Exp6      600001            0.0005                Exp7      600001            0.0005                Exp8      600001            0.0005                Exp9      600001            0.0005                Exp10     600001            0.0005                                                               Outputs      Unit (if specified)                     y1                                                y2                                                                                               Inputs       Unit (if specified)                     u1                                                                                               

Identificación del modelo

Queremos identificar un modelo dinámico cuya respuesta de frecuencia coincida con la de la aeronave real lo más cerca posible. Un modelo dinámico encapsula una relación matemática entre las entradas y salidas del sistema como una ecuación diferencial o de diferencia. Ejemplos de modelos dinámicos son las funciones de transferencia y los modelos de espacio de estado. En System Identification Toolbox, los modelos dinámicos se encapsulan mediante objetos como (para funciones de transferencia), (para modelos AR, ARMA) y (para modelos de espacio de estado).idtfidpolyidss Los modelos dinámicos se pueden crear ejecutando comandos de estimación, como y comandos en datos medidos, ya sea en el dominio de tiempo o en el dominio de frecuencia.tfestssest

En este ejemplo, primero convertimos los datos de dominio de tiempo medidos en datos de respuesta de frecuencia mediante estimación de la función de transferencia empírica mediante el comando.etfe El FRF estimado se utiliza entonces para identificar un modelo de espacio estatal de la respuesta de vibración de la aeronave. Es posible utilizar directamente los datos de dominio de tiempo para la identificación dinámica del modelo. Sin embargo, la forma de datos FRF permite la compresión de grandes conjuntos de datos en menos muestras, así como ajustar más fácilmente el comportamiento de estimación a los rangos de frecuencia relevantes. Los FRF se encapsulan mediante objetos.idfrd

Calcule una función de respuesta de frecuencia de dos salidas/una entrada (FRF) para cada experimento de datos. No use ventanas. Utilice 24.000 puntos de frecuencia para el cálculo de la respuesta.

G = cell(1, 10); N = 24000; for k = 1:10    % Convert time-domain data into a FRF using ETFE command    Data_k = getexp(Data, k);    G{k} = etfe(Data_k, [], N); % G{k} is an @idfrd object end

Concatene todos los FRF en un solo FRF de 20 salidas/una entrada.

G = cat(1, G{:});     % concatenate outputs of all estimated FRFs G.OutputName = 'y';   % name outputs 'y1', 'y2', ..., 'y20' G.InterSample = 'bl';

Para obtener una idea de la respuesta de frecuencia estimada, trace la amplitud de las salidas 1 y 15 (seleccionadas arbitrariamente). Amplíe el rango de frecuencia de interés (4–35 Hz).

opt = bodeoptions;           % plot options opt.FreqUnits = 'Hz';        % show frequency in Hz opt.PhaseVisible = 'off';    % do not show phase OutputNum = [1 15];          % pick outputs 1 and 15 for plotting clf bodeplot(G(OutputNum, :), opt)   % plot frequency response xlim([4 35]) grid on

El FRF muestra al menos 9 frecuencias resonantes. Para el análisis queremos centrarnos en el intervalo de frecuencia de 6-35 Hz donde se encuentran los modos de flexión flexibles más críticos de la aeronave. Por lo tanto, reduzca el FRF a esta región de frecuencia.

f = G.Frequency/2/pi;           % extract frequency vector in Hz (G stores frequency in rad/s) Gs = fselect(G, f>6 & f<=32)    % "fselect" selects the FRF in the requested range (6.5 - 35 Hz)
Gs = IDFRD model. Contains Frequency Response Data for 20 output(s) and 1 input(s). Response data is available at 624 frequency points, ranging from 37.96 rad/s to 201.1 rad/s.   Sample time: 0.0005 seconds Output channels: 'y(1)', 'y(2)', 'y(3)', 'y(4)', 'y(5)', 'y(6)', 'y(7)', 'y(8)', 'y(9)', 'y(10)', 'y(11)', 'y(12)', 'y(13)', 'y(14)', 'y(15)', 'y(16)', 'y(17)', 'y(18)', 'y(19)', 'y(20)' Input channels: 'u1' Status:         Estimated model 

por lo tanto, contiene las mediciones de respuesta de frecuencia en las 20 ubicaciones de medición.Gs A continuación, identifique un modelo de espacio de estado que se ajuste.Gs El estimador subespacial proporciona una estimación no iterativa rápida.n4sid La estructura del modelo de espacio de estado se configura de la siguiente manera:

  1. Estimamos un modelo de tiempo continuo de 18o orden. La orden se encontró después de algunas pruebas con varias órdenes y comprobar el ajuste resultante del modelo a la FRF.

  2. El modelo contiene un término de avance (matriz no es cero).D Esto se debe a que estamos limitando nuestro análisis a

    <math display="inline">
    <mrow>
    <mo></mo>
    </mrow>
    </math>
    35 Hz mientras que el ancho de banda del ala es significativamente mayor que eso (la respuesta es significativa en 35 Hz).

  3. Para acelerar el cálculo, suprimimos el cálculo de la covarianza de parámetros.

  4. La magnitud del FRF varía significativamente a lo largo del rango de frecuencia. Con el fin de garantizar que las amplitudes bajas reciban el mismo énfasis que las amplitudes más altas, aplicamos una ponderación personalizada que varía inversamente como la raíz cuadrada de la 11a respuesta. La elección de la salida 11 es algo arbitraria, pero funciona ya que los 20 FRF tienen perfiles similares.

Configuramos las opciones de estimación para usar .n4sidn4sidOptions

FRF = squeeze(Gs.ResponseData); Weighting = mean(1./sqrt(abs(FRF))).'; n4Opt = n4sidOptions('EstimateCovariance',false,...    'WeightingFilter',Weighting,...    'OutputWeight',eye(20)); sys1 = n4sid(Gs,24,'Ts',0,'Feedthrough',true,n4Opt); Fit = sys1.Report.Fit.FitPercent'
Fit = 1×20

   57.0200   57.9879   85.0160   86.3815   80.4879   80.4430   58.2216   45.2692   61.5057   76.7612   84.7305   86.2600   86.4266   85.0251   76.9208   82.1191   74.7982   79.6837   67.9078   76.7249

Los números "Fit" muestran el porcentaje de ajuste entre los datos ( ) y la respuesta de frecuencia del modelo ( ) utilizando una medida de bondad de ajuste de raíz media cuadrada (NRMSE) normalizada.Gssys1 Los ajustes más pobres y mejores se trazan a continuación.

[~,Imin] = min(Fit); [~,Imax] = max(Fit); clf bodeplot(Gs([Imin, Imax],:), sys1([Imin, Imax],:), opt); xlim([6 32]) title('Worst (top) and best (bottom) fits between data and model') grid on legend('Data', 'Model')

Los ajustes logrados con el modelo se pueden mejorar aún más mediante el refinamiento iterativo no lineal de los mínimos cuadrados de los parámetros del modelo.sys1 Esto se puede lograr usando el comando.ssest Configuramos las opciones de estimación para usar el comando.ssestssestOptions Esta vez también se estima la covarianza del parámetro.

ssOpt = ssestOptions('EstimateCovariance',true,...    'WeightingFilter',n4Opt.WeightingFilter,...    'SearchMethod','lm',...     % use Levenberg-Marquardt search method    'Display','on',...    'OutputWeight',eye(20)); sys2 = ssest(Gs, sys1, ssOpt);  % estimate state-space model (this takes several minutes) Fit = sys2.Report.Fit.FitPercent'
Fit = 1×20

   89.7225   89.5185   89.7260   90.4986   88.5522   88.8727   81.3225   83.5975   75.9215   83.1763   91.1358   89.7310   90.6844   89.8444   89.6685   89.1467   87.8532   88.0385   84.2898   83.3578

Como antes de trazar lo peor y lo mejor. También visualizamos la región de confianza 1-sd para la respuesta de frecuencia del modelo.

[~, Imin] = min(Fit); [~, Imax] = max(Fit); clf h = bodeplot(Gs([Imin, Imax],:), sys2([Imin, Imax],:), opt); showConfidence(h, 1) xlim([6.5 35]) title('Worst (top) and best (bottom) fits between data and refined model') grid on legend('Data', 'Model (sys2)')

El refinado modelo de espacio-estado se adapta bastante bien a los FRF en la región de 7-20 Hz.sys2 Los límites de incertidumbre son estrechos alrededor de la mayoría de los lugares resonantes. Estimamos un modelo de orden 24, lo que significa que hay como máximo 12 modos oscilatorios en el sistema.sys2 Utilice el comando para capturar las frecuencias naturales en Hz para estos modos.modalfit

f = Gs.Frequency/2/pi;     % data frequencies (Hz) fn = modalfit(sys2, f, 12); % natural frequencies (Hz) disp(fn)
    7.7721     7.7953     9.3147     9.4009     9.4910    15.3463    19.3291    23.0219    27.4164    28.7256    31.7014    63.3034 

Los valores en sugieren dos frecuencias muy cerca de 7.8 Hz y tres cerca de 9.4 Hz.fn Una inspección de las respuestas de frecuencia cerca de estas frecuencias revela que la ubicación de los picos cambia un poco a través de las salidas. Estas discrepancias pueden eliminarse mediante un mejor control sobre el proceso del experimento, realizando una identificación directa de dominio de tiempo con ancho de banda de entrada limitado a un rango estrecho centrado en estas frecuencias, o ajustando un único modo oscilatoico a la frecuencia respuesta en torno a estas frecuencias. Estas alternativas no se exploran en este ejemplo.

Cálculo de parámetros modales

Ahora podemos usar el modelo para extraer los parámetros modales.sys2 Una inspección de los FRF indica alrededor de 10 frecuencias modales, aproximadamente alrededor de las frecuencias [5 7 10 15 17 23 27 30] Hz. Podemos hacer esta evaluación más concreta mediante el comando que comprueba la estabilidad de los parámetros modales a medida que se cambia el orden del modelo subyacente.modalsd Esta operación lleva mucho tiempo. Por lo tanto, el trazado resultante se inserta directamente como una imagen. Ejecute el código dentro del bloque de comentarios a continuación para reproducir la figura.

FRF = permute(Gs.ResponseData,[3 1 2]); f = Gs.Frequency/2/pi; %{  figure  pf = modalsd(FRF,f,fs,'MaxModes',20,'FitMethod','lsrf'); %}

La inspección de la gráfica y los valores sugiere una lista refinada de frecuencias naturales reales:pf

Freq = [7.8 9.4 15.3 19.3 23.0 27.3 29.2 31.7];

Usamos los valores en vector como guía para elegir los modos más dominantes del modelo.Freqsys2 Esto se hace usando el comando.modalfit

[fn,dr,ms] = modalfit(sys2,f,length(Freq),'PhysFreq',Freq);

son las frecuencias naturales (en Hz), los coeficientes de amortiguación correspondientes y los residuos normalizados como vectores de columna (una columna para cada frecuencia natural).fndrms En el proceso de extracción de estos parámetros modales, sólo se utilizan los polos estables y subamplificados del modelo. Las columnas contienen datos solo para los polos con parte imaginaria positiva.ms

Visualización de formas de modo

Para visualizar los diferentes modos de flexión, necesitamos los parámetros modales extraídos anteriormente. Además necesitamos información sobre la ubicación de los puntos de medición. Estas posiciones (valores x-y) se registran para cada acelerómetro en la matriz:AccePos

AccelPos = [...   % see figure 2    16.63 18.48;   % nearest right of center    16.63 24.48;    27.90 22.22;    27.90 28.22;    37.17 25.97;    37.17 31.97;    46.44 29.71;    46.44 35.71;    55.71 33.46;    55.71 39.46;  % farthest right    -16.63 18.48; % nearest left of center    -16.63 24.18;    -27.90 22.22;    -27.90 28.22;    -37.17 25.97;    -37.17 31.97;    -46.44 29.71;    -46.44 35.71;    -55.71 33.46;    -55.71 39.46]; % farthest left

Las formas de modo están contenidas en la matriz donde cada columna corresponde a la forma para una frecuencia.ms Animar el modo superponiendo amplitudes de forma de modo sobre las coordenadas del sensor y variando las amplitudes a la frecuencia natural del modo. La animación muestra la flexión sin amortiguación. Como ejemplo, considere el modo a 15,3 Hz.

AnimateOneMode(3, fn, ms, AccelPos); 

Conclusiones

En este ejemplo se muestra un enfoque basado en la identificación del modelo paramétrico para la estimación de parámetros modales. Utilizando un modelo de espacio de estado de orden 24, se extrajeron 8 modos oscilatorcidos estables en la región de frecuencia 6-32 Hz. La información modal se combinó con el conocimiento de las posiciones del acelerómetro para visualizar los diversos modos de flexión. Algunos de estos modos se muestran en las figuras siguientes.

Referencias

[1] Gupta, Abhineet, Peter J. Seiler y Brian P. Danowsky. "Pruebas de vibración terrestre en un avión flexible de ala voladora." (AIAA 2016-1753).AIAA Atmospheric Flight Mechanics Conference, AIAA SciTech Forum.

function AnimateOneMode(ModeNum, fn, ModeShapes, AccelPos) % Animate one mode. % ModeNum: Index of the mode.      % Reorder the sensor locations for plotting so that a continuous,  % non-intersecting curve is traced around the body of the aircraft. PlotOrder = [19:-2:11, 1:2:9, 10:-2:2, 12:2:20, 19]; Fwd = PlotOrder(1:10); Aft = PlotOrder(20:-1:11);  x = AccelPos(PlotOrder,1); y = AccelPos(PlotOrder,2); xAft = AccelPos(Aft,1); yAft = AccelPos(Aft,2); xFwd = AccelPos(Fwd,1); yFwd = AccelPos(Fwd,2);  wn = fn(ModeNum)*2*pi;  % Mode frequency in rad/sec T = 1/fn(ModeNum);      % Period of modal oscillation Np  = 2.25;             % Number of periods to simulate tmax = Np*T;            % Simulate Np periods Nt = 100;               % Number of time steps for animation t = linspace(0,tmax,Nt); ThisMode = ModeShapes(:,ModeNum)/max(abs(ModeShapes(:,ModeNum))); % normalized for plotting z0 = ThisMode(PlotOrder); zFwd = ThisMode(Fwd); zAft = ThisMode(Aft);  clf col1 = [1 1 1]*0.5; xx = reshape([[xAft, xFwd]'; NaN(2,10)],[2 20]); yy = reshape([[yAft, yFwd]'; NaN(2,10)],[2 20]); plot3(x,y,0*z0,'-', x,y,0*z0,'.', xx(:), yy(:), zeros(40,1),'--',...    'Color',col1,'LineWidth',1,'MarkerSize',10,'MarkerEdgeColor',col1); xlabel('x') ylabel('y') zlabel('Amplitude') ht = max(abs(z0)); axis([-100  100  10  40 -ht ht]) view(5,55) title(sprintf('Mode %d. FN = %s Hz', ModeNum, num2str(fn(ModeNum),4)));  % Animate by updating z-coordinates.   hold on col2 = [0.87 0.5 0]; h1 = plot3(x,y,0*z0,'-', x,y,0*z0,'.', xx(:), yy(:), zeros(40,1),'--',...    'Color',col2,'LineWidth',1,'MarkerSize',10,'MarkerEdgeColor',col2); h2 = fill3(x,y,0*z0,col2,'FaceAlpha',0.6); hold off  for k = 1:Nt    Rot1 = cos(wn*t(k));    Rot2 = sin(wn*t(k));    z_t = real(z0)*Rot1 - imag(z0)*Rot2;    zAft_t = real(zAft)*Rot1 - imag(zAft)*Rot2;    zFwd_t = real(zFwd)*Rot1 - imag(zFwd)*Rot2;    zz = reshape([[zAft_t, zFwd_t]'; NaN(2,10)],[2 20]);    set(h1(1),'ZData',z_t)    set(h1(2),'ZData',z_t)    set(h1(3),'ZData',zz(:))    h2.Vertices(:,3) = z_t;    pause(0.1) end  end