Vibraciones libres no amortiguadas en barras y chapas delgadas. Método de los elementos finitos
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 aEsta 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
comodonde 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 obtieneque substituyendo en la ecuación de
movimiento reducida, se llega aReordenando, 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.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.- Funciones
de interpolación - 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 quedacon 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
quedacon 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.- 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
voladizoEl 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 barraFigura 3. Función 2 de
interpolación lineal para la
barraUna 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
voladizoSe 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ínculosSe 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
barraFigura 7. Función 2 de
interpolación cúbica para la barraFigura 8. Función 3 de
interpolación cúbica para la
barraFigura 9. Función 4 de
interpolación cúbica para la barraUna 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
voladizoSe 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úbicasEl error porcentual se calculó comparando con
la expresión exacta de frecuencias
naturalescon 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. - 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
voladizoFigura 13. Función 1 de
interpolación para la
chapaFigura 14. Función 2 de
interpolación para la
chapaFigura 15. Función 3 de
interpolación para la
chapaFigura 16. Función 4 de
interpolación para la
chapaUna 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 HzFigura 18. Desplazamientos
nodales para la chapa en voladizo para la segunda frecuencia
natural de 3178 HzFigura 19. Desplazamientos
nodales para la chapa en voladizo para la tercera frecuencia
natural de 5296 HzFigura 20. Desplazamientos
nodales para la chapa en voladizo para la cuarta frecuencia
natural de 5902 HzSe 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ínculosResolviendo 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 HzFigura 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 HzSe 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 voladizoFigura 27. 4º frecuencia
natural en función del número de elementos para
chapa en voladizoFigura 28. 1º frecuencia
natural en función del número de elementos para
chapa con 4 vínculosFigura 29. 4º frecuencia
natural en función del número de elementos para
chapa con 4 vínculosA 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
densidadPor ú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 HzFigura 33. Curvas de nivel
para la chapa con aluminio con 4 vínculos a una
frecuencia de 4990 HzEn 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. - Chapa
- 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
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