Monografias.com > Ingeniería
Descargar Imprimir Comentar Ver trabajos relacionados

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




Enviado por vigliano



    1. Barra
    2. Chapa
    3. Bibliografía
      consultada
    4. 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.

      1. 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.

      2. Funciones
        de interpolación
      3. 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
    1. 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

    Nota al lector: es posible que esta página no contenga todos los componentes del trabajo original (pies de página, avanzadas formulas matemáticas, esquemas o tablas complejas, etc.). Recuerde que para ver el trabajo en su versión original completa, puede descargarlo desde el menú superior.

    Todos los documentos disponibles en este sitio expresan los puntos de vista de sus respectivos autores y no de Monografias.com. El objetivo de Monografias.com es poner el conocimiento a disposición de toda su comunidad. Queda bajo la responsabilidad de cada lector el eventual uso que se le de a esta información. Asimismo, es obligatoria la cita del autor del contenido y de Monografias.com como fuentes de información.

    Categorias
    Newsletter