Problema de Flameo Aeroelástico con método de Galerkin

En este artículo voy a presentar el desarrollo y solución de un problema de Aeroelasticidad (específicamente el fenómeno de Flameo) correspondiente a un examen parcial de la asignatura homónima del Máster de Ingeniería Aeronáutica de la Universidad de Sevilla. A continuación se presenta el enunciado:

Enunciado de problema aeroelástico de flameo
Enunciado de problema aeroelástico de flameo

Datos adicionales

Para las características de la placa se emplea un aluminio 6061-T6 cuyas propiedades son las siguientes:

\rho = 2700 \ kg/m^3

E = 68.9 \ GPa

G = 26 \ GPa

El momento de inercia de área y el módulo de torsión se calcula, respectivamente, como:

I = \dfrac{1}{12} 2bt^{3}

J = \dfrac{1}{3} 2bt^{3}

La masa e inercia con respecto al eje y por unidad de envergadura se calcula, respectivamente, como:

m = 2 b t \rho

i_{yG} = \dfrac{1}{12} m (2 b)^{2}

Finalmente se toma la densidad del aire en condiciones normales:

\rho_{\infty} = 1.225 kg/m^{3}

Planteamiento del problema

Para resolver el problema emplearemos el modelo de barra continua como aproximación. De esta manera, analizando una rebanada, podremos calcular los esfuerzos en la sección.

Una barra tiene infinitos grados de libertad. Sin embargo, cada sección presenta 6 grados de libertad (3 traslaciones y 3 rotaciones). No obstante, debido a las condiciones de contorno del problema y las fuerzas externas, sólo son no nulos los desplazamientos verticales en el eje z y los giros alrededor del eje x. Los esfuerzos los concentramos en el eje elástico, que es el punto en el que se aplica la Torsión, Flexión y Axil.

Como la placa tiene densidad constante la posición del centro de gravedad debe encontrarse a una distancia x_{G} = b del borde de ataque, y debido a las condiciones de contorno podemos decir que el eje elástico tambièn se encuentra en x_{e} = b.

Definimos el desplazamiento vertical del eje elástico positivo con el eje z como una función que depende de la sección y el tiempo w(y,t) = w_{G}(y,t) y el giro alrededor de dicho eje positivo con el eje y de igual forma \alpha(y,t) donde dicho giro coincide con el ángulo de ataque de la placa.

A continuación se muestran los esfuerzos aplicados en la rebanada, considerando la sustentación aplicada en el centro de gravedad junto con un momento aerodinámico en dicho punto, las inercias en desplazamiento vertical y rotación alrededor del eje y, y el peso:

Rebanada de la placa con fuerzas y esfuerzos internos
Rebanada de la placa con fuerzas y esfuerzos internos

El sumatorio de fuerzas y momentos nos proporcionan las siguientes ecuaciones:

\displaystyle \sum F_{z} = 0 \rightarrow \cancel{V} + dV - \cancel{V} + l \cdot dy - mg \cdot dy - m \ddot w_{G} \cdot dy = 0

\displaystyle \sum  M_{x}^{P} = 0 \rightarrow \cancel{M} + dM - \cancel{M} + V \cdot dy + \cancelto{0}{ (\dots) \cdot \dfrac{dy^{2}}{2} }= 0

\displaystyle \sum  M_{y}^{P} = 0 \rightarrow \cancel{T}  + dT - \cancel{T}  + m_{ac} \cdot dy - i_{yG} \ddot \alpha \cdot dy = 0

Obsérvese que en el sumatorio de momentos alrededor del eje x hemos despreciado los momentos generados por las cargas distribuidas al ser de un orden de magnitud muy inferior (diferencial al cuadrado) así como que no hemos considerado la inercia de cada sección alrededor de dicho eje puesto que se asume que la energía cinética asociada al desplazamiento vertical de la sección es mucho mayor que la asociada al giro.

De la segunda ecuación obtenemos:

\boxed {V = - \dfrac{dM}{dy}}

Introduciendo dicho resultado en la primera ecuación obtenemos:

\dfrac{dV}{dy} = \boxed{-  \dfrac{d^{2}M}{dy^{2}} =  mg  + m \ddot w_{G} - l}

La última ecuación resulta:

\boxed{\dfrac{dT}{dy} = - m_{ac}+ i_{yG} \ddot \alpha}

El modelo de barras nos proporciona una relación entre los esfuerzos y las deformaciones. Despreciando la deformación por cortante e incluyendo los datos de la placa (módulos e inercia) tenemos la siguiente relación:

\dfrac{M}{EI} = \dfrac{d^{2} w}{dy^{2}}

\dfrac{T}{GJ} = \dfrac{d \alpha}{dy}

Sustituyendo los esfuerzos obtenemos las ecuaciones en derivadas parciales para los dos grados de libertad:

EI \dfrac{\partial^{4} w (y,t)}{\partial y^{4}} + mg + m \dfrac{\partial ^{2} w (y,t)}{\partial t^{2}} - l = 0

GJ \dfrac{\partial ^{2} \alpha}{\partial y^{2}} + m_{ac} - i_{yG} \dfrac{\partial ^{2} \alpha}{\partial t^{2}} = 0

Para resolver las ecuaciones se empleará el método de Galerkin. Definimos cada variable w(y,t), \alpha(y,t) como un sumatorio de funciones independientes tal que:

\displaystyle w(y,t) = \sum_{i=1}^{n} F_{i}(y) \cdot W_{i}(t)

\displaystyle \alpha(y,t) = \sum_{i=1}^{m} H_{i}(y) \cdot A_{i}(t)

Donde las funciones F_{i}(y),H_{i}(y) representan el modo de deformación que adquiere la barra en cada sección y W_{i}(t),A_{i}(t) la magnitud de cada uno de los modos de deformación en cada instante de tiempo.

Condiciones de contorno

En y=0

La placa en dicho extremo se encuentra apoyada en muelles en las esquinas. Sustituyendo los muelles por la reacción y visualizando los esfuerzos internos se tiene el siguiente esquema:

Problema de Flameo Aeroelástico con método de Galerkin 1

Los desplazamientos en las esquinas se componen de desplazamiento del eje elástico + giro de la sección:

w_{1} = w(0,t) + \alpha(0,t) \cdot b

w_{1} = w(0,t)-  \alpha(0,t) \cdot b

1) Flector nulo

M(0,t) = 0 = EI \dfrac{d^{2} w(0,t)}{dy^{2}} \iff w''(0,t) = 0 \iff \boxed{ F_{i}''(0) = 0}

2) Cortante equivalente

V(0,t) = k (w_{1} + w_{2}) = 2kw(0,t) = -EI w'''(0,t) \iff \boxed{F_{i}(0) = -\dfrac{EI}{2k} F_{i}'''(0)}

3) Torsor equivalente

T(0,t) = k w_{1} b - k w_{2} b = 2 k b^{2} \alpha(0,t) = GJ \alpha'(0,t) \iff \boxed{H_{i}(0) = \dfrac{GJ}{2kb^{2}} H_{i}'(0)}

En y=s

La placa en dicho extremo se encuentra apoyada en el eje elástico. Sustituyendo el apoyo por la reacción se tiene el siguiente esquema:

Problema de Flameo Aeroelástico con método de Galerkin 2

1) Desplazamiento nulo

w(s,t) = 0 \iff \boxed{F_{i}(s) = 0}

2) Flector nulo

M(s,t) = 0 = EI \dfrac{d^{2} w(s,t)}{dy^{2}} \iff w''(s,t) = 0 \iff \boxed{ F_{i}''(s) = 0}

3) Torsor nulo

T(s,t) = 0 = GJ \alpha'(s,t) \iff \boxed{H_{i}'(s) = 0}

Funciones de aproximación

Como funciones de aproximación se emplearan funciones polinómicas tal que:

 F_{i} = a_{0} + a_{1} \left (\dfrac{y}{s} \right ) +a_{2} \left(\dfrac{y}{s} \right)^{2} + a_{3} \left(\dfrac{y}{s} \right)  ^{3} +  \left(\dfrac{y}{s} \right)^{3+i}

Con i \in [1,n]

H_{i} = b_{0} + b_{1} \left(\dfrac{y}{s} \right) +\left (\dfrac{y}{s} \right)  ^{1+i}

Con i \in [1,m]

Sustituyendo dichas funciones en las condiciones de contorno (6 incógnitas, 6 condiciones) podemos obtener el valor de las constantes:

a_{0} = \dfrac{3EI}{ks^{3}} \cdot \dfrac{(3+i)\cdot (2+i)}{6}

a_{1} = \dfrac{(3+i)(2+i)(ks^3-3EI) - 6ks^3}{6ks^{3}}

a_{2} = 0

a_{3} = - \dfrac{(3+i)\cdot(2+i)}{6}

b_{0} = -\dfrac{GJ(1+i)}{2kb^{2}s}

b_{1} = -(1+i)

Modelo aerodinámico

La solución temporal a dichas ecuaciones diferenciales debe ser una función exponencial compleja tal que:

W_{i}(t) = \bar{W_{i}} e^{\nu t}

A_{i}(t) = \bar{A_{i}} e^{\nu t}

Donde el exponente es un número complejo tal que:

\nu = \gamma + i \cdot \omega

Para definir el modelo aerodinámico y así poder resolver las ecuaciones, debemos emplear la aproximación de Theodorsen que define la sustentación y el momento aerodinámico en el eje elástico en el caso de que el movimiento fuera oscilatorio sin término de amortiguamiento (es decir, w = \bar{w} e^{i \cdot \omega}, \alpha= \bar{\alpha} e^{i \cdot \omega} ) y usar dichos valores en el problema completo.

Los valores de sustentación y momento aerodinámico en el eje elástico según Theodorsen resultan:

l=\pi \rho_{\infty} b^{3} \omega^{2} [l_{w}(k) \cdot \dfrac{\bar{w}}{b} + l_{\alpha}(k) \bar{\alpha}]

m_{ac}=\pi \rho_{\infty} b^{4} \omega^{2} [m_{w}(k)\cdot\dfrac{\bar{w}}{b} + m_{\alpha}(k) \bar{\alpha}]

Cuyos coeficientes se extraen directamente del libro de Hodges [Dewey H Hodges, G Alvin Pierce – Introduction to Structural Dynamics and Aeroelasticity 2011, página 216 ecuación 5.134]:

l_{w} = 1 - 2i\dfrac{C(k)}{k}

l_{\alpha} = a + \dfrac{i}{k} \left  [1+2(1/2 + a) \cdot C(k)  \right ] + \dfrac{2 \cdot C(k)}{k^{2}}

m_{w} = a - \dfrac{2i(1/2+a) \cdot C(k)}{k}

m_{\alpha} = a^{2} + \dfrac{1}{8} - \dfrac{(1/2-a)[1-2(1/2+a) \cdot C(k)] \cdot i}{k} + \dfrac{2(1/2 + a) \cdot C(k)}{k^{2}}

Con la función de Theodorsen definida como:

C(k) = \dfrac{K_{1}(ik)}{K_{0}(ik)+K_{1}(ik)}

Donde K_{0},K_{1} son las Funciones de Bessel modificadas de segunda especie y k la frecuencia reducida.

Sustitución ecuaciones

Recordemos que las ecuaciones a sustituir son las siguientes:

EI \dfrac{\partial^{4} w (y,t)}{\partial y^{4}} + mg + m \dfrac{\partial ^{2} w (y,t)}{\partial t^{2}} - l = 0

GJ \dfrac{\partial ^{2} \alpha}{\partial y^{2}} + m_{ac} - i_{yG} \dfrac{\partial ^{2} \alpha}{\partial t^{2}} = 0

Sustituyendo los grados de libertad por el sumatorio de funciones que hemos definido, el modelo aerodinámico por el modelo según Theodorsen e integrando las ecuaciones por cada uno de los modos de deformación (para así aplicar Galerkin correctamente y repartir el error) obtenemos las siguientes ecuaciones integrales:

\displaystyle \int_{0}^{s} \left [EI \sum_{j=1}^{n} F_{j}^{IV}\bar{W}_{j}+ mg + m \nu^{2} \sum_{j=1}^{n} F_{j} \bar{W}_{j} -  \pi \rho_{\infty} b^{3} \omega^{2} \left ( l_{w} \sum_{j=1}^{n} F_{j} \dfrac{\bar{W}_{j}}{b} +  l_{\alpha} \sum_{j=1}^{m} H_{j} \bar{A}_{j} \right ) \right ] \cdot F_{i} dy = 0

\displaystyle \int_{0}^{s} \left [GJ \sum_{j=1}^{m} H_{j}^{II} \bar{A}_{j} - i_{yG} \nu^{2} \sum_{j=1}^{m} H_{j}\bar{A}_{j} + \pi \rho_{\infty} b^{4} \omega^{2} \left ( m_{w} \sum_{j=1}^{n} F_{j} \dfrac{\bar{W}_{j}}{b} + m_{\alpha} \sum_{j=1}^{m} H_{j} \bar{A}_{j} \right ) \right ] \cdot H_{i} dy = 0

Hemos obtenido un sistema de n ecuaciones (F_{1} \dots F_{n}) + m ecuaciones (H_{1} \dots H_{m}) para los n+m grados de libertad (W_{1} \dots W_{n}, A_{1} \dots A_{m}) . Matricialmente es un sistema de ecuaciones tal que:

 \left \{[A] \nu^{2} + [B]  \right \} \cdot \left \{ \bar{q} \right \} = [C]

Donde [C] representa la matriz de cargas externas (en este caso, el peso propio). Tendremos una solución particular + homogénea. El flameo se encuentra en la solución homogénea, por lo que podemos eliminar el peso propio para estudiar su comportamiento.

Por otra parte, se está integrando toda la placa de 0 a s. No obstante, un primer paso para adimensionalizar las ecuaciones es adimensionalizar la integración de forma que definimos una nueva variable de integración \psi = \dfrac{y}{s} tal que:

\dfrac{d^{n}}{dy^{n}} = \dfrac{1}{s^{n}} \cdot \dfrac{d^{n}}{d\psi ^{n}}

\displaystyle \int_{0}^{s} dy = s \cdot \int_{0}^{1} d\psi

Además, definimos los siguientes términos adimensionales:

\mu = \dfrac{m}{\pi \rho_{\infty} b^{2}}

r^{2} = \dfrac{i_{yG}}{m b^{2}}

\dfrac{GJ}{i_{yG} s^{2}} = \omega_{T}^{2}

\dfrac{EI}{ms^{4}} = \omega_{F}^{2} = \sigma^{2} \omega_{T}^{2}

\dfrac{b^{2} \omega_{T}^{2}}{U_{\infty}^{2}} = \dfrac{1}{V^{2}}

\dfrac{b^{2} \nu^{2}}{U_{\infty}^{2}} = p^{2}

\dfrac{b^{2} \omega^{2}}{U_{\infty}^{2}} = k^{2}

Adimensionalizando las ecuaciones empleando estas magnitudes adimensionales llegamos a las siguientes ecuaciones:

\displaystyle \int_{0}^{1} \left [\mu \cdot \sigma^{2} \cdot \omega_{T}^{2} \sum_{j=1}^{n} F_{j}^{IV} \dfrac{\bar{W}_{j}}{b} + \mu \cdot \nu^{2} \sum_{j=1}^{n} F_{j} \dfrac{\bar{W}_{j}}{b} - \omega^{2} l_{w} \sum_{j=1}^{n} F_{j} \dfrac{\bar{W}_{j}}{b} - \omega^{2} l_{\alpha} \sum_{j=1}^{m} H_{j} \bar{A}_{j} \right ] \cdot F_{i} d\psi = 0

\displaystyle \int_{0}^{1} \left [\mu \cdot r^{2} \cdot \omega_{T}^{2} \sum_{j=1}^{m} H_{j}^{II} \bar{A}_{j} - \mu \cdot r^{2} \cdot \nu^{2} \sum_{j=1}^{m} H_{j}\bar{A}_{j} + \omega^{2} m_{w} \sum_{j=1}^{n} F_{j} \dfrac{\bar{W}_{j}}{b} +\omega^{2} m_{\alpha} \sum_{j=1}^{m} H_{j} \bar{A}_{j} \right ] \cdot H_{i} d\psi = 0

Definiendo las integrales como:

\displaystyle I_{1}(i,j) = \int_{0}^{1} \left ( F_{j}^{IV} \cdot F_{i} \right ) d\psi

\displaystyle I_{2}(i,j) = \int_{0}^{1} \left ( F_{j} \cdot F_{i}  \right ) d\psi

\displaystyle I_{3}(i,j) = \int_{0}^{1} \left ( H_{j} \cdot F_{i} \right )  d\psi

\displaystyle I_{4}(i,j) = \int_{0}^{1}\left (  H_{j}^{II} \cdot H_{i} \right )  d\psi

\displaystyle I_{5}(i,j) = \int_{0}^{1}\left (  H_{j} \cdot H_{i} \right )  d\psi

\displaystyle I_{6}(i,j) = \int_{0}^{1} \left ( F_{j} \cdot H_{i}  \right ) d\psi

Multiplicando las ecuaciones por \dfrac{b^{2}}{U_{\infty}^{2}} y agrupando las integrales resulta el siguiente sistema matricial:

\left \{[MM] p^{2} + [KK] \dfrac{1}{V^{2}} + [AA] k^{2} \right \} \cdot \left \{ \bar{q} \right \} = 0

Con las matrices, de dimensiones {(n+m)x(n+m)} definidas como:

[MM] = \begin{bmatrix} \mu \cdot I_{2}(1,1) & \hdots & \mu \cdot I_{2}(1,n) & 0 & \hdots & 0\\ \vdots & \ddots & \vdots & \vdots & \ddots  & \vdots\\ \mu \cdot I_{2}(n,1) & \hdots & \mu \cdot I_{2}(n,n) & 0 & \hdots & 0\\ 0 & \hdots & 0 & -\mu \cdot r^{2} \cdot I_{5}(1,1) & \hdots & -\mu \cdot r^{2} \cdot I_{5}(1,m)\\ \vdots & \ddots & \vdots & \vdots & \ddots & \vdots\\ 0 & \hdots & 0 & -\mu \cdot r^{2} \cdot I_{5}(m,1) & \hdots & -\mu \cdot r^{2} \cdot I_{5}(m,m)\\ \end{bmatrix}

[KK] = \begin{bmatrix} \mu \cdot \sigma^{2} \cdot I_{1}(1,1) & \hdots & \mu \cdot \sigma^{2} \cdot I_{1}(1,n) & 0 & \hdots & 0\\ \vdots & \ddots & \vdots & \vdots &  \ddots & \vdots\\ \mu \cdot \sigma^{2} \cdot I_{1}(n,1) & \hdots & \mu  \cdot \sigma^{2} \cdot I_{1}(n,n) & 0 & \hdots & 0\\ 0 & \hdots & 0 & \mu \cdot r^{2} \cdot I_{4}(1,1) & \hdots & \mu \cdot r^{2} \cdot I_{4}(1,m)\\ \vdots & \ddots & \vdots & \vdots & \ddots & \vdots\\ 0 & \hdots & 0 & \mu \cdot r^{2} \cdot I_{4}(m,1) & \hdots & \mu \cdot r^{2} \cdot I_{4}(m,m)\\ \end{bmatrix}

[AA] = \begin{bmatrix} -l_{w} \cdot I_{2}(1,1) & \hdots & -l_{w} \cdot I_{2}(1,n) & -l_{\alpha} \cdot I_{3}(1,1) & \hdots & -l_{\alpha} \cdot I_{3}(1,m)\\ \vdots & \ddots & \vdots & \vdots & \ddots & \vdots\\ -l_{w} \cdot I_{2}(n,1) & \hdots &  -l_{w} \cdot I_{2}(n,n) & -l_{\alpha} \cdot I_{3}(n,1) & \hdots & -l_{\alpha} \cdot I_{3}(n,m)\\ m_{w} \cdot I_{6}(1,1) & \hdots & m_{w} \cdot I_{6}(1,n) & m_{\alpha} \cdot I_{5}(1,1) & \hdots & m_{\alpha} \cdot I_{5}(1,m)\\ \vdots & \ddots & \vdots & \vdots & \ddots & \vdots\\ m_{w} \cdot I_{6}(m,1) & \hdots & m_{w} \cdot I_{6}(m,n) & m_{\alpha} \cdot I_{5}(m,1) & \hdots & m_{\alpha} \cdot I_{5}(m,m)\\ \end{bmatrix}

Para encontrar una solución a la frecuencia compleja adimensional debemos hacer el determinante nulo, es decir:

\|  [MM] p^{2} + [KK] \dfrac{1}{V^{2}} + [AA] k^{2} \| = 0

Lo cual es un problema de autovalores generalizados de la forma:

\boxed{\| [A] - \lambda [B] \| = 0} con \sqrt{\lambda} = p, [A]= [KK] \dfrac{1}{V^{2}} + [AA] k^{2} y [B]= -[MM]

Es decir, estamos interesados en obtener los valores de la frecuencia compleja adimensional p = g + i \cdot k que hacen el determinante nulo. No obstante, el determinante también depende de la velocidad adimensional V y de la propia frecuencia reducida k. Para resolver el problema emplearemos el método p-k siguiendo el siguiente procedimiento:

  1. Se hace un barrido de la velocidad adimensional. Para cada V:
    1. Se hace un barrido de diferentes valores de k
    2. Para cada valor de k se resuelve el problema de autovalores generalizados y obtenemos \lambda \rightarrow p
    3. Sólo son válidas las soluciones en las que la parte imaginaria de p coinciden con el k barrido para que la solución sea coherente con la definición de fuerzas aerodinámicas.
  2. Para cada V hay asociadas tantas soluciones p como modos de deformación considerados, teniendo una parte compleja y una parte real.
  3. La primera velocidad dimensional U_{\infty} para la cual se tiene Re(p) > 0 es la velocidad de flameo.
  4. En dicho punto obtenemos la frecuencia reducida de flameo k_{F}, la velocidad de flameo U_{F} y la frecuencia de flameo \omega_{F}

Solución

Solución primer apartado

1 Función de aproximación para flexión/torsión

Solución de flameo para 1 (2 en total) función de aproximación en torsión/flexión
Solución de flameo para 1 (2 en total) función de aproximación en torsión/flexión

La velocidad de flameo es = 23.288289 m/s

La frecuencia de flameo es = 12.021141 rad/s

5 Funciones de aproximación para flexión/torsión

Solución de flameo para 5 (10 en total) funciones de aproximación en torsión/flexión
Solución de flameo para 5 (10 en total) funciones de aproximación en torsión/flexión

La velocidad de flameo es = 23.344114 m/s

La frecuencia de flameo es = 11.957095 rad/s

Solución segundo apartado

Para resolver este apartado se usa la lógica: con muelles de infinita rigidez el problema se resumiría al de un ala empotrada. En dicho caso, la solución sería la que se muestra a continuación:

Solución placa empotrada
Solución placa empotrada

Es decir, se obtendría una velocidad de flameo de 35.355709 m/s y una frecuencia de flameo de 14.746335 rad/s.

Dicha velocidad es un 51.82% superior con respecto a la velocidad en el problema inicial empleando 1 función de aproximación. Dada la similitud con lo que se propone en el último apartado y teniendo en cuenta que la solución es aproximada, se supone que el límite estaría alrededor del 50%.

A continuación se muestra el cálculo de la velocidad y frecuencia de flameo para diferentes valores de rigidez donde puede observarse claramente cómo se encuentra una asíntota al aumentar la rigidez de los muelles y el resultado tiende a lo obtenido con el ala empotrada:

Evolución de valores de flameo en función de la rigidez de los muelles
Evolución de valores de flameo en función de la rigidez de los muelles

Código de MATLAB

A continuación se presenta el código de MATLAB empleado para resolver el problema del flameo en base a las ecuaciones expuestas en este artículo para un número indefinido de funciones de aproximación (bastaría cambiar los valores de m y n):

%%% Metodo PK en ala continua

clear all
clc
close all

%% INICIALIZACIÓN
    % Datos
a=0;

% Características de la placa
b=0.5;
s=4;
t=0.01;

rho=2700;
E=70.5e9;

E=68.9e9;
G=26e9;
I=(2*b)*t^3/12;
J=1/3*2*b*t^3;

m=2*b*t*rho;
igy=(m*(2*b)^2)/12;

%Características del problema
rhoinf=1.225;
K=10000;

%Parámetros adimensionales
mu=m/(pi*rhoinf*b^2);
r=sqrt(igy/(m*b^2));
wf=sqrt(E*I/(m*s^4));
wt=sqrt(G*J/(igy*s^2));
sigma=wf/wt;





%% Funciones de aproximación
    %Funciones y sus derivadas con respecto a y, expresando phi=y/s
    nn=1; %n funciones aproximacion flexion
    nm=1; %m funciones aproximacion torsion

    
    [email protected](i) (3*E*I/(K*s^3))*(3+i)*(2+i)/6; 
    [email protected](i) ((3+i)*(2+i)*(K*s^3-3*E*I)-(6*K*s^3))/(6*K*s^3);
    a2=0;
    [email protected](i) -((3+i)*(2+i))/6;
    
            %Flexión
        F   = @(phi,i)   a0(i) + a1(i).*phi + a2.*phi.^2 + a3(i).*phi.^3 + phi.^(3+i);
        F_I = @(phi,i)   a1(i) + 2*a2.*phi + 3*a3(i).*phi.^2 + (3+i)*phi.^(2+i);    %Primera derivada
        F_II = @(phi,i)  2*a2 + 6*a3(i).*phi + (3+i)*(2+i)*phi.^(1+i);              %Segunda derivada
        F_III = @(phi,i)  6*a3(i) + (3+i)*(2+i)*(1+i)*phi.^(i);                     %Tercera derivada
        F_IV = @(phi,i)  (3+i)*(2+i)*(1+i)*i*phi.^(i-1);                            %Cuarta derivada
        
        
    [email protected](i) -(1+i);
    [email protected](i) -(G*J/(2*K*b^2*s))*(1+i);
    
             %Torsión
        H   = @(phi,i)   b0(i) + b1(i).*phi + phi.^(1+i);
        H_I = @(phi,i) b1(i) + (1+i)*phi;     %Primera derivada
        H_II = @(phi,i) (i+1).*i.*phi.^(i-1); %Segunda derivada
        
        
        %% Comprueba que las funciones de aproximación cumplen
        %Las condiciones de contorno del problema
        %Comprueba flexión
        residuof=0;
        for i=1:nm
            residuof=residuof+abs(F_II(0,i));
            residuof=residuof+abs(F_II(1,i));
            residuof=residuof+abs(s^3*F(0,i)+(E*I)/(2*K)*F_III(0,i));
            residuof=residuof+abs(F(1,i));
        end
        %residuof
        
        %Comprueba flexión
        residuot=0;
        for i=1:nm
            residuot=residuot+abs(s*H(0,i)-(G*J)/(2*K*b^2)*H_I(0,i));
            residuot=residuot+abs(H_I(1,i));
        end
        %residuot
        
        
    %% Calcula las integrales

for i=1:nn
    for j=1:nn
        I1(i,j)=quadgk(@(phi) F(phi,i).*F_IV(phi,j),0,1);
    end
end

    
for i=1:nn
    for j=1:nn
        I2(i,j)=quadgk(@(phi) F(phi,i).*F(phi,j),0,1);
    end
end

for i=1:nn
    for j=1:nm
        I3(i,j)=quadgk(@(phi) F(phi,i).*H(phi,j),0,1);
    end
end


for i=1:nm
    for j=1:nm
        I4(i,j)=quadgk(@(phi) H(phi,i).*H_II(phi,j),0,1);
    end
end


for i=1:nm
    for j=1:nm
        I5(i,j)=quadgk(@(phi) H(phi,i).*H(phi,j),0,1);
    end
end


for i=1:nm
    for j=1:nn
        I6(i,j)=quadgk(@(phi) F(phi,j).*H(phi,i),0,1);
    end
end
        
        

%% Cálculo de cortes
    %Matriz de MASA/INERCIA (no depende de k ni V)
    for i=1:nn
        for j=1:nn
            MM1(i,j)=mu*I2(i,j);
        end
    end
    
    for i=1:nn
        for j=1:nm
            MM2(i,j)=0;
        end
    end
    
    for i=1:nm
        for j=1:nn
            MM3(i,j)=0;
        end
    end
    
    for i=1:nm
        for j=1:nm
            MM4(i,j)=-r^2*mu*I5(i,j);
        end
    end
    
    MM=[MM1,MM2;MM3,MM4];
        
    
    % Genera valores para recorrer
    VV=linspace(0.01,5,200);  % Velocidades adimensionales del fluido
    kk=logspace(-4,3,500); % Rango de frecuencias reducidas (muchas cercanas al 0-1)
%     
    

    %Recorre cada velocidad adimensional
    for i=1:length(VV)
        V=VV(i);
        %Recorre cada frecuencia
        for j=1:length(kk)
            k=kk(j);
                
            %Cálculo aerodinámica
                C = besselk(1,1i*k)/(besselk(0,1i*k)+besselk(1,1i*k));        % Función de Theodorsen
                k2lw =  k^2 * (1-2.*1i.*C/k);                             % k^2 * lw en C
                k2la =  k^2 * (a+(1i/k).*(1+2.*(0.5-a).*C)+2.*C./k.^2); % k^2 * la en C
                k2mw =  k^2 * (a-(2.*1i.*(0.5+a).*C)./k);                    % k^2 * mw en C
                k2ma =  k^2 * (a^2+1/8-((0.5-a)*(1-2*(0.5+a)*C))*1i/k+(2*(0.5+a)*C)/k^2); % k^2 * ma en C
    
            %Matriz de rigidez
    for z=1:nn
        for x=1:nn
            KK1(z,x)=sigma^2*mu*I1(z,x)/V^2;
        end
    end
    
    for z=1:nn
        for x=1:nm
            KK2(z,x)=0;
        end
    end
    
    for z=1:nm
        for x=1:nn
            KK3(z,x)=0;
        end
    end
    
    for z=1:nm
        for x=1:nm
            KK4(z,x)=mu*r^2*I4(z,x)/V^2;
        end
    end
    
    KK=[KK1,KK2;KK3,KK4];
            
            %Matriz aerodinámica
    for z=1:nn
        for x=1:nn
            AA1(z,x)=-k2lw*I2(z,x);
        end
    end
    
    for z=1:nn
        for x=1:nm
            AA2(z,x)=-k2la*I3(z,x);
        end
    end
    
    for z=1:nm
        for x=1:nn
            AA3(z,x)=k2mw*I6(z,x);
        end
    end
    
    for z=1:nm
        for x=1:nm
            AA4(z,x)=k2ma*I5(z,x);
        end
    end
    
    AA=[AA1,AA2;AA3,AA4];
            
            
        % Determinante a resolver
        % | MM*p^2 + KK + AA | = 0
        % Lo cual se puede expresar como un problema de autovalores
        % | (KK + AA) -lambda*(-MM) | = 0 -> lambda=p^2
            A=KK+AA;
            B=-MM;
            
            lambdas=eig(A,B);
            p(j,:)=sqrt(lambdas);

                %Tomo la solución con parte imaginaria positiva
                for f=1:(nn+nm)
                    if imag(p(j,f))<0
                        p(j,f)=-p(j,f);
                    end
                end

           %Ordeno las soluciones
           
        ptemp = -1i*p(j,:); % Cambiar parte imaginaria por compleja para ordenar
        ptemp = sort(ptemp,'ComparisonMethod','real'); % Ordeno los valores según k (parte imag de p)
        p(j,:) = 1i*ptemp; %Vuelvo a generar el vector pero ordenado

        end
        
        %Intersección de k con las soluciones p
        %En psol defino las soluciones para cada velocidad

        for f=1:(nn+nm) %para cada solucion preliminar
        for j = 1:(length(kk) - 1)
        k   = kk(j);
        kmas1 = kk(j+1);
            if (imag(p(j,f)) - kk(j))*(imag(p(j+1,f)) - kk(j+1)) < 0 %Hay corte con k
                 psol(i,f) = p(j,f) - (p(j+1,f) - p(j,f))*(imag(p(j,f))-kk(j))/(imag(p(j+1,f)) - kk(j+1) - imag(p(j,f)) + kk(j));
                 break
            end
        end
        end
        
        
        
        

    end
    
%     Ya están todas las soluciones de p para cada velocidad adimensional V
%  %% Buscamos punto de flameo
    Vf=1e5; %Toma valor inicial alto
    
    
    for i=1:length(VV)-1
        V=VV(i);
        
        %Mira si en alguna solucion la parte real se anula
        %Sólo recorre la longitud de psol porque a veces no siempre
        %converge k
        for f=1:length(psol(i,:)) %para cada solucion válida
        if real(psol(i,f))*real(psol(i+1,f))<0 %Cambia el signo, Bolzano
            V_candidato(f)=VV(i) - real(psol(i,f))*(VV(i+1) - VV(i))/(real(psol(i+1,f))-real(psol(i,f)));
            k_candidato(f)=imag(psol(i,f)) - real(psol(i,f))*(imag(psol(i+1,f)) - imag(psol(i,f)))/(real(psol(i+1,f)) - real(psol(i,f)));
                  if V_candidato(f) < Vf
                    Vf = V_candidato(f);
                    kf = k_candidato(f);
                  end
        end
        end
           
    end
    
    
 %% Representación de solución

figure
subplot(2,1,1)
plot(VV,VV(:).*real(psol(:,1)),'DisplayName','Solución 1','LineWidth',2)
 hold on
 for f=2:length(psol(1,:))
      plot(VV,VV(:).*real(psol(:,f)),'DisplayName',sprintf('Solución %d',f),'LineWidth',2)
 end
title('Amortiguamiento','FontSize',15)
xlabel('$\displaystyle V = \frac{U}{b \omega_{T}}$','interpreter','latex','FontSize',15)
ylabel('$\displaystyle V*Re(p) = \frac{\gamma}{\omega_{T}}$','interpreter','latex','FontSize',15)
grid on
legend('Location','Best')
if Vf<10
plot(Vf,0,'ok','DisplayName','Flameo','Markersize',5,'MarkerFaceColor','k')
text(Vf,0.05,{['V_f = ',num2str(Vf)]},'FontSize', 12);
end
set(gca,'FontSize',14)
if Vf<10
ylim=get(gca,'ylim');
xlim=get(gca,'xlim');
text(xlim(1),ylim(2),{['k_f = ',num2str(kf)]},'FontSize', 12)
end

subplot(2,1,2)
plot(VV,VV(:).*imag(psol(:,1)),'DisplayName','Solución 1','LineWidth',2)
hold on
 for f=2:length(psol(1,:))
     plot(VV,VV(:).*imag(psol(:,f)),'DisplayName',sprintf('Solución %d',f),'LineWidth',2)
 end
title('Frecuencia','FontSize',15)
xlabel('$\displaystyle V = \frac{U}{b \omega_{T}}$','interpreter','latex','FontSize',15)
ylabel('$\displaystyle V*Im(p) = \frac{\omega}{\omega_{T}}$','interpreter','latex','FontSize',15)
grid on
legend('Location','Best')
%ylim([0 3])
if Vf<10
plot(Vf,Vf*kf,'ok','DisplayName','Flameo','Markersize',5,'MarkerFaceColor','k')
text(Vf,Vf*kf,{['V_f = ',num2str(Vf)];['\omega_f / \omega_T = ',num2str(Vf*kf)]},'FontSize', 12);
end
set(gca,'FontSize',14)

%set(gcf, 'Position',  [200, 200, 1000, 800])

fprintf('La velocidad de flameo es = %f m/s \n',Vf*b*wt);
fprintf('La frecuencia de flameo es = %f rad/s \n',Vf*kf*wt);

Deja un comentario