ANÁLISIS ESPECTRAL

Código Guía de Programación

Vibraciones Mecánicas

Considere el siguiente problema

En el estudio de las vibraciones en válvulas hidráulicas de control, la válvula y su vástago son modeladas como un sistema masa-resorte-amortiguador. En adición a la fuerza ejercida por el resorte y por el amortiguador, existe una fuerza ocasionada por la presión del fluido a la entrada y a la salida. Asumiendo que la rigidez es de 2500[N/m], la constante de amortiguamiento es 10[N-s/m], la masa de 250[gramos] y la tubería de 50[mm] de diámetro. Encontrar la respuesta en estado permanente de la válvula, cuando la presión varía de acuerdo a la siguiente gráfica.

Graficando la presión en la válvula.

Se supone que la presión se comporta como un diente de sierra, es decir aumenta gradualmente y se despresuriza súbitamente.

Pmax=50 (10^(3 ))(*[Pa]*);
τ0=2(*[s]*);
(*AQUÍ SE CAMBIA EL TIPO DE SEÑAL*)
Pt=SawtoothWave[{0,Pmax},(1/(τ0))t];
Plot[Pt,{t,0,20},PlotLabel->"Variación de la presión en la válvula",FrameLabel->{"t[s]","P(t)[Pa]"},Frame->True,GridLines->Automatic,PlotStyle->Orange,ExclusionsStyle->Dashed,ImageSize->Full]

Obteniendo la fuerza de excitación.

A partir de conocer el comportamiento de la presión, se procede por establecer la ecuación de la fuerza periódica.

Diametro=50 (10^(-3));
b0=Pmax(π((Diametro^2)));
(*AQUÍ SE CAMBIA DEL TIPO DE SEÑAL*)
bt=SawtoothWave[{0,b0},(1/(τ0))t];
Plot[bt,{t,0,20},PlotLabel->"Variación de la fuerza de presión en la válvula",FrameLabel->{"t[s]","b(t)[N]"},Frame->True,GridLines->Automatic,PlotStyle->Green,ExclusionsStyle->Dashed]

Datos del sistema mecánico.

A partir de lo establecido en el enunciado, se procede por definir los parámetros dinámicos del sistema, como las frecuencias y los factores de amortiguamiento.

m=0.25(*[kg]*);
k=2500(*[N/m]*);
c=10(*[N-s/m]*);
ω=(Sqrt[k/m]);
ω//N
ζ=(c/(2 ω m));
ζ//N
ωd=ω(Sqrt[1-(ζ^2)]);

Coeficientes de Fourier de la ecuación de Fuerza

Se toma la ecuación de la señal de entrada (o de la fuerza aplicada) y se obtienen sus coeficientes de Fourier

fe=(1/(τ0))
ωe=2π fe;
ω0=ωe;

a0=(2/(τ0))(∫_0^(τ0))bt\[DifferentialD]t

aj=(2/(τ0))(∫_0^(τ0))bt Cos[j ω0 t]\[DifferentialD]t

bj=(2/(τ0))(∫_0^(τ0))bt Sin[j ω0 t]\[DifferentialD]t

Expansión de la fuerza en series de Fourier

Una vez calculados los coeficientes, se sustituyen en la descripción en series de Fourier. Se definen cuántos armónicos se van a calcular (grado de discretización).

Elementos=300;
btFourier=(a0/2)+(∑_(j=1)^Elementos)aj Cos[j ω0 t]+(∑_(j=1)^Elementos)bj Sin[j ω0 t];

Análisis de la convergencia

Antes de proseguir, es importante cerciorarse de que la discretización es convergente para asegurar la precisión de la solución numérica.

GraficaFourier=Plot[{bt,btFourier},{t,0,6},PlotStyle->{{Orange,Dashed},{Green,Dotted}},PlotLegends->{"b(t)","b(t)Fourier"},Frame->True,PlotLabel->"Convergencia de la fuerza de excitación por series de Fourier",FrameLabel->{"t[s]","b(t)[s]"},GridLines->Automatic]

Magnificaciones y razones de frecuencia

Una vez validada la convergencia de la discretización, se procede por obtener las razones de frecuencia de cada combinación de armónicos y su respectiva magnificación.

Ωj=Table[((j ω0)/ω),{j,1,Elementos}];
HΩj=Table[(1/(Sqrt[(((1-((Ωj[[j]])^2)))^2)+(((2 ζ Ωj[[j]]))^2)])),{j,1,Elementos}];

Amplitudes de Fourier

Con base en las magnificaciones, y las deflexiones estáticas en términos de los coeficientes de Fourier, se procede a calcular las amplitudes magnificadas en esta descripción.

Raices=Table[(Sqrt[(aj^2)+(bj^2)]),{j,1,Elementos}];
cj=Table[HΩj[[j]] Raices[[j]],{j,1,Elementos}];
cjHΩj=Table[cj[[j]]HΩj[[j]],{j,1,Elementos}];

Espectro de frecuencias

Una vez calculadas las amplitudes, se grafican para contrastarlas con las deflexiones estáticas. Se recomienda poner en el fondo al factor de magnificación para observar la correlación quen existe entre los parámetros

Espectro=ListPlot[{cj,cjHΩj,HΩj},ImageSize->Large,AspectRatio->0.5,FillingStyle->{{Green,Dashed},{Orange,Dotted}},Filling->{1->Axis,2->Axis},PlotLegends->Placed[{Style["Entradas",FontFamily->"Montserrat Light"],Style["Salidas",FontFamily->"Montserrat Light"]},Above],PlotMarkers->None,Joined->{False,False,True},PlotStyle->{Green,Orange,Magenta},PlotRange->{{0,Elementos/2},{0,3}},Frame->True,FrameLabel->{{"(Subscript[c, j])/(Subscript[c, j])H((Subscript[Ω, j]))","H(Ω)"},{Style["Armónicos ((Subscript[jω, 0]))",FontFamily->"Montserrat Light"],Style["Espectro de frecuencias de la válvula a excitación periódica @ 0.5[Hz]",FontFamily->"Montserrat Light"]}}]

POSTPROCESO

La gráfica anterior es irrelevante si no se interpreta y se obtienen conclusiones de la misma. Quizá la cuestión aquí sería ¿Cómo interpertarla y tomar decisiones con base en ella?.

Para este caso en concreto...

Se observa que el sistema mecánico se comporta como un filtro pasa altas. Por lo que resultaría conveniente proponer una solución a bajas frecuencias si lo que se busca es disminuir la amplitud de la vibración actual.

1/Ωj[[1]]

Ωj[[1]]

Manipulate[
Solωe=Solve[Ωpropuesta==((ω0p)/ω),ω0p];
ωepropuesta=ω0p/.Solωe[[1]];
fepropuesta=((ωepropuesta)/(2π));
τpropuesta=(1/fepropuesta)
,{Ωpropuesta,0.1,100}]

Respuesta en el tiempo

c0=(a0/2);
(*tg inversa para Ψj, sólo está definida cuando j es impar, por lo tanto, todo debe adecuarse a números impares de j*)
Ψj=Table[-ArcTan[aj,bj]+(π/2),{j,1,Elementos,2}];
φΩjss=Table[-ArcTan[2 ζ Ωj[[j]],1-((Ωj[[j]])^2)],{j,1,Elementos,2}];
Cj=Table[cj[[j]],{j,1,Elementos,2}];
ΩΩj=Table[Ωj[[j]],{j,1,Elementos,2}];

tt=ω t;
ut=(1/k)(c0+(∑_(j=1)^(Elementos/2))Cj[[j]]Sin[ΩΩj[[j]]tt-φΩjss[[j]]+Ψj[[j]]]);

Graficaut=Plot[{(ut/(b0/k))(100)-50,(bt/b0)(100)},{t,0,6},PlotStyle->{{Orange},{Green}},PlotLegends->{"((u (t))/((Subscript[b, 0])/k))","((b (t))/(Subscript[b, 0]))"},Frame->True,PlotLabel->"Respuesta de la válvula a excitación periódica",FrameLabel->{"t[s]","Amplitud[%]"},GridLines->Automatic]