Vibraciones libres no amortiguadas en barras y chapas delgadas. Método de los elementos finitos

2. Barra
3. Chapa
4. Bibliografía consultada
5. Apéndice

1. 1.1. Teoría de vibraciones

   Cuando una estructura elástica está sujeta a una carga dinámica, el campo de desplazamiento varía con el tiempo debido a cuatro fuerzas: una es por la fuerza de inercia de la estructura, otra es por fricciones internas o externas, la tercera es por fuerzas elásticas actuando en el sistema, y la cuarta es por cualquier fuerza externa aplicada. Como el sistema está en equilibrio estático, la fuerza externa aplicada debe ser igual a la suma de todas las demás fuerzas. Expresado matemáticamente,

   siendo M la masa del sistema, C el coeficiente de amortiguamiento, K la constante elástica, F la fuerza externa, x el desplazamiento y t el tiempo. Es la llamada ecuación de movimiento.

   Si se toma el caso ideal en el cual el sistema no tiene ni amortiguamiento ni fuerzas externas, la ecuación de movimiento se reduce a

   Esta ecuación expresa la condición de natural vibración, en el cual en cualquier instante las fuerzas que hacen restaurar al sistema son las de inercia y las elásticas. Con esto se deduce que para barras y chapas delgadas, las vibraciones sólo dependen del módulo de Young y de la densidad del material. Los estados de vibración natural se llaman modos naturales, y las frecuencias de vibración son las frecuencias naturales.

   Para hallar dichos modos y frecuencias naturales se asume que el desplazamiento puede ser expresado como

   donde x’ es la parte del desplazamiento nodal que se asume independiente del tiempo, i es el número imaginario, y w es la frecuencia natural.

   Derivando la ecuación anterior dos veces respecto del tiempo, se obtiene

   que substituyendo en la ecuación de movimiento reducida, se llega a

   

   Reordenando, se obtiene

   

   Como eiwt no es cero,

   

   La ecuación anterior es un arreglo de ecuaciones lineales homogéneas, que tiene una solución no trivial sólo si el determinante de la matriz de coeficientes de x’ es cero. Es decir,

   

   que es un problema de autovectores y autovalores, cumpliendo la siguiente ecuación:

   

   con vi el autovector asociado al autovalor l i, que es igual a wi2.

   2. Un objeto de gran importancia en el análisis de elementos finitos es la elección de apropiadas funciones de interpolación dentro de cada elemento. Dichas funciones de interpolación N se pueden expresar de la siguiente manera:

      

      siendo f i el valor de la función a aproximar en el nodo local i.

      Las funciones de interpolación no se pueden elegir arbitrariamente, requieren continuidad en la función, y en algunos casos, en la derivada, entre los elementos. Se elige que sean polinomios, debido a su facilidad para integrar y derivar. Además, tiene que cumplir que la función Ni debe valer 1 en el nodo local i. Es decir,

      

      Más adelante se verá cómo quedan las funciones de interpolación en una dimensión lineal y cúbica, y en dos dimensiones.

   3. Funciones de interpolación
   4. Matrices locales

   La matriz de masa local se expresa de la siguiente manera:

   

   siendo N las funciones de interpolación y r la densidad, integrado en el volumen del elemento.

   Si tomamos la densidad constante en el caso unidimensional, la fórmula anterior queda

   

   con a el largo del elemento.

   Análogamente en el caso bidimensional,

   

   siendo a el largo del elemento en la dirección x y b el ancho en la dirección y.

   La matriz elástica local se expresa de la siguiente manera:

   

   con E el módulo de Young, integrado en el volumen del elemento.

   Si tomamos el módulo de Young constante, en el caso unidimensional, la fórmula anterior queda

   

   con a el largo del elemento.

   Análogamente en el caso bidimensional,

   

   siendo a el largo del elemento en la dirección x y b el ancho en la dirección y.

2. Introducción

   2.1. En voladizo

   Para poder realizar los cálculos por elementos finitos, la barra se dividió en segmentos iguales o elementos de longitud L, como se observa en la Figura 1. Se observa que está en voladizo, pues contiene una condición de borde en x = 0, en el cual no hay desplazamiento nodal, siempre tiene valor cero (y = 0). Los números recuadrados representan el número de elemento, los números en negrita sobre los nodos la numeración global de dichos nodos, y los números inclinados por debajo de los nodos la numeración local.

   

   Figura 1. Representación en elementos de una barra en voladizo

   El paso siguiente es encontrar las funciones de interpolación (Figuras 2 y 3), en este caso lineales, con todos los requisitos que se indicaron en la sección 1.2., que resultaron ser las siguientes:

   

   

   Figura 2. Función 1 de interpolación lineal para la barra

   

   

   Figura 3. Función 2 de interpolación lineal para la barra

   Una vez que ya se obtuvieron las funciones de interpolación, se pueden hallar las matrices locales como se explicó en la sección 1.3., quedando de la siguiente manera:

   

   Ahora se deben obtener las matrices globales de masa y elástica. Como las funciones de interpolación deben ser continuas y de valor 1 entre los elementos, y observando la forma que tienen en las Figuras 2 y 3, se llega a la conclusión que la función de interpolación 2 del elemento anterior debe empalmar con la función de interpolación 1 del elemento siguiente. Es decir, en la matriz global se deben sumar las matrices locales, pero superponiéndose la segunda fila y columna del elemento anterior con la primera fila y columna del elemento siguiente.

   Gráficamente, para la matriz de masa,

   

   Los cuadrados dentro de la matriz representan los elementos finitos, que al superponerse se suman los elementos de la matriz local. Los elementos de la matriz global que no tienen valor son cero.

   Análogamente, para la matriz elástica,

   

   Como condición de borde se necesita que el desplazamiento nodal sea cero en el nodo 1 del elemento 1. Para que ello suceda se pide que la función de interpolación 1 del elemento 1 valga cero, así la función desplazamiento siempre vale cero en el nodo 1 global, pues la función de interpolación 2 ya vale cero. Entonces lo que se hace es colocar un 1 en la diagonal de la fila 1 y columna 1 de las matrices globales, y todo lo restante de dicha fila y columna cero. Las matrices globales quedan de la siguiente manera:

   

   

   Al ya haber obtenido las matrices globales M y K, se está en condiciones de resolver la ecuación de autovectores y autovalores .

   Tomando una barra delgada de acero (E = 210.109 N/m2, r = 7800 kg/m3) de 1 m de largo con 500 elementos, los resultados fueron los siguientes:

   

   Figura 4. Desplazamientos nodales y frecuencias naturales para la barra en voladizo

   Se observan los desplazamientos nodales para cada frecuencia normal de vibración a lo largo de toda la barra. En x = 0, por condición de borde, siempre hay un nodo, mientras que en el extremo derecho siempre hay un antinodo. Sólo se graficaron los primeros cinco autovectores, con sus frecuencias naturales asociadas, pues hay tantos autovectores y autovalores como grados de libertad, que en este caso es la cantidad de elementos mas uno (cantidad de nodos).

   2.2. 2 vínculos

   En este caso, la barra se dividió de la misma manera que la anterior, pero esta vez contiene dos condiciones de borde: en x = 0 y en x = nL (n es el número de elementos) no hay desplazamiento nodal, y = 0. Como el caso anterior, los números recuadrados representan el número de elemento, los números en negrita sobre los nodos la numeración global de dichos nodos, y los números inclinados por debajo de los nodos la numeración local.

   

   Figura 5. Representación en elementos de una barra con dos vínculos

   Se realizaron dos funciones de interpolación distintas: lineal y cúbica. Como las funciones de interpolación lineal son las mismas que en el caso anterior, sólo se mostrarán las funciones de interpolación cúbicas (Figuras 6, 7, 8 y 9), que cumplen con todos los requisitos que se indicaron en la sección 1.2., que son las siguientes:

   

   

   Figura 6. Función 1 de interpolación cúbica para la barra

   

   

   Figura 7. Función 2 de interpolación cúbica para la barra

   

   

   Figura 8. Función 3 de interpolación cúbica para la barra

   

   

   Figura 9. Función 4 de interpolación cúbica para la barra

   Una vez que ya se obtuvieron las funciones de interpolación, se pueden hallar las matrices locales como se explicó en la sección 1.3., quedando de la siguiente manera:

   

   Ahora se deben obtener las matrices globales de masa y elástica. Como las funciones de interpolación deben ser continuas en función y en derivada entre los elementos, y observando la forma que tienen en las Figuras 6, 7, 8 y 9, se llega a la conclusión que las funciones de interpolación 3 y 4 del elemento anterior deben empalmar con las funciones de interpolación 1 y 2 del elemento siguiente, respectivamente. Es decir, en la matriz global se deben sumar las matrices locales, pero superponiéndose la tercera y cuarta fila y columna del elemento anterior con la primera y segunda fila y columna del elemento siguiente, respectivamente.

   Gráficamente, para la matriz de masa,

   

   Los cuadrados dentro de la matriz representan los elementos, que al superponerse se suman los elementos de la matriz local. Los elementos de la matriz global que no tienen valor son cero.

   Análogamente, para la matriz elástica,

   

   Como condición de borde se necesita que el desplazamiento nodal sea cero en el nodo 1 del elemento 1 y en el nodo 2 del último elemento. Para que ello suceda se pide que la función de interpolación 1 del elemento 1 y la función de interpolación 3 del último elemento valgan cero, así la función desplazamiento siempre vale cero en el primer y último nodo global, pues las otras funciones de interpolación ya valen cero. Entonces lo que se hace es colocar un 1 en la diagonal de la primera y penúltima fila y columna de las matrices globales, y todo lo restante de dichas filas y columnas cero. Las matrices globales quedan de la siguiente manera:

   

   

   Al ya haber obtenido las matrices globales M y K, se está en condiciones de resolver la ecuación de autovectores y autovalores .

   Tomando una barra delgada de acero (E = 210.109 N/m2, r = 7800 kg/m3) de 1 m de largo con 500 elementos, los resultados fueron los siguientes:

   

   Figura 10. Desplazamientos nodales y frecuencias naturales para la barra en voladizo

   Se observan los desplazamientos nodales para cada frecuencia normal de vibración a lo largo de toda la barra. En x = 0 y en x = nL (n es el número de elementos), por condición de borde, siempre hay un nodo. Sólo se graficaron los primeros cinco autovectores, con sus frecuencias naturales asociadas.

   Como se dijo anteriormente, además se resolvió el problema con funciones de interpolación lineales, en el cual es igual a la barra en voladizo con la diferencia en las condiciones de borde, en el cual se necesita que el desplazamiento nodal sea cero en el último nodo global. Para que ello suceda se pide que, además de que la función de interpolación 1 del elemento 1 valga cero, que la función de interpolación 2 del último elemento también valga cero, así la función desplazamiento siempre vale cero en el primer y último nodo global, pues las otras funciones de interpolación ya valen cero. Entonces la única diferencia con el caso de la barra en voladizo lineal es que se debe colocar un 1 en la diagonal de la última fila y columna de las matrices globales, y todo lo restante de dicha filas y columna cero.

   La gran ventaja de usar funciones de interpolación cúbicas respecto de las lineales es la rápida convergencia, es decir, para una misma cantidad de elementos, se obtiene un error porcentual mucho menor. Esto se puede observar en la Figura 11.

   

   Figura 11. Error porcentual en función del número de elementos, para distintas frecuencias naturales, con funciones de interpolación lineales y cúbicas

   El error porcentual se calculó comparando con la expresión exacta de frecuencias naturales

   

   con fn la n frecuencia natural de vibración en Hz y l el largo de la barra. n sólo toma números naturales: 1, 2, 3 …

   En el Apéndice se muestran para todos los casos el código fuente en el lenguaje Fortran.

3. Barra

   3.1. En voladizo

   Como se tomó una chapa cuadrada, se dividió en áreas iguales o elementos cuadrados de ancho y longitud L, como se observa en la Figura 12. Como está en voladizo, contiene una condición de borde en x = 0, en el cual no hay desplazamiento nodal, siempre tiene el valor cero (z = 0). Como en los casos anteriores, los números recuadrados representan el número de elemento, los números en negrita sobre los nodos la numeración global de dichos nodos, y los números inclinados la numeración local.

   Ahora nuevamente como antes se deben encontrar las funciones de interpolación, que en este caso no son curvas sino superficies, pues se está trabajando en dos dimensiones. Se debe cumplir con todos los requisitos que se indicaron en la sección 1.2., resultando las funciones de interpolación que se observan en las Figuras 13, 14, 15 y 16.

   

   Figura 12. Representación en elementos de una chapa cuadrada en voladizo

   

   

   Figura 13. Función 1 de interpolación para la chapa

   

   

   Figura 14. Función 2 de interpolación para la chapa

   

   

   Figura 15. Función 3 de interpolación para la chapa

   

   

   Figura 16. Función 4 de interpolación para la chapa

   Una vez halladas las funciones de interpolación, se pueden encontrar las matrices locales como se explicó en la sección 1.3., quedando de la siguiente manera:

   

   Ahora se deben obtener las matrices globales de masa y elástica. Como las funciones de interpolación deben ser continuas y de valor 1 entre los elementos, se debe llegar a una coherencia entre los nodos locales y globales. Es decir, y observando la Figura 12, el llenado de la matriz global debe ser teniendo en cuenta la relación entre nodos locales y globales. Tomando el elemento 1 de la Figura 12, los nodos locales 1, 2, 3 y 4 se corresponden con los nodos globales 1, 2, 5 y 6 respectivamente, entonces la fila y columna 1, 2, 3 y 4 de la matriz local se coloca en la fila y columna 1, 2, 5 y 6 de la matriz global. Si luego se toma el elemento 2, los nodos locales 1, 2, 3 y 4 se corresponden con los nodos globales 2, 3, 6 y 7 respectivamente, entonces la fila y columna 1, 2, 3 y 4 de la matriz local va en la fila y columna 2, 3, 6 y 7 de la matriz global, y así sucesivamente, quedando de la siguiente manera:

   

    

   

   Como condición de borde se necesita que el desplazamiento nodal sea cero en x = 0, es decir, en este caso de la Figura 12, los nodos globales 1, 5, 9 y 13. Para que ello suceda se pide que las funciones de interpolación que valen 1 en dichos nodos valgan cero, así la función desplazamiento siempre vale cero en aquellos nodos globales, porque las demás funciones de interpolación ya valen cero. Entonces se coloca un 1 en la diagonal de las filas y columnas 1, 5, 9 y 13 de las matrices globales, y todo lo restante de dichas filas y columnas cero.

   Al obtener las matrices globales M y K, se resuelve la ecuación de autovectores y autovalores .

   Tomando una chapa delgada de acero (E = 210.109 N/m2, r = 7800 kg/m3) de 1 m de lado con 10 por 10 elementos, los resultados fueron los siguientes:

   

   Figura 17. Desplazamientos nodales para la chapa en voladizo para la primera frecuencia natural de 1059 Hz

    

   

   Figura 18. Desplazamientos nodales para la chapa en voladizo para la segunda frecuencia natural de 3178 Hz

   

   Figura 19. Desplazamientos nodales para la chapa en voladizo para la tercera frecuencia natural de 5296 Hz

   

   Figura 20. Desplazamientos nodales para la chapa en voladizo para la cuarta frecuencia natural de 5902 Hz

   Se observan los desplazamientos nodales para las primeras cuatro frecuencias normales de vibración en toda la superficie de la chapa. En x = 0, por condición de borde, siempre hay un nodo.

   3.2. 4 vínculos

   El análisis para la chapa cuadrada delgada con 4 vínculos (Figura 21) es igual al caso anterior en voladizo, sólo cambian las condiciones de borde, en el cual en x = 0, x = nL, y = 0, e y = nL (con n el número de elementos) no hay desplazamientos nodales (z = 0). Entonces, tomando como ejemplo la Figura 21, en los nodos globales 1, 2, 3, 4, 5, 8, 9, 12, 13, 14, 15 y 16 las funciones de interpolación que valen 1 en dichos nodos tienen que valer cero, así la función desplazamiento siempre vale cero en dichos nodos globales, porque las demás funciones de interpolación ya valen cero. Entonces se coloca un 1 en la diagonal de las filas y columnas 1, 2, 3, 4, 5, 8, 9, 12, 13, 14, 15 y 16 de las matrices globales, y todo lo restante de dichas filas y columnas cero.

   

   Figura 21. Representación en elementos de una chapa cuadrada con 4 vínculos

   Resolviendo nuevamente la ecuación de autovectores y autovalores, y tomando los parámetros anteriores, los resultados fueron los siguientes:

   

   Figura 22. Desplazamientos nodales para la chapa con 4 vínculos para la primera frecuencia natural de 3350 Hz

   

   Figura 23. Desplazamientos nodales para la chapa con 4 vínculos para la segunda frecuencia natural de 5607 Hz (degenerada)

   

   Figura 24. Desplazamientos nodales para la chapa con 4 vínculos para la segunda frecuencia natural de 5607 Hz (degenerada)

   

   Figura 25. Desplazamientos nodales para la chapa con 4 vínculos para la tercera frecuencia natural de 6875 Hz

   Se puede ver que hay una degeneración en la segunda frecuencia natural de vibración, es decir, que para dos desplazamientos nodales diferentes, tienen la misma frecuencia natural. Esto se debe a la simetría que contiene el sistema.

   3.3. Convergencia

   En las Figuras 26, 27, 28 y 29 se analizaron las convergencias de las frecuencias naturales según la cantidad de elementos utilizados, para la chapa en voladizo y con 4 vínculos, para la primera y cuarta frecuencia natural de vibración. Notar que la escala del número de elementos es logarítmica, indica que rápidamente converge hacia el número exacto.

   

   Figura 26. 1º frecuencia natural en función del número de elementos para chapa en voladizo

   

   Figura 27. 4º frecuencia natural en función del número de elementos para chapa en voladizo

   

   Figura 28. 1º frecuencia natural en función del número de elementos para chapa con 4 vínculos

   

   Figura 29. 4º frecuencia natural en función del número de elementos para chapa con 4 vínculos

   A modo informativo, se hallaron los errores porcentuales para la cantidad de mil elementos, valor límite en el programa ABAQUS versión estudiante.

   3.4. Cambio de densidad

   Por último, a la chapa con 4 vínculos, en la parte superior derecha se cambió el material. En vez de usar acero se usó aluminio, material con menor densidad. Para ser tomado en los cálculos, lo que se hizo fue multiplicar por 0,346 (que es 2700 / 7800, la densidad del aluminio dividido la densidad del acero) las matrices locales de masa que pertenecen a los elementos que son de aluminio. Los resultados fueron los siguientes:

   

   Figura 30. Curvas de nivel para la chapa de acero con 4 vínculos con una frecuencia de 5607 Hz (degenerada)

   

   Figura 31. Curvas de nivel para la chapa de acero con 4 vínculos con una frecuencia de 5607 Hz (degenerada)

   

   Figura 32. Curvas de nivel para la chapa con aluminio con 4 vínculos a una frecuencia de 5622 Hz

    

   Figura 33. Curvas de nivel para la chapa con aluminio con 4 vínculos a una frecuencia de 4990 Hz

   En las figuras anteriores se puede observar que se rompió la degeneración por cambiar el material en una parte que no es en el centro de la chapa. Esto sucede porque la chapa dejó de ser simétrica, y para que haya degeneración tiene que haber simetría.

   En el Apéndice se muestran para todos los casos de la chapa el código fuente en el lenguaje Fortran.

4. Chapa
5. Bibliografía consultada

• K. H. Heubner, D. L. Dewhirst, D. E. Smith, T. G. Byrom, The Finite Element Method for Engineers, Fourth Edition, John Wiley and Sons, 2001
• D. L. Logan, A First Course in the Finite Element Method, Third Edition, Wadsworth Group, Brooks/Cole, 2002
• G. R. Buchanan, Finite Element Analysis, Schaum’s Outline Series, McGraw-Hill, 1995

2. Apéndice

5.1. Código fuente para la barra en voladizo

program EF

use imsl

implicit none

real(8)::deltax,l

real(8),allocatable::K(:,:),M(:,:),eval(:),evec(:,:)

integer(4)::elem,i

write(*,*)"Ingrese cantidad de elementos y el largo"

read(*,*)elem,l

allocate(K(elem+1,elem+1),M(elem+1,elem+1),eval(elem+1),evec(elem+1,elem+1))

deltax=l/elem

K(elem+1,elem+1)=1.

K(1,1)=1.

M(elem+1,elem+1)=2.

M(1,1)=1.

do i=2,elem

K(i,i)=2.

K(i,i+1)=-1.

K(i+1,i)=-1.

M(i,i)=4.

M(i,i+1)=1.

M(i+1,i)=1.

enddo

K=K*6.*210.E9/(7800.*deltax**2)

call dgvcsp(elem+1,K,elem+1,M,elem+1,EVAL,EVEC,elem+1)

write(*,*)"Modos de vibracion"

write(*,'(F10.3)')dsqrt(eval)/(2.*3.141592654)

pause

write(*,*)"Desplazamientos nodales"

call dwrrrl('',elem+1,elem+1,evec,elem+1,0,'(F6.3)','','')

end program EF

5.2. Código fuente para la barra con dos vínculos lineal

program EF

use imsl

implicit none

real(8)::deltax,l,E,d,j

real(8),allocatable::K(:,:),M(:,:),eval(:),evec(:,:),frec(:)

integer(4)::elem,i

write(*,*)"Ingrese cantidad de elementos y el largo"

read(*,*)elem,l

E=210.E9

d=7800.

allocate(K(elem+1,elem+1),M(elem+1,elem+1),eval(elem+1),evec(elem+1,elem+1),frec(elem+1))

deltax=l/elem

K(elem+1,elem+1)=1.

K(1,1)=1.

M(elem+1,elem+1)=1.

M(1,1)=1.

do i=2,elem

K(i,i)=2.

K(i,i+1)=-1.

K(i+1,i)=-1.

M(i,i)=4.

M(i,i+1)=1.

M(i+1,i)=1.

enddo

K(1,2)=0.

K(2,1)=0.

K(elem,elem+1)=0.

K(elem+1,elem)=0.

M(1,2)=0.

M(2,1)=0.

M(elem,elem+1)=0.

M(elem+1,elem)=0.

K=K*6.*E/(d*deltax**2)

call dgvcsp(elem+1,K,elem+1,M,elem+1,EVAL,EVEC,elem+1)

write(*,*)"Modos de vibracion"

frec=dsqrt(eval)/(2.*3.141592654)

write(*,'(F8.2)')frec

pause

write(*,*)"Error porcentual"

do j=1.,4.

write(*,'(F12.5)')j*dsqrt(E/d)/(2.*l)

write(*,'(F12.5)')(dabs(frec(elem+2-j)-j*dsqrt(E/d)/(2.*l))/(j*dsqrt(E/d)/(2.*l)))*100.

write(*,*)

enddo

pause

write(*,*)"Desplazamientos nodales"

call dwrrrl('',elem+1,elem+1,evec,elem+1,0,'(F6.3)','','')

end program EF

5.3. Código fuente para la barra con dos vínculos cúbica

program EF

use imsl

implicit none

real(8)::deltax,l,E,d,j

real(8),allocatable::K(:,:),M(:,:),eval(:),frec(:),evec(:,:)

integer(4)::elem,i

write(*,*)"Ingrese cantidad de elementos y el largo"

read(*,*)elem,l

E=210.E9

d=7800.

allocate(K(2*elem+2,2*elem+2),M(2*elem+2,2*elem+2),eval(2*elem+2),frec(2*elem+2),evec(2*elem+2,2*elem+2))

deltax=l/elem

do i=1,elem

M(2*i,2*i+1)=13.*deltax

M(2*i+1,2*i)=13.*deltax

M(2*i,2*i+2)=-3.*deltax**2

M(2*i+2,2*i)=-3.*deltax**2

K(2*i,2*i+1)=-3.*deltax

K(2*i+1,2*i)=-3.*deltax

K(2*i,2*i+2)=-(deltax**2)

K(2*i+2,2*i)=-(deltax**2)

enddo

do i=2,elem

M(2*i-1,2*i-1)=312.

M(2*i,2*i)=8.*deltax**2

M(2*i-1,2*i+1)=54.

M(2*i+1,2*i-1)=54.

M(2*i-1,2*i+2)=-13.*deltax

M(2*i+2,2*i-1)=-13.*deltax

K(2*i-1,2*i-1)=72.

K(2*i,2*i)=8.*deltax**2

K(2*i-1,2*i+1)=-36.

K(2*i+1,2*i-1)=-36.

K(2*i-1,2*i+2)=3.*deltax

K(2*i+2,2*i-1)=3.*deltax

enddo

M(1,1)=1.

M(2*elem+1,2*elem+1)=1.

M(2,2)=4.*deltax**2

M(2*elem+2,2*elem+2)=4.*deltax**2

M(2*elem+1,2*elem+2)=-22.*deltax

M(2*elem+2,2*elem+1)=-22.*deltax

K(1,1)=1.

K(2*elem+1,2*elem+1)=1.

K(2,2)=4.*deltax**2

K(2*elem+2,2*elem+2)=4.*deltax**2

K(2*elem+1,2*elem+2)=-3.*deltax

K(2*elem+2,2*elem+1)=-3.*deltax

do i=2,2*elem+1

M(1,i)=0.

M(i,1)=0.

M(i-1,2*elem+1)=0.

M(2*elem+1,i-1)=0.

M(2*elem+1,2*elem+2)=0.

M(2*elem+2,2*elem+1)=0.

K(1,i)=0.

K(i,1)=0.

K(i-1,2*elem+1)=0.

K(2*elem+1,i-1)=0.

K(2*elem+1,2*elem+2)=0.

K(2*elem+2,2*elem+1)=0.

enddo

K=K*14.*E/(d*deltax**2)

call dgvcsp(2*elem+2,K,2*elem+2,M,2*elem+2,EVAL,EVEC,2*elem+2)

write(*,*)"Modos de vibracion"

frec=dsqrt(eval)/(2.*3.141592654)

write(*,'(F8.2)')frec

pause

write(*,*)"Error porcentual"

do j=1.,4.

write(*,'(F12.5)')j*dsqrt(E/d)/(2.*l)

write(*,*)(dabs(frec(2*elem+3-j)-j*dsqrt(E/d)/(2.*l))/(j*dsqrt(E/d)/(2.*l)))*100.

write(*,*)

enddo

pause

write(*,*)"Desplazamientos nodales"

call dwrrrl('',2*elem+2,2*elem+2,evec,2*elem+2,0,'(F6.3)','','')

end program EF

5.4. Código fuente para la chapa en voladizo

program EF

use imsl

implicit none

real(8)::delta,l

real(8),dimension(4,4)::a,b

real(8),allocatable::K(:,:),M(:,:),evec(:,:),eval(:)

integer(4)::elem,i,nex,ney,j,p,q

write(*,*)"Ingrese cantidad de elementos y el lado"

read(*,*)elem,l

delta=l/elem

allocate(K((elem+1)**2,(elem+1)**2),M((elem+1)**2,(elem+1)**2),eval((elem+1)**2),evec((elem+1)**2,(elem+1)**2))

do i=1,4

a(i,i)=4.

b(i,i)=4.

enddo

do i=1,3

a(i,i+1)=-1.

a(i+1,i)=-1.

b(i,i+1)=2.

b(i+1,i)=2.

enddo

do i=1,2

a(i,i+2)=-2.

a(i+2,i)=-2.

b(i,i+2)=1.

b(i+2,i)=1.

enddo

a(1,4)=-1.

a(4,1)=-1.

b(1,4)=2.

b(4,1)=2.

do nex=1,elem

do ney=1,elem

do i=1,4

do j=1,4

if (i.eq.1) then

p=nex+(ney-1)*(elem+1)

elseif (i.eq.2) then

p=nex+1+(ney-1)*(elem+1)

elseif (i.eq.3) then

p=nex+(elem+1)+(ney-1)*(elem+1)

elseif (i.eq.4) then

p=nex+(elem+2)+(ney-1)*(elem+1)

endif

if (j.eq.1) then

q=nex+(ney-1)*(elem+1)

elseif (j.eq.2) then

q=nex+1+(ney-1)*(elem+1)

elseif (j.eq.3) then

q=nex+(elem+1)+(ney-1)*(elem+1)

elseif (j.eq.4) then

q=nex+(elem+2)+(ney-1)*(elem+1)

endif

K(p,q)=K(p,q)+a(i,j)

M(p,q)=M(p,q)+b(i,j)

enddo

enddo

enddo

enddo

do i=1,(elem+1)**2,elem+1

do j=1,(elem+1)**2

K(i,j)=0.

K(j,i)=0.

M(i,j)=0.

M(j,i)=0.

enddo

K(i,i)=1.

M(i,i)=1.

enddo

K=K*(210.E9)*6./(7800.*delta**2)

call dgvcsp((elem+1)**2,K,(elem+1)**2,M,(elem+1)**2,EVAL,EVEC,(elem+1)**2)

write(*,*)"Modos de vibracion"

write(*,'(F15.3)')dsqrt(eval)/(2.*3.141592654)

pause

write(*,*)"Error porcentual"

write(*,*)"1 modo",(dabs(dsqrt(eval((elem+1)**2))/(3.141592654*2.)-1059.)/1059.)*100.

write(*,*)"2 modo",(dabs(dsqrt(eval((elem+1)**2-2))/(3.141592654*2.)-3178.)/3178.)*100.

write(*,*)"3 modo",(dabs(dsqrt(eval((elem+1)**2-4))/(3.141592654*2.)-5296.)/5296.)*100.

write(*,*)"4 modo",(dabs(dsqrt(eval((elem+1)**2-6))/(3.141592654*2.)-5902.)/5902.)*100.

pause

write(*,*)"Desplazamientos nodales"

call dwrrrl('',(elem+1)**2,(elem+1)**2,evec,(elem+1)**2,0,'(F6.3)','','')

open(1,file="placa.dat")

write(1,'(F12.5)')evec(:,118)

end program EF

5.5. Código fuente para la chapa con 4 vínculos

program EF

use imsl

implicit none

real(8)::delta,l

real(8),dimension(4,4)::a,b

real(8),allocatable::K(:,:),M(:,:),evec(:,:),eval(:)

integer(4)::elem,i,nex,ney,j,p,q

write(*,*)"Ingrese cantidad de elementos y el lado"

read(*,*)elem,l

delta=l/elem

allocate(K((elem+1)**2,(elem+1)**2),M((elem+1)**2,(elem+1)**2),eval((elem+1)**2),evec((elem+1)**2,(elem+1)**2))

do i=1,4

a(i,i)=4.

b(i,i)=4.

enddo

do i=1,3

a(i,i+1)=-1.

a(i+1,i)=-1.

b(i,i+1)=2.

b(i+1,i)=2.

enddo

do i=1,2

a(i,i+2)=-2.

a(i+2,i)=-2.

b(i,i+2)=1.

b(i+2,i)=1.

enddo

a(1,4)=-1.

a(4,1)=-1.

b(1,4)=2.

b(4,1)=2.

do nex=1,elem

do ney=1,elem

do i=1,4

do j=1,4

if (i.eq.1) then

p=nex+(ney-1)*(elem+1)

elseif (i.eq.2) then

p=nex+1+(ney-1)*(elem+1)

elseif (i.eq.3) then

p=nex+(elem+1)+(ney-1)*(elem+1)

elseif (i.eq.4) then

p=nex+(elem+2)+(ney-1)*(elem+1)

endif

if (j.eq.1) then

q=nex+(ney-1)*(elem+1)

elseif (j.eq.2) then

q=nex+1+(ney-1)*(elem+1)

elseif (j.eq.3) then

q=nex+(elem+1)+(ney-1)*(elem+1)

elseif (j.eq.4) then

q=nex+(elem+2)+(ney-1)*(elem+1)

endif

K(p,q)=K(p,q)+a(i,j)

M(p,q)=M(p,q)+b(i,j)

enddo

enddo

enddo

enddo

do i=1,elem+2

do j=1,(elem+1)**2

K(i,j)=0.

K(j,i)=0.

M(i,j)=0.

M(j,i)=0.

enddo

K(i,i)=1.

M(i,i)=1.

enddo

do i=elem+2,elem**2+elem,elem+1

do j=1,(elem+1)**2

K(i,j)=0.

K(j,i)=0.

M(i,j)=0.

M(j,i)=0.

enddo

K(i,i)=1.

M(i,i)=1.

enddo

do i=(elem+1)*2,elem**2+elem,elem+1

do j=1,(elem+1)**2

K(i,j)=0.

K(j,i)=0.

M(i,j)=0.

M(j,i)=0.

enddo

K(i,i)=1.

M(i,i)=1.

enddo

do i=elem**2+elem,(elem+1)**2

do j=1,(elem+1)**2

K(i,j)=0.

K(j,i)=0.

M(i,j)=0.

M(j,i)=0.

enddo

K(i,i)=1.

M(i,i)=1.

enddo

K=K*210.E9*6./(7800.*delta**2.)

call dgvcsp((elem+1)**2,K,(elem+1)**2,M,(elem+1)**2,EVAL,EVEC,(elem+1)**2)

write(*,*)"Modos de vibracion"

write(*,'(F15.3)')dsqrt(eval)/(2.*3.141592654)

pause

write(*,*)"Error porcentual"

write(*,*)"1 modo",(dabs(dsqrt(eval((elem+1)**2))/(3.141592654*2.)-3350.)/3350.)*100.

write(*,*)"2 modo",(dabs(dsqrt(eval((elem+1)**2-2))/(3.141592654*2.)-5607.)/5607.)*100.

write(*,*)"3 modo",(dabs(dsqrt(eval((elem+1)**2-4))/(3.141592654*2.)-6875.)/6875.)*100.

write(*,*)"4 modo",(dabs(dsqrt(eval((elem+1)**2-6))/(3.141592654*2.)-8218.)/8218.)*100.

pause

write(*,*)"Desplazamientos nodales"

call dwrrrl('',(elem+1)**2,(elem+1)**2,evec,(elem+1)**2,0,'(F6.3)','','')

open(1,file="placa.dat")

write(1,'(F12.5)')evec(:,118)

end program EF

5.6. Código fuente para la chapa con 4 vínculos con cambio en densidad

program EF

use imsl

implicit none

real(8)::delta,l

real(8),dimension(4,4)::a,b,b2

real(8),allocatable::K(:,:),M(:,:),evec(:,:),eval(:)

integer(4)::elem,i,nex,ney,j,p,q

write(*,*)"Ingrese cantidad de elementos y el lado"

read(*,*)elem,l

delta=l/elem

allocate(K((elem+1)**2,(elem+1)**2),M((elem+1)**2,(elem+1)**2),eval((elem+1)**2),evec((elem+1)**2,(elem+1)**2))

do i=1,4

a(i,i)=4.

b(i,i)=4.

enddo

do i=1,3

a(i,i+1)=-1.

a(i+1,i)=-1.

b(i,i+1)=2.

b(i+1,i)=2.

enddo

do i=1,2

a(i,i+2)=-2.

a(i+2,i)=-2.

b(i,i+2)=1.

b(i+2,i)=1.

enddo

a(1,4)=-1.

a(4,1)=-1.

b(1,4)=2.

b(4,1)=2.

do nex=1,elem

do ney=1,elem

b2=b

if (((nex.eq.(elem-1)).or.(nex.eq.(elem-2)).or.(nex.eq.(elem-3))).and.((ney.eq.(elem- 1)).or.(ney.eq.(elem-2)).or.(ney.eq.(elem-3)))) then

b2=b*0.346

endif

do i=1,4

do j=1,4

if (i.eq.1) then

p=nex+(ney-1)*(elem+1)

elseif (i.eq.2) then

p=nex+1+(ney-1)*(elem+1)

elseif (i.eq.3) then

p=nex+(elem+1)+(ney-1)*(elem+1)

elseif (i.eq.4) then

p=nex+(elem+2)+(ney-1)*(elem+1)

endif

if (j.eq.1) then

q=nex+(ney-1)*(elem+1)

elseif (j.eq.2) then

q=nex+1+(ney-1)*(elem+1)

elseif (j.eq.3) then

q=nex+(elem+1)+(ney-1)*(elem+1)

elseif (j.eq.4) then

q=nex+(elem+2)+(ney-1)*(elem+1)

endif

K(p,q)=K(p,q)+a(i,j)

M(p,q)=M(p,q)+b2(i,j)

enddo

enddo

enddo

enddo

do i=1,elem+2

do j=1,(elem+1)**2

K(i,j)=0.

K(j,i)=0.

M(i,j)=0.

M(j,i)=0.

enddo

K(i,i)=1.

M(i,i)=1.

enddo

do i=elem+2,elem**2+elem,elem+1

do j=1,(elem+1)**2

K(i,j)=0.

K(j,i)=0.

M(i,j)=0.

M(j,i)=0.

enddo

K(i,i)=1.

M(i,i)=1.

enddo

do i=(elem+1)*2,elem**2+elem,elem+1

do j=1,(elem+1)**2

K(i,j)=0.

K(j,i)=0.

M(i,j)=0.

M(j,i)=0.

enddo

K(i,i)=1.

M(i,i)=1.

enddo

do i=elem**2+elem,(elem+1)**2

do j=1,(elem+1)**2

K(i,j)=0.

K(j,i)=0.

M(i,j)=0.

M(j,i)=0.

enddo

K(i,i)=1.

M(i,i)=1.

enddo

K=K*210.E9*6./(7800.*delta**2.)

call dgvcsp((elem+1)**2,K,(elem+1)**2,M,(elem+1)**2,EVAL,EVEC,(elem+1)**2)

write(*,*)"Modos de vibracion"

write(*,'(F15.3)')dsqrt(eval)/(2.*3.141592654)

pause

write(*,*)"Desplazamientos nodales"

call dwrrrl('',(elem+1)**2,(elem+1)**2,evec,(elem+1)**2,0,'(F6.3)','','')

open(1,file="placa.dat")

write(1,'(F12.5)')evec(:,36)

end program EF

Matías Daniel Vigliano

Noviembre 2005

COMISIÓN NACIONAL DE ENERGÍA ATÓMICA

UNIVERSIDAD NACIONAL GRAL SAN MARTIN

"Instituto de Tecnología Prof. Jorge Sabato"

Ingeniería en Materiales

Modelización