MODELAMIENTO COMPUTACIONAL POR ELEMENTOS FINITOS PARA EL ESTUDIO DE RECIPIENTES DE PLÁSTICO REFORZADO CON FIBRA DE VIDRIO SOMETIDOS A CARGAS IMPULSIVAS CARLOS HERRERA VILLA YURY DANIEL MORA GONZÁLEZ UNIVERSIDAD TECNOLÓGICA DE BOLÍVAR FACULTAD DE INGENIERÍA MECÁNICA CARTAGENA DE INDIAS DT Y C. 2003 MODELAMIENTO COMPUTACIONAL POR ELEMENTOS FINITOS PARA EL ESTUDIO DE RECIPIENTES DE PLÁSTICO REFORZADO CON FIBRA DE VIDRIO SOMETIDOS A CARGAS IMPULSIVAS CARLOS HERRERA VILLA YURY DANIEL MORA GONZÁLEZ Trabajo de grado, presentado para optar al título de Ingeniero Mecánico DIRECTOR JAIRO F. USECHE VIVERO ING. MECANICO UNIVERSIDAD TECNOLÓGICA DE BOLÍVAR FACULTAD DE INGENIERÍA MECÁNICA CARTAGENA DE INDIAS DT Y C. 2003 _____________________________________ _____________________________________ _____________________________________ _____________________________________ _____________________________________ _____________________________________ _____________________________________ Firma del presidente del jurado _____________________________________ Firma de l jurado _____________________________________ Firma del jurado Cartagena, de 2003 A Dios, a mis padres, hermanos y a todas aquellas personas de confiaron en que esta meta pudo realizarse. Yury Daniel Mora G. A Jehová, a mis padres y amigos que brindaron su apoyo y confianza para alcanzar esta meta en mi vida Carlos A. Herrera Villa AGRADECIMENTOS Especial agradecimiento al ingeniero Jairo Useche Vivero por su gran colaboración durante todo el transcurso de este proyecto. También agradecemos a Rotofibra Ltda. por toda su participación en este trabajo. CONTENIDO Pág. INTRODUCCION 1. CONCEPTOS BASICOS 4 1.1 CONDICIONES DE EQUILIBRIO 4 1.2 ESFUERZO 5 1.2.1 Componentes de Esfuerzo 6 1.3 RESULTANTE DE FUERZA INTERNA 8 1.4 ECUACIONES DIFERENCIALES DE EQUILIBRIO 10 1.5 PRINCIPIO DE SUPERPOSICIÓN 12 1.6 RELACIONES DESPLAZAMIENTO – DEFORMACION 12 1.6.1 Ecuaciones de compatibilidad. 16 2. TEORÍA DE PLACAS 18 2.1 COMPORTAMIENTO GENERAL DE PLACAS 19 2.2 RELACIONES DE DEFORMACIÓN – CURVATURA 20 2.3 ESFUERZOS Y RESULTANTES DE ESFUERZOS 22 2.4 VARIACIÓN DEL ESFUERZO DENTRO DE UNA PLACA 24 2.5 ECUACIÓN PARA LA DEFLEXIÓN DE PLACAS 27 2.6 ENERGIA DE DEFORMACION DE LA PLACA 28 2.7 PRINCIPIO DE TRABAJO VIRTUAL 29 2.8 PRINCIPIO DE LA ENERGÍA POTENCIAL MINIMA 30 2.9 PLACAS ANISOTRÓPICAS 31 2.9.1 Relaciones Básicas 32 2.9.2 Determinación de la matriz de elasticidad. 34 2.9.3 Placas de multicapas isotrópicas 37 2.9.4 Placas compuestas laminadas 39 3. FORMULACION POR ELEMENTOS FINITOS 42 3.1 ANALISIS DE PLACAS EN EL PLANO 42 3.2 ANALISIS DE PLACAS A FLEXION 44 3.3 MATRIZ COMBINADA 47 3.4 ELEMENTOS FINITOS 48 3.4.1 Transformación isoparamétrica. 49 3.4.2 Elementos conformes 52 3.5 TRANSFORMACIÓN DE COORDENADAS GLOBALES 55 3.5.1 Elementos orientados en forma arbitraria 57 3.6 INTEGRACIÓN NUMÉRICA 58 3.6.1 Integrales bidimensionales 60 3.7 AJUSTE POR MINIMOS CUADRADOS PARA UN CUADRILATERO DE CUATRO NODOS 61 3.8 ANÁLISIS DINÁMICO 63 3.8.1 Comportamiento Dinámico de las estructuras elásticas con amortiguamiento lineal. 63 3.8.2 Método de integración del tiempo de Newmark 64 4 VALIDACION DEL MODELO 67 4.1 DESCRIPCION DEL MODELO COMPUTACIONAL 67 4.2 VALIDACION ANALITICA 67 4.2.1 Viga en voladizo cargada puntualmente en el extremo libre 68 4.2.2 Viga en voladizo con carga distribuida uniforme 70 4.2.3 Viga en voladizo sometida a carga axial 72 4.2.4 Cilindro recto sometido a compresión 74 4.2.5 Placa sujeta a una carga uniforme 76 4.2.6 Placa laminada sujeta a una carga uniforme 78 4.2.7 Análisis de resultados 87 4.3 VALIDACION DINAMICA CON UN MODELO REAL 88 4.3.1 Validación dinámica en una lámina de PRFV 89 4.3.2 Calibración de los coeficientes de Rayleigh 91 5 APLICACIÓN DEL MODELO 94 5.1 DESCRIPCION DEL ENSAYO 94 5.2 MATERIALES DEL TANQUE 94 5.2.1 VFG-Woven Roving 800 95 5.2.2 Mat 450 95 5.2.3 Cristalán 805 95 5.3 SIMULACION ESTATICA 95 5.3.1. Resultados del análisis. 97 5.4 SIMULACION DINAMICA 99 5.4.1 Análisis de resultados 101 5.5 ESTUDIO COMPARATIVO 102 5.5.1 Análisis de resultados 102 CONCLUSIONES 104 RECOMENDACIONES 106 BIBLIOGRAFIAS 107 ANEXOS 109 LISTA DE FIGURAS Pág. Figura 1. Diagrama de cuerpo libre de la sección cortada 6 Figura 2. Esfuerzos en un elemento en 3D 6 Figura 3. Representación en dos dimensiones del elemento de esfuerzos 7 Figura 4. Resultante de las fuerzas internas y componentes de dF 8 Figura 5. Elementos con esfuerzos y fuerzas de cuerpo 10 Figura 6. Esfuerzo normal en una barra 13 Figura 7. Deformaciones en el elemento 2D 14 Figura 8. Comportamiento de una placa 19 Figura 9. Resultante de esfuerzos positivos y cargas en el elemento placa 25 Figura 10. Lamina ortotrópica 34 Figura 11. Placa laminada simétricamente construida 37 Figura 12. Laminado compuesto 39 Figura 13. Elemento Transformado 50 Figura 14. Coordenadas locales y globales 55 Figura 15. Coordenadas locales y globales para un elemento cuadrilátero transformado 57 Figura 16. Cuadratura Gaussiana en dos dimensiones 61 Figura 17. Ajuste por mínimos cuadrados para un cuadrilátero 61 Figura 18. Viga en voladizo cargada puntualmente en el extremo libre. 68 Figura 19. Número de Elementos vs % de Error 69 Figura 20. Número de Elementos vs Deflexión 69 Figura 21. Viga en voladizo con carga distribuida uniforme 70 Figura 22. Número de Elementos vs % de Error 71 Figura 23. Número de Elementos vs Deflexión 71 Figura 24. Viga en voladizo sometida a carga axial 72 Figura 25. Número de Elementos vs % de Error 73 Figura 26. Número de Elementos vs Deflexión 73 Figura 27. Cilindro recto 74 Figura 28. Desplazamientos en dirección Y con enmallado de 512 elementos 75 Figura 29. Desplazamientos totales con enmallado de 512 elementos 75 Figura 30. Desplazamientos en dirección Z 77 Figura 31. Giro totales 77 Figura 32. Placa laminada sujeta a una carga uniforme 78 Figura 33. Desplazamiento en Z 83 Figura 34. Giros Totales 83 Figura 35. Esfuerzos en dirección X en la capa 1 84 Figura 36. Esfuerzos en dirección X en la capa 5 84 Figura 37. Esfuerzos en dirección Y en la capa 1 85 Figura 38. Esfuerzos en dirección Y en la capa 5 85 Figura 39. Esfuerzo cortante XY en la capa 1 86 Figura 40. Esfuerzo cortante XY en la capa 5 86 Figura 41. Estado inicial de la placa 88 Figura 42. Comportamiento dinámico de la placa 88 Figura 43. Adquisición de datos 89 Figura 44. Variación del desplazamiento en el tiempo. 91 Figura 45. Comparación de resultados. 93 Figura 46. Distribución de fuerzas en tanque de 10000 Litros 96 Figura 47. Desplazamientos en dirección X y Y 97 Figura 48. Desplazamientos totales 97 Figura 49. Esfuerzos normales 98 Figura 50. Esfuerzos cortantes 98 Figura 51. Desplazamiento vs tiempo en la zona inferior 99 Figura 52. Desplazamiento vs tiempo en la zona media 100 Figura 53. Desplazamiento vs tiempo en la zona superior 100 LISTA DE TABLAS Pág. Tabla 1. Funciones de interpolación de Hermite φ 53 Tabla 2. Puntos de Gauss y Pesos para cuadratura Gaussiana 59 Tabla 3. Resultados del análisis por elementos finitos 68 Tabla 4. Resultados del análisis por elementos finitos 70 Tabla 5. Resultados del análisis por elementos finitos 72 Tabla 6. Comparación de resultados 74 Tabla 7. Aproximación de la ecuación de la deflexión máxima 76 Tabla 8. Resultados utilizados diferentes mallas 77 Tabla 9. Valores de los coeficientes α como una función de n y m 79 Tabla 10. Valores de los coeficientes β y γ como función de m y n 79 Tabla 11. Variación de esfuerzos en el centro de la placa. 81 Tabla 12. Resultados del análisis por elementos finitos 82 Tabla 13. Datos de la lámina 90 Tabla 14. Resultados de la prueba 90 Tabla 15. Resultados de la pruebas 93 Tabla 16. Datos del tanque 96 Tabla 17. Comparación de resultados 102 LISTA DE ANEXOS Pág. ANEXO A. TRANSFORMACION DEL DIFERENCIAL DE AREA (dA) 109 ANEXO B. SOLUCION DE NAVIER PARA UNA PLACA RECTANGULAR SOPORTADA SIMPLEMENTE Y CARGADA UNIFORMEMENTE 113 ANEXO C. SOLUCION DE UNA PLACA LAMINADA SOPORTADAS SIMPLEMENTE Y SUJETA A UNA CARGA UNIFORME 116 ANEXO D. CODIGO DE PROGRAMCION EN MATLAB 120 RESUMEN La formulación por elementos finitos para el estudio de tanques fabricados en plásticos reforzados en fibra de vidrio se comenzará desarrollando un elemento de placa que tenga la capacidad de analizar fuerzas en las tres direcciones x, y, z. Este elemento es un cuadrilátero de cuatro nodos que puede ser no simétrico, el cual tiene la capacidad de calcular desplazamiento y giros en las tres direcciones de los ejes coordenados. Conociendo que las paredes del tanque están fabricadas con material compuesto y que estas se conforman por capas, las cuales tienen diferentes propiedades, el estudio se hace por medio de un elemento laminado y la matriz elástica de este se halla de forma diferente a la de un elemento no laminado. No obstante este elemento sufre una transformación de coordenadas para que funcione en cualquier geometría. Con los anteriores análisis se pueden calcular la matriz de rigidez para cada elemento y con esta conformar la matriz de rigidez global. Y así, poder calcular desplazamiento, giros y esfuerzos que se presentan en toda la geometría. Los análisis anteriores hacen parte del análisis estático, sin embargo para el estudio dinámico se hará uso de la matriz de rigidez descrita anteriormente y adicionalmente las matrices de masa y amortiguamiento, las cuales depende directamente de las funciones de formas del elemento. Con estas matrices y el vector de fuerzas se pueden calcular la variación del desplazamiento, velocidad y aceleración en función del tiempo. Estos cálculos se logran gracias a método de integración del tiempo de Newmark. INTRODUCCION Hoy en día los Plásticos Reforzados con Fibra de Vidrio (PRFV) desempeñan un papel fundamental en la fabricación de piezas y recipientes para la industria, ya que esto presentan las ventajas de poco peso respecto a los metales, son más resistentes a la corrosión y buena resistencia mecánica. Estas y muchas más ventajas son lo que hacen que los PRFV sean una excelente alternativa a la hora de fabricar diferentes productos. Aunque los PRFV presentan muchas ventajas en su fabricación y a lo largo de su vida útil, su diseño convencional puede dejar ciertas incertidumbres al diseñador ya que este es un material compuesto que combina diferentes propiedades. Además la gran mayoría de las empresas no cuentan con herramientas o procedimientos que le permitan optimizar y evaluar sus diseños. En vista de la falta de herramientas que permitan evaluar y optimizar los diseños se hace necesario incorporar en la línea de producción de estas empresas una herramienta que lo permita. En la actualidad existen dichas herramientas como son softwares especializados que permiten al diseñador evaluar su diseño desde un punto de vista ingenieril, además traen consigo ahorro de tiempo y dinero. El análisis estructural de estos softwares se basa en métodos numéricos, como es el caso del método de elementos finitos. El método de elementos finitos es un método numérico que facilita la solución de problemas de ingeniería. Analizando un problema específico como es el caso de tanque de paredes delgadas, que es el problema a analizar en este trabajo de grado. El cual pretende estudiar el comportamiento elástico en tanque fabricados en PRFV cuando son sometidos a una carga estática o dinámica, teniendo presente que se omitirán los efectos por temperatura ya que las cambio de temperatura durante la operación de estos tanques son mínimas y se mantienen en un rango de temperaturas (temperatura ambiente) donde las propiedades del material no cambian, por tanto sus efectos son despreciables para el estudio. El estudio comprende el desarrollo de desplazamientos, giros y esfuerzos en los tanques cuando son expuestos a las cargas ya mencionadas. El tener esta información es de gran importancia para el diseñador ya que le permite observar el comportamiento real de su diseño, y así efectuar cambios si son necesarios de acuerdo a los resultados. Sin embargo, los resultados arrojados por modelamiento numérico están ligado a un porcentaje de error, que dependerán de la cantidad de factores de diseño omitidos en el análisis, el tipo de elemento, método de soluciones de las ecuaciones y cantidad de elementos. No obstante, es necesario aclarar que estos errores son bastante pequeños y no desvirtúan los resultados de manera considerable, de tal forma la solución es confiable. 1. CONCEPTOS BASICOS 1.1 CONDICIONES DE EQUILIBRIO Fuerzas externas, son las que actúan sobre el cuerpo y pueden clasificarse en fuerzas de superficie y fuerzas de cuerpo. Las fuerzas de cuerpo, están asociadas con el cuerpo en estudio y se distribuyen en toda la amplitud del mismo. No son consecuencia de un contacto directo con otros cuerpos. Se especifican en términos de fuerzas por unidad de volumen y entre ellas podemos citar las fuerzas gravitacionales, las de inercia, las magnéticas, entre otras. Las componentes de la intensidad de estas fuerzas según los ejes coordenados, las denominaremos. Fx, Fy, y Fz. Las fuerzas de superficie, son una consecuencia del contacto físico entre dos cuerpos. Las fuerzas internas, son las fuerzas que están junto con las partículas que conforman el elemento. Cuando en un sistema de fuerzas, que actúa sobre un cuerpo tiene como resultante cero, el cuerpo se dice que esta en equilibrio. El equilibrio de fuerzas es el estado en el cual las fuerzas aplicadas son balanceadas. En un problema tridimensional, las condiciones de equilibrio requiere el cumplimiento de las siguientes ecuaciones de estática: ∑ Fx = 0 ∑Fy = 0 ∑ Fz = 0 1.1 ∑M x = 0 ∑M y = 0 ∑M z = 0 En un problema en el plano, donde todas las fuerzas actúan sobre el plano xy, hay solo tres ecuaciones de estática: ∑ Fx = 0 ∑ Fy = 0 ∑M A = 0 1.2 donde A, B y C son puntos sobre el plano. Reemplazando una sumatoria de fuerzas por un equivalente sumatoria de momento en la ecuación anterior, obtenemos otra alternativa a la condición de equilibrio: ∑ Fx = 0 ∑M A = 0 ∑M B = 0 1.2a Asumiendo que la línea de conectividad entre los puntos A y B no es perpendicular al eje x: ∑M A = 0 ∑M B = 0 ∑MC = 0 1.2b donde los puntos A, B, y C no son colineales. Una estructura se dice que es estáticamente determinada cuando todas las fuerzas en los miembros pueden obtenerse usando solo las condiciones de equilibrio. En el caso contrario, una estructura es llamada estáticamente indeterminada. 1.2 ESFUERZO Considere un cuerpo en equilibrio, sujeto al sistema de fuerzas externas. Con el fin de investigar las fuerzas interiores en algún punto Q, se realiza un corte en Q por un plano imaginario que divide el cuerpo en dos porciones. El equilibrio de fuerzas que actúan en una porción requiere la presencia de fuerzas interiores en el lugar de la sección. Estas fuerzas interiores aplicadas a ambas porciones, son continuamente distribuidas sobre la superficie cortada. Figura 1. Diagrama de cuerpo libre de la sección cortada Fuerzas Externas Fuerzas Internas ∆A x Q ∆F y z Sobre la superficie cortada, en el punto Q se localiza un ∆A, la cual es influenciada por una fuerza ∆F. En general ∆A no queda a lo largo de “x”, “y” y “z”. Donde sus componentes son ∆Fx normal, ∆Fy, ∆Fz tangente a ∆A. Entonces el esfuerzo normal σx, y el esfuerzo cortante τxy, τxz pueden definirse como: ∆Fx dFx ∆Flim ; lim y dFy ; lim ∆Fz dFσ x = = τ x y = = τ = = z 1.3 ∆A→0 ∆A dA ∆A→0 ∆A dA x z ∆A→0 ∆A dA 1.2.1 Componentes de Esfuerzo Figura 2. Esfuerzos en un elemento en 3D σy y x τ yx z τ xy τ xz τ yz τ σzx x τ zy σz Un cubo de dimensiones infinitesimales aislado de un sólido expondría el caso general del estado tridimensional de esfuerzos. Se considera que el esfuerzo mostrado es el mismo en las caras mutuamente paralelas, y uniformemente distribuida en cada cara. Sin embargo, en general los esfuerzos variarían de una cara a una cara paralela y también variarían encima de una cara particular. Un total de nueve componentes de esfuerzo escalar, definiendo el estado de esfuerzo a un punto, puede escribirse en la siguiente forma: τ x x τ x y τ x z  σ x τ x y τ x z  τ = τ τ τ  =  i j  y x y y y z  τ y x σ y τ y z  1.4 τ τ τ  z x z y z z τ z x τ z y σ z  Esta es una representación en matriz del esfuerzo. La anotación del subíndice doble se interpreta como sigue: El primer subíndice denota la dirección de una normal a la cara en que la componente del esfuerzo actúa; el segundo denota la dirección del esfuerzo. En el caso de esfuerzo bidimensional, sólo las caras de “x” e “y” del elemento está sujeto a los esfuerzos, y todo los esfuerzos actúan paralelo al eje “xy”. Para conveniencia, dibujamos a menudo una vista bidimensional del elemento de esfuerzo plano. Cuando sólo dos esfuerzos normales están presentes, el estado de esfuerzo se denomina en término biaxial. Figura 3. Representación en dos dimensiones del elemento de esfuerzos y σ y x τ xy τ xy = τ yx σ σ x x τ xy τ yx σ y Convención de signos. Cuando la componente de esfuerzo actúa en una cara positiva en una dirección de la coordenada positiva, el componente de esfuerzo es positiva, también cuando actúa en una cara negativa en una dirección de la coordenada negativa. Pero un esfuerzo es considerado negativo cuando actúa en una cara positiva en dirección negativa de coordenada (o viceversa). De acuerdo, los esfuerzos de tensión siempre son positivos y los esfuerzos de compresión siempre son negativos. 1.3 RESULTANTE DE FUERZA INTERNA Las fuerzas distribuidas dentro de un miembro pueden representarse por una fuerza estáticamente equivalente y un vector del momento que actúan a cualquier punto arbitrario de una sección. Éstas resultantes de fuerzas internas, también llamadas resultantes de esfuerzo, consisten en fuerza axial, cortante y momentos. Estas son expuestas usualmente en un corte imaginario del plano incluyendo el centroide C a través del elemento y resultado en componentes normal y tangencial a la sección cortada. Figura 4. Resultante de las fuerzas internas y componentes de dF V y y C P T dA z x M y dFx dF Vz x y M z dFz z Cada fuerza interior y componente del momento refleja un efecto diferente de la carga aplicada en un miembro. Estos efectos son como sigue, con referencia a la figura 4: La fuerza axial P tiende a alargar o acortar al miembro. Fuerzas cortantes Vy y Vz tienden a cortar una parte del miembro relativa a la parte adyacente. El torque T es responsable para torcer al miembro. Los momentos de flexión M y y M z causan la flexión del miembro. La convención de signos adoptado para el esfuerzo también se aplica a la fuerza y componentes del momento. Para establecer la relación entre las componentes de los esfuerzos y las resultantes de fuerza internas, considere un área infinitesimal dA de la sección cortada mostrada en la figura 4. Este área ha actuado en las componentes de dF, como dFx = σ xdA , dFy = τ x ydA y dFz = τ x ydA Obviamente, los componentes de esfuerzos en la sección cortada causan las resultantes de fuerza interna en esa sección. Por consiguiente, la suma de las fuerzas increméntales en las direcciones de “x”, “y” y “z” son: P = ∫σ xdA Vy = ∫τ xydA Vz = ∫τ xzdA 1.5 Similarmente, las sumas de los momentos de las mismas fuerzas sobre los ejes “x”, “y” y “z” da como resultado: T = ∫ (τ x z y −τ x y z )dA M y = ∫σ x zdA M z = −∫σ x ydA 1.6 1.4 ECUACIONES DIFERENCIALES DE EQUILIBRIO Figura 5. Elementos con esfuerzos y fuerzas de cuerpo ∂σ σ yy + dy ∂y ∂τ τ y x y x + dy∂y ∂τ τ x y + x y dx Fy ∂x σ dy ∂σ σ + x dx x xF ∂xx τ x y dx y τ y x x σ y Para un caso bidimensional, la acción de esfuerzos en un elemento de lados, dx y dy , y de unidad de espesor se muestra en la figura 5. Las componentes “x” y “y” de las fuerzas del cuerpo por unidad de volumen, Fx y Fy , son independientes de la componente “z” de la fuerza del cuerpo Fz = 0 . Si establecemos el equilibrio de fuerzas en el sentido del eje “x”, considerando un valor unitario de “z”, tendremos:  ∂σ  ∂τ  Fx dx dy + σ x + x dx dy −σ dy + τ y x x  y z + dy ∂x −τ y xdx = 0  ∂x   ∂y  que una vez simplificada se convierte en: ∂σ x ∂τ  + y z + Fx  dx dy = 0  ∂x ∂y  Como el producto “dxdy” no es cero, debe serlo la expresión entre corchetes, por lo que, finalmente, obtenemos: ∂σ x ∂τ+ y z + Fx = 0 1.7 ∂x ∂y Si hacemos el mismo análisis para la dirección del eje Y: ∂σ y ∂τ+ x y + Fy = 0 1.8 ∂y ∂x Estas ecuaciones de equilibrio se pueden generalizar considerando el plano tridimensional ∂σ ∂τ y z ∂τx + + z x + Fx = 0∂x ∂y ∂z ∂σ y ∂τ ∂τ+ x y + z y + F = 0 1.9 ∂x ∂y ∂z y ∂σ z ∂τ x z ∂τ+ + y z + F = 0 ∂x ∂y ∂z z En consecuencia, para un cuerpo en equilibrio los esfuerzos varían de punto a punto según las ecuaciones anteriores. Podemos aplicar la condición de equilibrio ∑M z = 0  ∂σ y dy dx dx −  ∂σ x  ∂τ   ∂τ     dx dy  dy  + x y τ x y + dx  dy dx − τ y x + x y dy dx dy  ∂y  2  ∂x  2  ∂x    ∂y  +Fydx dy dx dy − Fxdy dx = 02 2 Despreciando los términos que contengan triples productos de dx y dy, la ecuación anterior se convierte en: τ x y =τ y x Considerando el caso tridimensional y tomando momentos respecto a cada uno de los ejes veríamos que: τ x y = τ y x τ x z = τ z x τ y z = τ z y por lo que podemos decir que sólo seis de las nueve componentes de la esfuerzo. 1.5 PRINCIPIO DE SUPERPOSICIÓN Este expresa que el estado de equilibrio debido a varias acciones exteriores es igual a la superposición de las soluciones que corresponden a cada uno de los estados si cada acción exterior actuara independientemente. Esta regla es válida siempre que la cantidad (deformación o esfuerzo) obtenida sea directamente proporcional a las cargas aplicadas. 1.6 RELACIONES DESPLAZAMIENTO – DEFORMACION Se dice que un cuerpo está deformado cuando las posiciones relativas de sus puntos han cambiado. En el movimiento como sólido rígido, dichas posiciones permanecen estables. Cuando se aplican fuerzas exteriores a un cuerpo, la posición de cada punto, en general, se modifica. Se define el desplazamiento de un punto como el vector que une el punto original con el desplazado. Se denomina a las componentes “x”, “y”, “z” del desplazamiento con las letras u, v, w respectivamente. Por tanto, un punto que estuviera inicialmente en la posición (x, y, z) se moverá al punto (x+u, y+v, z+w). En general, u, v, w serán función de x, y, z. Se considerará un modelo unidimensional para entender claramente el concepto de deformación. Figura 6. Esfuerzo normal en una barra x dx A B (a) Estado Original u + ∂u dx ∂x u A’ B’ (b) Estado De Deformación du ∂u+ dx ∂x En la figura 6 puede verse una barra sometida a una fuerza axial. Inicialmente, los puntos A y B están separados una distancia dx. Dichos puntos, bajo el efecto de la fuerza se desplazan a los puntos A' y B', y vemos que la distancia entre ambos ha aumentado ligeramente. Si definimos la deformación como el cambio unitario de longitud, se tiene: (∂u ∂x) dx ∂u ε x = = dx ∂x Se considerará un cuerpo en un estado de deformación plana, que se define por: u = u ( x, y) v = v ( x, y) w = 0 En este caso, todos los puntos del plano “xy”, permanecen en el plano después de la deformación. La deformación costa de dos tipos normal y tangencial. La deformación normal ε en una dirección dada se define como el cambio unitario de longitud de una línea que estaba originalmente orientada según la mencionada dirección. Es positiva si el cambio en la longitud consiste en un alargamiento, y negativa si se trata de un acortamiento. Figura 7. Deformaciones en el elemento 2D y ∂u dy ∂y C ' D ' ∂v dy ∂y C β D λ B ' A ' dy v ∂vθ dx∂x A B dx ∂u dx dx ∂x y x u x La deformación tangencial está asociada con dos direcciones y se define como el cambio en el ángulo recto original entre dos ejes (en radianes). Es positiva si el ángulo original decrece. El signo que se le da depende del sistema de coordenadas. A 'B '− AB A ' B '− dx ε x = =AB dx A 'D '− AD A ' D '− dy ε y = = AD dy π γ x y = − β = θ − λ2 donde el signo negativo de λ se basa en el hecho de que, para los giros, son positivos aquéllos que van en sentido contrario a las agujas del reloj. Si las componentes del desplazamiento del punto A son “u” y “v”, el punto B se desplazará u + (∂u ∂x) dx y v + (∂v ∂x) dx ya que “y” es constante a lo largo de la línea AB. Análogamente, las componentes del desplazamiento del punto D son u + (∂u ∂y) dy y v + (∂v ∂y) dy . En consecuencia se puede escribir: 2 2 ( A 'B ')2 = dx (1+ ε x ) 2  dx ∂u   ∂v =   +  +  dx   ∂x   ∂x  de forma que: ∂u ∂u 2 ∂v 2 ε 2 + 2ε +1 = 1+ 2  +  x     +    ∂x ∂x     ∂x  Pero como se consideran deformaciones infinitesimales, los términos elevados al cuadrado son despreciables (en un orden de magnitud) y se puede escribir finalmente: ∂u ∂v ε x = ε y = x y 1.10 ∂ ∂ Por otro lado, también se podrá escribir de acuerdo con la figura 6 (∂u ∂x) dx θ = dx + (∂u ∂x) dx ya que para desplazamientos infinitesimales se considerará la tangente igual al ángulo, es decir tanθ = θ . Por otro lado el segundo término del denominador también se puede despreciar, ya que ∂u ∂x << 1. Como lo mismo es cierto en el cálculo de λ, podemos escribir: ∂u ∂u θ = λ = − ∂x ∂y con lo que la deformación tangencial quedará: ∂u ∂v γ x y = + ∂y ∂x donde las dos derivadas parciales son positivas si AB y AD giran hacia adentro como es el caso de la figura. En el caso de un sistema tridimensional en el que el elemento original es un prisma rectangular, las componentes de la deformación serían: ∂u ∂u ∂v ε x = γ = +∂x x y ∂y ∂x ∂v ∂v ∂w ε y = γ y z = + 1.11 ∂y ∂z ∂y ∂w ∂w ∂u ε z = γ z x = +∂z ∂x ∂z donde se observa también el efecto simétrico con los esfuerzos: γ x y = γ y x γ y z = γ z y γ z x = γ x z 1.6.1 Ecuaciones de compatibilidad. Se analizará en detalle las ecuaciones 1.11 que se dedujeron. Se trata de un conjunto de seis ecuaciones para cada componente de la deformación, pero en función de sólo tres componentes del desplazamiento. Si se especifican estas últimas en función de x, y, z, se obtienen las deformaciones. Sin embargo se pude razonar a la inversa y partir del hecho de que las seis deformaciones sean funciones dadas de las mismas variable x, y, z.. En este caso se tienen seis ecuaciones para el cálculo de tres incógnitas u, v, w.. Ese sistema de ecuaciones sería indeterminado y no tendría solución para los desplazamientos, a no ser que las componentes de la deformación estuvieran relacionadas de alguna otra manera. En otras palabras, las componentes de las deformaciones no se pueden definir arbitrariamente, si se desea encontrar funciones de desplazamiento únicas y continuas. Por tanto, deben existir al menos tres ecuaciones adicionales que nos permitan establecer la relación biunívoca entre deformaciones y desplazamientos. Si se deriva la primera de las ecuaciones 1.11 dos veces con respecto a “y”, y la segunda dos veces respecto a “x” y sumamos los resultados, se tiene: ∂2ε 2x ∂ ε y ∂ 3u ∂3v 2 + 2 =∂y ∂x ∂y2 + ∂x ∂x2 ∂y Y si se deriva la cuarta ecuación respecto a “x” y respecto a “y, se tendrán: ∂2γ x y ∂2  ∂u ∂v  = + ∂x∂y ∂x∂y  ∂y ∂x   y como el orden de la derivación es indiferente: ∂2ε 2 ∂2γx ∂ ε+ y = x y 2 2 ∂y ∂x ∂x∂y De forma similar se pueden obtener cinco ecuaciones más. Las seis ecuaciones las se escriben a continuación: ∂2ε ∂2x ε y ∂ 2γ 2 + 2 = x y ∂y ∂x ∂x∂y ∂2ε ∂2ε ∂2y γ+ z = y z ∂z2 ∂y2 ∂y∂z ∂2ε ∂2ε ∂2x γz 2 + x z ∂x ∂z2 = ∂z∂x 1.12 2 ∂ 2ε x ∂  ∂γ ∂γ ∂γ = − y z + x z + x y ∂y∂z ∂x  ∂x ∂y ∂z   ∂2ε 2 y ∂  ∂γ = y z ∂γ − x z ∂γ x y + ∂z∂x ∂y  ∂x ∂y ∂z   ∂2ε z ∂  ∂γ ∂γ ∂γ 2 = y z + x z + x y ∂x∂y ∂z  ∂x ∂y ∂z   Estas ecuaciones se denominan: ecuaciones de compatibilidad de Saint-Venant . Las componentes de las deformaciones deben cumplir estas seis ecuaciones para que exista solución en desplazamientos. 2. TEORÍA DE PLACAS Las placas son los elementos estructurales inicialmente planos, con un espesor mucho más pequeño que las otras dimensiones. Se tratarán las placas de tal forma que sea común dividir el espesor t en las mitades de un mismo plano paralelo a las caras. Este plano recibe el nombre de plano medio de la placa. El espesor del plano es moderado en una dirección normal al medio plano. Las placas pueden ser clasificadas en tres grupos: placas delgadas con las desviaciones pequeñas, placas delgadas con las desviaciones grandes y placa de espesor. Según el criterio a menudo aplicado para definir una placa delgada (para los propósitos de cálculos técnicos) la proporción del espesor a la longitud del plano debe ser menos que 1/20. Consideraremos sólo desviaciones pequeñas de placas delgadas, unas simplificaciones consistentes con la magnitud de deformación normalmente encontrada en las estructuras de las placas. Como es notado previamente, a menos que por otra parte especificado, se asumirán que la placa y materiales de la cubierta son homogéneos e isotropicos. 2.2 COMPORTAMIENTO GENERAL DE PLACAS Figura 8. Comportamiento de una placa m x0 A t x x n z w t 2 y0 A m ∂w n ∂x z y r u z ∂wx = − ∂x (a) Placa de espesor constante (b) Parte de la placa antes y después de la flexión. Considere una placa libre de carga, mostrada en la figura 8a, en que el plano xy coincide con el medio plano y que la desviación de x es cero. Los componentes de desplazamiento a un punto, ocurriendo en las direcciones de x, y, y z, se denota por u v y w, respectivamente. Cuando debido a la carga lateral, la deformación tiene lugar, el medio plano a cualquier punto A(xa,ya) tiene desviación w (figura 8b). Los supuestos fundamentales de la teoría de pequeñas deformaciones doblamiento o llamada clásica o teoría de costumbre para isotrópicos, laminas homogéneas, elásticas, delgadas, esta basado en la geometría de deformaciones. Ellos pueden empezarse como: 1. La deformación del plano medio es pequeña comparada con el espesor de la placa. La inclinación de la superficie inclinada es muy pequeña y el cuadrado de la inclinación es despreciable comparado con la unidad. 2. El plano medio permanece libre subsecuente al pandeo. 3. La sección (mm) inicialmente normal al plano medio permanece plana y normal a la superficie después del pandeo. Esto quiere decir que las deformaciones angulares verticales γ x z y γ y z son despreciables. La deflexión de la placa esta asociada principalmente al pandeo de las deformaciones. De esto se deduce que la deformación normal ε z resultante de las cargas transversales también es omitida. 4. El esfuerzo normal al plano medio, σ z , es pequeño comparado con las otras componentes de esfuerzo y puede ser despreciable. Los supuestos anteriores se conocen como las hipótesis de Kirchhoff [1], son análogas a las deflexión de la teoría de vigas. En la mayoría de las aplicaciones de ingeniería, pueden encontrarse simplificaciones respecto al estado de deformación y esfuerzo. Esto resulta menos complejo, un problema tridimensional de placas reducido a dos dimensiones. Consecuentemente, las ecuaciones de gobierno de placas pueden ser derivadas de manera sencilla. 2.2 RELACIONES DE DEFORMACIÓN – CURVATURA Como consecuencia del supuesto 3, las relaciones de desplazamiento-deformación se reducen a: ∂u ∂v ∂u ∂v ε x = ε y = γ = +∂x ∂y x y ∂y ∂x 2.1 ∂w ε z = = 0 ∂v ∂w 0 ∂w ∂uγ y z = + = γ z x = + = 0∂z ∂z ∂y ∂x ∂z Integrando ε z , se obtiene: w = w(x, y) ____________ [1] UGURAL, Ansel C. Stresses in Plates and Shells. 2 ed. New Jersey: Mac Graw Hill, 1981. 502 p. Como las deflexiones laterales no varían sobre el espesor de la placa, integramos las expresiones γ x z y γ y z da como resultado u z ∂w= − + u0 ( x, y) v = −z ∂w + v ( x, y) ∂x ∂y 0 Donde u0 ( x, y) y v0 ( x, y) representan, respectivamente, los valores de u y v en el plano medio. Basándose en el supuesto 2 se puede concluir: u z ∂w ∂w= − v = −z x y 2.2 ∂ ∂ Sustituyendo al ecuación (2.2) en la (2.1), se tiene: 2 z ∂ w z ∂ 2w ∂2w ε x = − ε = − γ∂x2 y ∂y2 x y = −2z 2.3 ∂x∂y Estas ecuaciones son para cualquier punto en la placa. El radio de curvatura de un plano curvo es definido como la tasa de cambio del ángulo de declinación de la curva respecto a la distancia a lo largo de la curva. Tomando el supuesto 1 y derivando parcialmente la ecuación 2.3 representa las curvas de la placa. La curvatura κ (kappa) en la superficie media en el plano paralelo a los planos xz, yz, y xy son, respectivamente: 1 ∂  ∂w  1 ∂  ∂w  1 ∂  ∂w =   = κ x =  = κ = = κ r 2.4 x ∂x  ∂x  ry ∂y  ∂y  y  r  x y ∂x  ∂y  x y  donde κ x y = κ y x . La curvatura κ x y el radio de curvatura rx en la superficie media en el plano xz mostrado en la figura 8b. Similarmente, κ x y rr pueden estar en el plano yz. Las ecuaciones 3.4 en la declinación varía sobre la placa. La ultimas de estas expresiones también es llamada como torcedura del plano medio con respecto a los ejes “x” y “y”. Las relaciones de curvatura-deformación, están dadas por: ε x = −zκ x ε y = −zκ y τ x y = −2zκ x y 2.5 2.7 ESFUERZOS Y RESULTANTES DE ESFUERZOS En el caso de un estado tridimensional de esfuerzos, los esfuerzos y deformaciones relacionadas por la ley de Hooke generalizada, válida para materiales homogéneos isotrópico. 1 τ ε x yx = σ x −ν (σ y +σ z ) γ x y =E G 1 γ ε y = σ y −ν (σ x +σ z ) τ = x z E x z 2.6 G 1 τ ε z = σ z −ν (σ x +σ y ) τ y z E y z = G Donde τ ij = τ ji (i, j = x, y, z ) . Las constantes E, ν , G representan los módulos de elasticidad longitudinal o módulo de Young, radio de Possion y módulo de elasticidad transversal. G E= 2 (1+ν ) 2.7 Sustituyendo ε z = γ y z = γ x z = 0 en las ecuaciones 2.6, si obtiene: E E σ x = (ε +νε ) σ = (ε +νε ) τ = Gγ 1−ν 2 x y x 1−ν 2 y x x y x y Introduciendo las curvaturas de la placa, las ecuaciones 2.4 y 2.5 toma la forma: E z 2 2 σ x = − (κ +ν κ ) E z  ∂ w ∂ w = − +ν1−ν 2 x y 1−ν 2  ∂x2 ∂y2   E z E z  ∂2w ∂2( ) w σ y = − κ +ν κ = − +ν 2.8 1−ν 2 y x 1−ν 2  ∂y2 ∂x2   E z E z ∂2w τ x y = − 2 κ x y = −1−ν 1−ν 2 ∂x ∂y Los esfuerzos distribuidos sobre el espesor de la placa produce momentos flectores, momentos torsores y fuerzas cortantes. Hay momentos y fuerzas por unidad de longitud que también son llamados esfuerzos resultantes. t / 2 t / 2 ∫ zσ dy dz = dy zσ dz = M dy −t / 2 x ∫−t / 2 x Similarmente, las expresiones para los otros resultantes son derivados:  M x  σ x    t / 2  M y  =   ∫ σ y z dz 2.9 M  −t / 2    x y  τ x y  donde M x y = M y x y Qx  t / 2 γ   =Q  ∫ x z t / 2  dz 2.10  −y  γ y z  Es importante mencionar que la teoría de láminas delgadas omite los efectos de las componentes de deformación γ x z = τ x z G y γ y z = τ y z G en la flexión. Las fuerzas verticales Qx y Qy o son despreciables. Las ecuaciones de momentos flectores y torsores en términos de la curvatura y la deflexión, se obtiene sustituyendo la ecuación 2.8 en 2.9 M D ( ) D  ∂ 2w ∂2w  x = − κ x +ν κ y = −  +ν∂x2 ∂y2   2 2 M y = −D (κ y +ν κ x ) D ∂ w ∂ w = −  +ν∂y2 ∂x2  2.11   2 M x y = −D (1 ∂ w −ν )κ x y = −D (1−ν ) ∂x∂y donde 3 D Et= 12(1 2 ) 2.12 −ν es la rigidez flexional de la placa. Los signos negativos está de acuerdo con la convención de signos para el momento y curvatura. La fuerza de cortante vertical Qx y Qy se relacionan en w, y en la derivación de las ecuaciones de equilibrio. Cabe anotar que si un elemento placa, de ancho la unidad y con su eje paralelo al eje x, es libre de mover la cara bajo la carga transversal, la cima y superficie del fondo se deforma en forma de silla o antielástica, aparecería una curvatura κ y . La rigidez flexionante sería entonces Et3 12 . El resto de la lámina previene la curvatura, sin embargo; debido a esta acción, una parte de la lámina presenta rigidez mayor que una viga por un factor 1 (1−ν 2 ) , aproximadamente 10 por ciento. Las componentes de esfuerzo en dos dimensiones se encuentran ecuación 2.8 sustituyendo por 2.11 y empleando 2.12. De esta manera obtenemos 12M 12M z σ = x z 12M z x 3 σ y = y x y 3 τ x y = 3 2.13 t t t Los esfuerzos máximos ocurren en el fondo y superficie de la cima (en z = ± t 2 ) de la lámina. 2.8 VARIACIÓN DEL ESFUERZO DENTRO DE UNA PLACA Las componentes de esfuerzos generalmente varían de un punto a otro en la placa cargada. Estas variaciones son gobernadas por las condiciones de equilibrio de la estática. El cumplimiento de estas condiciones establece las relaciones conocidas como ecuaciones de equilibrio. Podríamos reducir el sistema de ecuaciones a una relación expresada en términos de momento. Figura 9. Resultante de esfuerzos positivos y cargas en el elemento placa Q M y y M x y Qx M ∂Mx M x + x dx ∂x M p x y ∂M ∂M M + x y dx M + yy dy x y ∂x ∂y ∂Q Qx + x dx ∂x ∂M M x yx y + dy ∂y ∂QQy + y dy ∂y Considerando un elemento dxdy de la placa sujeto a una distribución de cargas uniformes por unidad de área p. Asumiendo el peso de la placa como una cantidad pequeña, esto no afecta la exactitud de los resultados. Como el elemento es muy pequeño, para simplificar las componentes de fuerzas y momento puede considerarse una distribución uniforme en cada cara. En la figura 11 es mostrado por un vector que representa el valor de la carga aplicada en cada cara. Con cambiar la localización, como por ejemplo, del la esquina superior izquierda a la esquina inferior derecha, una de las componentes de momento, por ejemplo, Mx, actuando en la cara negativa x varía el valor relativo a la cara positiva x. Esta variación con la posición puede ser expresada por una expansión truncada de Taylor: M ∂M x dx x + ∂x La derivada parcial es usada porque Mx es una función de “x” y “y”. Tratando todas las componentes similarmente, el estado de resultante de esfuerzos mostrado en la figura es obtenido. La sumatoria de fuerzas en la dirección “z” es: ∂Qx ∂Q dx dy + y dx dy + p dx dy = 0 ∂x ∂y de la cual ∂Qx ∂Q+ y y + p = 0 2.14 ∂x ∂y El equilibrio de momento al eje “x” esta gobernado por: ∂M x y ∂M dx dy + y dx dy + Q dx dy = 0 ∂x ∂y y ∂M x y ∂M+ y + Qy = 0 2.15 ∂x ∂y Los términos de los productos infinitesimales tal como el momento de p y el momento debido al cambio en Qy han sido obtenidos. Similarmente, del equilibrio de momentos al eje y, tenemos: ∂M x y ∂M+ x + Q = 0 2.16 ∂y ∂x x Finalmente, la introducción de las expresiones para Qx y Qy de las ecuaciones 2.14 y 2.16 en 2.15: ∂2M x ∂ 2M ∂2x y M 2 + 2 + y 2 = − p 2.17 ∂x ∂x∂y ∂y Esta es la ecuación diferencial de equilibrio para flexión en placas delgadas. Las expresiones para las fuerzas cortantes Qx y Qy pueden ser escritas en términos de la deflexión w , de las ecuaciones 2.15 y 2.16 con 2.11:  2 2  Q ∂ ∂ w ∂ w ∂ 2x = −D + = −D ∇ w∂x  2 2  ( ) ∂x ∂y  ∂x 2.18  2 Q D ∂ ∂ w ∂ 2w  ∂ y = −  2 + 2  = −D ∇ 2w ∂y  ∂x y ( ) ∂  ∂y donde ∂2 ∂2 ∇2 = ∂x2 + ∂y2 este es el operador Laplace.[2] 2.9 ECUACIÓN PARA LA DEFLEXIÓN DE PLACAS La ecuación diferencial básica para la deflexión de placas puede obtenerse fácilmente introduciendo en la ecuación 2.17 las expresiones de M x , M y , M x y dadas en 2.11 y se obtiene: ∂2κ ∂2κ 2x + x y ∂ κ + y p 2 2 = 2.19 ∂x ∂x ∂y ∂y D Lo anterior fue la ecuación de equilibrio en términos de las curvaturas. Una forma alternativa de la ecuación 2.19 es determinada por la inserción de la definición de curvatura de la ecuación 2.4 ∂4w ∂4w ∂4w p 4 + 2 + = 2.20 ∂x ∂x2∂y y ∂y4 D ____________ [2] UGURAL, Ansel C. Stresses in Plates and Shells. 2 ed. New Jersey: Mac Graw Hill, 1981. 502 p. Esta ecuación, la primera derivada por Laplace en 1811, también puede escribirse como: 4 p∇ = 2.21 D 2 en la cual ∇4 = ∇2∇2 = (∇2 ) . Cuando no hay cargas laterales actuando en la placa, ∂4w ∂4w ∂42 w+ ∂x4 ∂x2∂y y + 4 = 0 2.22 ∂y La ecuación 2.20 es la ecuación diferencial de gobierno para deflexión de placas delgadas 2.10 ENERGIA DE DEFORMACION DE LA PLACA La energía de deformación en un cuerpo elástico para un estado general de esfuerzo [3] esta expresado por: U 1= ∫∫∫ (σ ε2 x x +σ y ε y +σ z ε z +τ x y γ x y +τ x z γ x z +τ y z γ y z )dx dy dz 2.23 V omitiendo los valores de σ z ,τ x z ,τ y z U 1= ∫∫∫ (σ x ε x +σ y ε y +τ x y γ x y )dx dy dz 2 V y tomando los valores de E E σ x = 2 (ε x +νε y ) σ x = (ε +νε ) τ = Gγ 1−ν 1−ν 2 y x x y x y y reduciendo la expresión y dejándola en función de esfuerzos y las constantes elásticas ____________ [3] UGURAL, Ansel C. Stresses in Plates and Shells. 2 ed. New Jersey: Mac Graw Hill, 1981. 502 p. U 1  1 1= ∫∫∫  σ 2 +σ 2 − 2νσ σ + τ 2 dx dy dz 2 V  2E ( x y x y ) 2G x y  2.24  Para una placa de espesor uniforme la ecuación 2.24 puede ser escrita en términos de la deflexión w y usando las ecuaciones 2.8 y 2.12 se tiene:  2 21  ∂ w   ∂2 2 2 w  ∂2w ∂2U D w  ∂2w   = ∫∫ 2  ∂x2  +  ∂y2  + 2ν ∂x2 2 + 2 (1−ν )  dx dy 2.25 A     ∂y  ∂x∂y    o 1  ∂2w ∂2 2 2 w   2 2U D ∂ w ∂ w  ∂2w   = ∫∫  2 + 2  − 2(1−ν )  2 2 −   dx dy 2.26 2 A  ∂x ∂y   ∂x ∂y  ∂x∂y   donde A es el área de la superficie. La ecuación 2.26 es conocida como curvatura Gaussiana. Se observa que la energía de deformación es no-lineal en función de la deformación o el esfuerzo. El principio de superposición no es valida para la energía de deformación. 2.10 PRINCIPIO DE TRABAJO VIRTUAL Suponga un cuerpo elástico el cual experimenta un incremento de desplazamiento virtual. Este desplazamiento no necesita realmente ocurrir y no necesita ser infinitesimal. Cuando el desplazamiento es tomado infinitesimal, como es frecuentemente hecho, es razonable considerar el sistema de fuerza actuando en el cuerpo como constante. El trabajo virtual [4] hecho por las fuerzas de superficie p por unidad de área en el cuerpo del estado inicial al estado de equilibrio es expresado como: ________________ [4] ORTIZ BERROCAL, Luis. Elasticidad. 3 ed. Madrid: Mac Graw Hill, 1998. 549 p. . δ w = ∫ ( pxδu + pyδ v + pzδ w)dA A donde A es la superficie de frontera y δ u, δ v, δ w son los desplazamientos virtuales en las direcciones x, y, z. La notación δ denota la cantidad de variación. La energía de deformación requerida de un cuerpo de volumen V como un resultado de la deformación virtual es: δU = ∫ (σ x δε x +σ y δε y +σ z δε z +τ x y δγ x y +τ x z δγ x z +τ y z δγ y z )dV 2.27 V El trabajo virtual hecho durante los desplazamientos virtuales son ceros: δW −δU = 0 . El principio de trabajo virtual para un cuerpo elástico es representado: δW = δU 2.11 PRINCIPIO DE LA ENERGÍA POTENCIAL MINIMA Considerando que el desplazamiento virtual no altera la forma del cuerpo y las fuerzas de superficie son estimadas como constantes, la ecuación puede escribirse como: δΠ = δ (U −W ) = 0 2.28 Esta expresión Π = U −W 2.29 Se muestra la energía del cuerpo. La ecuación 2.28 representa la condición de la energía potencial estacionaría del sistema. Esta puede mostrase para equilibrio estable la energía potencial es un mínimo. Para todo el desplazamiento satisfacen las condiciones de frontera y de equilibrio, la energía potencial será asumida un valor mínimo. Esto es referido como el principio de la energía potencial mínima. La energía potencial almacenada en la placa bajo una carga lateral distribuida p(x,y) es 1 Π = ∫∫∫ (σ xε x +σ yε y +τ x yγ x y ) dx dy dz − ∫∫ ( pw) dx dy 2 2.30 V A Para el caso de una placa de espesor constante 1 Π = 2 ∫∫ (M xκ x + M yκ y + M x yκ x y ) dx dy − ∫∫ ( pw) dx dy 2.31 A A Una explicación física de los términos de U en esta expresión es: como ∂2w ∂x2 = κ x representa la curvatura de la placa en el plano xz, el ángulo correspondiente el momento M 2 2y dy es igual a −(∂ x ∂x ) dx . La energía de deformación o el trabajo hecho por los momentos M y dy es 1 − M xκ x dx dy . La 2 energía de deformación debida a M y dx y M x y dy son interpretadas similarmente. El principio de la energía potencial, referido en la ecuación 2.31 es expresada en la forma: δΠ = −∫∫ (M xδκ x + M yδκ y + M x yδκ x y ) dx dy − ∫∫ ( pw)dx dy = 0 2.32 A A 2.12 PLACAS ANISOTRÓPICAS Las placas analizadas han sido asumidas de ser homogéneas y de material isotrópico. De manera que, las placas de materiales anisotrópicos tienen importantes aplicaciones debido a su excepcional rigidez. Un material no isotrópico o anisotrópico es quien las propiedades dependen de la dirección. El más simple entre ellos, el ortotrópico, es aquel en el que las propiedades del material difieren en dos direcciones mutuamente perpendiculares. Algunos materiales como la madera pueden ser modelados por propiedades ortotrópicas asumiendo que los efectos locales son ignorados. Ejemplos de estos son cubiertas de metal enrollada y corrugada, filtros en construcciones de placas sándwichs, madera laminada, compuestos de fibras reforzadas, concreto reforzado. Un material que tiene dos o más constituyentes es considerando un compuesto. Usualmente, los componentes consisten de materiales reforzados de alta fuerza (acero, vidrio, fibra) incrustado alrededor de un material (resina, concreto, nylon) llamado matriz. De este modo, un material compuesto tiene un relativo gran radio de fuerza a peso, tan bien como sus deseables características. Una placa de material compuesto u ortotrópico es ampliamente usado en muchas de las aplicaciones espaciales, en construcciones, en contenedores a presión en componentes de ingeniería tal como bombas diafragma y acoplamientos. 2.9.1 Relaciones Básicas. La solución de los problemas de flexión de placas ortotrópicas requiere la reformulación de las leyes de Hooke [5]. σ x = Ex ε x + Ex y ε y σ y = Ey ε y + Ex y ε x 2.33 τ x y = Gγ x y En forma matricial es: σ x   Ex Ex y 0   ε x    σ y  = E E 0   ε  x y y   y  → σ = E ε 2.34 τ x y   0 0 G γ x y  Donde los cuatro módulos Ex , Ey , Ex y y G son todos independientes entre si. E ' 'E = x E E = yx 1−ν x ν y y 1−ν x ν y 2.35 E 'x ν E ' ν Ex y = y = y x 1−ν xν y 1−ν xν y Aquí ν x , ν y y E ' x , E ' y son los radios de Poisson efectivos y los módulos de elasticidad efectivos, respectivamente. Los subíndices “x” y “y” relaciona las direcciones. El modulo cortante de elasticidad G es el mismos para ambos materiales isotrópico y ortotrópico. _________________ [5] ORTIZ BERROCAL, Luis. Elasticidad. 3 ed. Madrid: Mac Graw Hill, 1998. 549 p. Las relaciones esfuerzo-deformación son basados en consideraciones geométricas y estas no cambian en las placas ortotrópicas. ∂2w ∂2w ∂2w ε x = −z 2 ε y = −z γ = −2z 2.36 ∂x ∂y2 x y ∂x∂y Los esfuerzos son obtenidos introduciendo la ecuación 2.36 en 2.33  ∂2w ∂2w  σ x = −z  Ex + E∂x2 x y ∂y 2    ∂2w ∂2w  σ y = −z  Ey 2 + Ex y 2  2.37  ∂y ∂x  ∂2w τ x y = −2Gz ∂x ∂y Para los momentos integramos la ecuación 2.10: M  2 D ∂ w ∂ 2 D w x = − x 2 + x y 2   ∂x ∂y   M D ∂ 2w ∂2w  y = − y + D 2.38  ∂y 2 x y ∂x2  ∂2M 2G wx y = − x y ∂x∂y donde t3E 3 t3E 3 D = x t E D y D x y G t Gx y = x y = x y = 12 12 12 12 Las expresiones para Dx , Dy , Dx y y Gx y representa la rigidez de flexional y la rigidez de torcional de una placa ortotrópica, respectivamente. Podemos obtener las fuerzas cortantes en la placa sustituyendo la ecuación 2.38 en 2.17 y 2.16  Q ∂ D ∂ 2w H ∂ 2w  x = − + ∂x  x 2 2  ∂x ∂y  2.39 ∂  2Q D ∂ w H ∂ 2w  y = − + ∂y  y ∂y 2 ∂x2  donde H = Dx y + 2Gx y La ecuación diferencial de gobierno de la deflexión para una placa ortotrópica es D ∂ 4w 2H ∂ 4w ∂4 x 4 + 2 2 + D w y 4 = p 2.40 ∂x ∂x ∂y ∂x 2.9.2 Determinación de la matriz de elasticidad. El cálculo de la matriz de elasticidad para materiales isotrópicos fue descrita en las secciones anteriores. Sin embargo, para materiales ortotrópicos en los cuales la orientación de los ejes del material no coincide con la de los ejes globales, es necesario calcular una nueva matriz elástica que relacione los ejes del material con los globales. Figura 10. Lamina ortotrópica y y 2 1 z σ y τ x y x τ x y σ x σ x τ x y θ τ x x y σ y τ 1 2 Area de la cara =dA σ 2 τ 1 2 σ 1 τ x y σ θ Area de la cara =dA x σ x τ x y θ τ x y τ x y σ σ y y (a) (b) Del diagrama de cuerpo libre de la figura 10 y sumando fuerzas en la dirección 1: ∑ F1 = 0 σ1dA −σ x (dAcosθ )cosθ −σ y (dAsenθ )senθ 2.41 −τ x y (dAcosθ )senθ −τ x y (dAsenθ )cosθ = 0 Del diagrama de cuerpo libre de la figura 10b y sumando la fuerzas en la dirección 2: ∑ F2 = 0 σ 2dA −σ x (dAsenθ )senθ −σ y (dAcosθ )cosθ 2.42 +τ x y (dAcosθ )senθ +τ x y (dAsenθ )cosθ = 0 Del diagrama de cuerpo libre de la figura 10b y sumando fuerzas en la dirección 1: ∑ F2 = 0 τ12dA +σ x (dAsenθ )cosθ −σ y (dAcosθ )senθ 2.43 −τ x y (dAcosθ )cosθ +τ x y (dAsenθ ) senθ = 0 Simplificando las ecuaciones 2.41, 2.42, 2.43 σ1 = σ x cos 2 θ +σ ysen 2θ + 2τ x y senθ cosθ , σ 2 = σ sen 2θ +σ cos2x y θ + 2τ x y senθ cosθ , 2.44 τ12 = σ xsenθ cosθ +σ ysenθ cosθ +τ 2 2 x y (cos θ − sen θ ) La ecuación 2.44 puede escribirse en forma matricial como: σ1   cos 2 θ sen2θ 2senθ cosθ  σ x   σ  2  =   sen 2θ cos2 θ −2senθ cosθ   σ  y  2.45 τ   2 2    12  −senθ cosθ senθ cosθ cos θ − sen θ  τ x y  La matriz de 3x3 en la ecuación 2.45 es llamada matriz de transformación y se denota por [T]. Esta también es usada para transformar deformaciones. Si se desea transformar del sistema de coordenadas 1-2 al sistema de coordenadas x-y, puede utilizarse la inversa de [T], y esta dada por:  cos2 θ sen2θ −2senθ cosθ  T −1 =  sen 2θ cos2 θ 2senθ cosθ  2.46 senθ cosθ −senθ cosθ cos 2 θ − sen2θ  por tanto, σ1  σ x  σ x  σ1  σ  = T σ  , σ  = T-1 σ  2   y   y   2  2.47 τ  τ  τ  τ 12 x y x y 12  Similarmente para deformaciones,  ε1   ε x   ε x   ε1   ε  2  = T  ε  ,  ε  = T-1 ε  y   y   2  2.48 γ       12  ε x y  γ x y  ε12  Sustituyendo en la ecuación 2.34 en la segunda de la ecuación 2.47: σ x   ε1  1 0 0  ε1  σ  = T-1 E  ε  = T-1 y   2  E  0 1 0    ε  2  2.49 τ  γ  0 0 2 ε x y 12 12  Ahora poniendo la primera parte de la ecuación 2.48 en la ecuación 2.49: σ x  1 0 0  ε x    -1     σ y  = T E 0 1 0 T  ε y  2.50 τ x y  0 0 2 ε x y  Por tanto la nueva matriz de elasticidad es: σ x  1 0 0  σ  y  = T -1 E 0 1 0   T 2.51 τ x y  0 0 2 las componentes de la matriz de elasticidad (E*) son: E* 4 2 2 411 = E11m + 2 (E12 + 2E66 ) m n + E22n E* = ( E + E − 4E 2 2 4 412 11 22 66 ) m n + E12 (m + n ) E* 4 2 2 422 = E11n + 2 (E12 + 2E66 ) m n + E22m E* 2.52 16 = ( E11 − E 312 − 2E66 ) m n + (E 312 − E22 + 2E66 ) mn E*26 = (E11 − E 3 312 − 2E66 ) n m + (E12 − E22 + 2E66 ) m n E* 2 2 4 466 = ( E11 + E22 − 2E12 − 2E66 )m n + E66 (m + n ) donde m = cosθ , n = senθ . Sustituyendo la ecuación 2.52 en 2.50, se tiene: σ  E* * *x 11 E12 E16   ε x  σ  = E* E* E*   ε  y   12 22 26   y  τ  E * E* 2E*  x y 16 26 66 ε x y  2.53 σ  E* E* E*x 11 12 16   ε x   σ  y  =  E * * *    12 E22 E26   ε y  τ  E * * *   E E  x y 16 26 66 γ x y  2.9.3 Placa de multicapas isotrópicas. Figura 11. Placa laminada simétricamente construida p 2n + 1 x Plano medio 2t1 = t n + 1 k Ek, vy k 2 1 2t 2t n+1 k z Una estructura compuesta de un número arbitrario de capas definidas pueden ser frecuentemente aproximadas a una placa laminada. Generalmente, el ensamble de estas estructura, cada capa posee diferente espesor, orientación de los ejes principales y propiedades anisotrópicas. En este momento se analizaran placas de capas isotrópicas. Placas con multicapas. Son fabricadas de tal manera que ellas actúan como un solo material. Las capas no pueden deslizarse sobre otra, y los desplazamientos permanecen continuos; las deformaciones varían linealmente sobre el espesor. Tal que, las relaciones de esfuerzo-deformación para una capa k, son expresadas así: ∂2w 2 2 ε (k ) = −z ε (k ) z ∂ w= − γ (k ) ∂ wx k 2 x k 2 x y = −2zk 2.54 ∂x ∂y ∂x∂y Las leyes de Hooke son: σ (k ) E = k ε (k ) +ν ε (k )x 1−ν 2  x k y  k σ (k ) E = k ε (k ) +ν ε (k )y 1−ν 2  y k x  2.55 k E τ (k )x y = k γ (k ) 2 (1+ν k ) x y Sustituyendo las deformaciones de la ecuación 2.54 en 2.55, integrando sobre cada capa y sumando los resultados, obtenemos los esfuerzos resultantes actuando en el plano medio.  M x  σ  x  M  z  y  = ∑ ∫ σ  y  z dz  z−1 M  k   x y  τ x y (k ) Los esfuerzos para la capa k son: z Ek  ∂ 2w ∂2w  σ x = − k +ν1−ν 2  2 k 2 k  ∂x ∂y  E  ∂2w ∂2z k w  σ x = − k 2  +ν 1−ν k  ∂y 2 k ∂x2  E ∂2w τ kx y = −zk 1+ν k ∂x ∂y El método general de derivar las ecuaciones de gobierno para placas multicapas sigue un patrón idéntico con el que fue descrito en la sección 2.4. Ahora la ecuación 2.21 es: ∇4w p= D 2.56 t donde Dt es la rigidez flexional transformada de placas laminadas Una placa de capas de simetría estructural cerca del plano medio son de importancia. Para una placa de 2n +1de capas isotrópicas simétricas (figura 11) la rigidez flexional transformada esta dada por 2  n E E t3D k 3 3 n+1 n+1  t = t − t3 ∑1 2 (−ν k t+1 ) + 2  2.57  k =1 k 1−ν n+1  2.9.4 Placa compuesta laminada Figura 12. Laminado compuesto z1 z2 z3 Plano medio z x t z 4 z z 5 7 6 En una lámina se presentan fuerzas y momentos (en el plano y en flexión) estos pueden expresarse integrando los esfuerzos a través del espesor h. N  h σ  σ    x  2  x   x   N = N y  = ∫ σ y  + σ y   dz 2.58 N  −h       2xy   τ xy  τp  xy b  M  h σ      x  2  x σ   x   M = M y  = ∫ z σ y  + σ y   dz 2.59    M  −h     xy  2  τ τ xy p  xy b  donde los subíndices p y b indican los estados en el plano y en flexión respectivamente. Las integrales en las ecuaciones 2.58 y 2.59 se realizan sobre cada capa de lámina para luego ser sumadas. Usando el esquema de laminado de la figura 12 las ecuaciones 2.58 y 2.59 pueden escribirse de la siguiente manera: N  σ  σ    x  n zk+1  x   x   N = N y  = ∑ ∫ σ y  + σ y   dz 2.60   k=1 z     N k xy  τ xy  τp xy b k M  σ      x  n zk+1  x σ   x   M = M y  = ∑ ∫ z σ y  + σ y   dz 2.61 M  k =1 z k τ  τ    xy   xy p  xy b k Simplificando las ecuaciones 2.60 y 2.61 y sustituyendo el vector de esfuerzo (σ = E ε ), siendo E la matriz elásticas. ε  ε   n zk+1  x   x   n zk +1 N = ∑E*k ∫ ε y  + ε y   dz = ∑ E*k ∫ (ε p + εb ) dz k=1 z γ     k =1 z  x y p  γ x y b       n z ε ε k +1  x   x   n zk +1 M = ∑ E*k ∫ z ε *y  + ε y   dz = ∑ Ek  ∫ z (ε p + εb )dz k =1 z     k =1 z  γ γ x y p  x y b  sustituyendo la ecuación 2.5, se obtiene: n zk+1 N = ∑E*k ∫ (ε + zκ ) dzk k=1 zk z n k +1 M = ∑ E*k ∫ z (ε + zκ )k dz k =1 zk Resolviendo las integrales: n N = ∑E* 1 2 2k {( zk+1 − zk )ε + 2 ( zk+1 − zk )κ} k=1 k n M = ∑ E* {1 ( z2 − z2 )ε + 1 ( z3 − z3k 2 k +1 k 3 k +1 k )κ} k k=1 los resultados de la integrar pueden escribirse como: n n n A = ∑E* ( z − z ) , B = ∑ E* 1 ( z2 2 * 1 3 3k k +1 k k 2 k+1 − zk ), C = ∑Ek 3 ( zk +1 − zk ) 2.62 k =1 k =1 k=1 y finalmente se tiene:  Nx   A11 A12 A16 B11 B12 B16   ε x   N      y  A12 A22 A26 B12 B22 B26   ε y   Nx y  A16 A26 A66 B16 B26 B66  γ x y    =     M B B B 2.63  x   11 12 16 C11 C12 C16   κ x   M  y B12 B22 B26 C C      12 22 C 26 κ  y  M x y  B16 B26 B66 C16 C 26 C 66  κ x y  o  N  A B ε   =M      2.64   B C κ  3. FORMULACION POR ELEMENTOS FINITOS 3.1 ANALISIS DE PLACA EN EL PLANO Se hará el análisis en el plano sin tener en cuenta las fuerzas inerciales y de amortiguamiento. Discretización del funcional de energía [6] mediante M.E.F. Π = ∫σ Tε dV − ∫ f Tu dV − ∫ STΩε dΩ 3.1 V V Ω En el caso bidimensional (para condiciones de frontera homogénea) las componentes de desplazamiento pueden ser aproximadas por: 4 u ≅ ∑Niui = N1u1 + N2u2 + N3u3 + N4u4 i=1 3.2 4 v ≅ ∑Nivi = N1v1 + N2v2 + N3v3 + N4v4 i=1 Las cuales pueden ser escritas en forma vectorial como: u1    v1  u2    N 0 N 0 N 0 N 0  v U ≅ 1 2 3 4  2  → U ≅ N∆ 3.3  0 N1 0 N2 0 N 0 N  3 4  u  3  v3  u   4  v4  _________________ [6] REDDY, J.N. An Introduction to the Finite Element Method. 2 ed. Mc. Graw Hill, 1993, 672 p. Para el caso de elasticidad bidimensional se tiene las siguientes relaciones:  ε x  ∂ ∂x 0     u  ε y  =  0 ∂ ∂y     → ε = TU 3.4 2γ  v x y  ∂ ∂y ∂ ∂x   σ = Eε → σ = ETU 3.5 reemplazando δ Π = ∫ ((ETU )T δ (EU ) − f TδU )dV − ∫ STΩδU dΩ V Ω Teniendo en cuenta u ≅ N∆ δ Π = ∫ ((ETN∆)T δ (TN∆) − f Tδ (N∆))dV − ∫ STΩδ (N∆)dΩ V Ω haciendo B = TN y tomando factor común δ∆ 0 = ∫ ((EB∆)T B − f T N )δ∆dV − ∫ STΩNδ∆dΩ V Ω 0 = ∫ (∆T BT E B − f T N )δ∆dV − ∫STΩNδ∆dΩ V Ω lo cual debe cumplirse para todo δ∆ entonces: 0 = ∫ (∆T BT E B − f T N )dV − ∫ STΩNdΩ V Ω 0 = ∫∆T BT EBTdV − ∫ f T NdV −∫STΓ NdΩ V V Ω  T t  t t ∫ B E BdΩ∫ dz ∆ =−t  ∫ N T fdΩ∫ dz + ST NdΓ dz   −t ∫ Γ ∫−t Ω Ω Γ   t∫BT E BdΩ ∆ = t∫ NT fdΩ +t∫STΓ NdΓ  Ω  Ω Γ Lo cual puede ser re-escrito como: K∆ = f 3.6 K = t∫ BT E BdΩ 3.7 Ω f = t∫ N T fdΩ + t∫ STΓ Nδ∆dΓ 3.8 Ω Γ donde K es de 8 x 8 y f es de 8 x 1 3.2 ANALISIS DE PLACAS A FLEXION Se tendrá en cuenta las mismas consideraciones que en el análisis de membrana: Discretización del funcional de energía mediante F.E.M. Π = ∫σ Tε dV − ∫ f Tu dV − ∫ STΩε dΩ V V Ω Las componentes de desplazamiento pueden ser aproximadas por: n w ≅ ∑φi wi = φ1 w1 + φ2 w2 +K+φn wn → w = φ ∆ 3.9 i=1 donde φi son las funciones de interpolación de Hermite. Las cuales pueden ser escritas en forma vectorial como:  w1     ∂w1   ∂x   ∂w   1   ∂y   2   ∂ w1  ∂x∂y  w ≅ [φ  11 φ21 φ31 φ41 L φ14 φ24 φ34 φ44 ] M   w4     ∂w4   ∂x   ∂w   4   ∂y   ∂2  w4  ∂x∂y  Para el caso de láminas delgadas se tiene las siguientes relaciones:   ∂w      ∂   ε   ∂x   ∂x x    ∂w  z z ∂  ε =  ε y  = −   = −  [w] → ε = −zT∆ 3.10  ∂y2γ     ∂y   x y   ∂2w   ∂2      ∂x∂y  ∂x∂y  σ = Eε → σ = −z ET ∆ 3.11 reemplazando δΠ = ∫ ((−z ETw)T δ (−zTw) − f Tδ w)dV − ∫STΩδ wdΩ V Ω Teniendo en cuenta w ≅ N∆ δΠ = ∫ ((−z ETφ∆)T δ (−zTφ∆) − f Tδ (φ∆))dV − ∫ STΩδ (φ∆)dΩ V Ω haciendo B = TN y tomando factor común δ∆ .. 0 = ∫ z2 (E B∆)T B − ρ(φ ∆)T φ − f Tφ  Tδ∆dV − S φδ∆dΩ V   ∫ Ω Ω 0 = ∫ ( z2∆T BT E B − f Tφ )δ∆dV − ∫ STΩφδ∆dΩ V Ω lo cual debe cumplirse para todo δ∆ entonces: 0 = ∫ (∆T BT E B − f Tφ )dV − ∫ STΩφdΩ V Ω 0 = ∫∆T BT DBdV − ∫ f TφdV −∫ STΩφdΩ V V Ω   ∫ BT E BdV  ∆ = ∫φT fdV +∫STΩφdΩ V  V Ω  T t 2 2  T t 2 t 2 ∫ B E BdΩ∫ z dz ∆ = ∫φ fdΩ∫ dz +∫STΓφdΓ∫ dz  −t 2  −t 2 −t 2Ω Ω Γ  t3 ∫ BT  T  E BdΩ ∆ = t∫φ fdΩ +t12 ∫S T ΓφdΓ  Ω  Ω Γ Lo cual puede ser re-escrito como: K∆ = f K t 3 = ∫ BT E BdΩ 3.12 12 Ω f = t∫φT fdΩ + t∫ STΓφdΓ 3.13 Ω Γ donde K es de 16 x 16 y f es de 16 x 1 3.5 MATRIZ COMBINADA Considérese una placa sometida simultáneamente a acciones en su plano y a flexión. Cada una de estas acciones (en su plano y a flexión) fueron descritas en las secciones 3.1 y 3.2 respectivamente. Los efectos en el plano dieron como rebultado una matriz de rigidez ( K ep ) de 8x8, un vector de fuerza ( f e ) de 8x1 y un vector de desplazamiento (∆) de 8x1.  f f x  u =  ∆ =f    3.14  y  v Lo anterior representan las fuerzas y desplazamientos en x y en y para un nodo. Los efectos a flexión dieron como rebultado una matriz de rigidez ( K eb ) de 16x16, un vector de fuerza ( f e ) de 16x1 y un vector de desplazamiento (∆) de 16x1.  f z   w   M    f =  x  θ ∆ =  x   M y  θ y  3.15 M    x y  θ x y  Lo anterior representan las fuerzas y desplazamientos en z, y momentos y giros en “x”, ”y” y “xy” para un nodo. La matriz de rigidez y el vector de fuerzas utilizados para el análisis de placas sometidas a cargas combinadas (en el plano y a flexión) son expresados:  0 L 0  K e  p M M  e  0 L 0K = 0 L 0  3.16   M M K e   b  0 L 0  T f =  f x1 f y1 L f x4 f y 4 M f z1 M x1 L M y 4 M z 4  3.17 3.4 ELEMENTOS FINITOS La elección las funciones de forma está condicionada no sólo por la geometría de los dominios locales o elementos finitos, sino también por el tipo de problema de campo que se intenta resolver (elasticidad, plasticidad, transferencia de calor, electromagnetismo, etc) y por la manera de abordarlo (número de campos independientes). Dentro de cada dominio elemental o elemento finito, el campo independiente se aproxima o interpola mediante funciones de forma que cumplen las condiciones siguientes: 1. Cada función de forma está asociada a un punto dentro del elemento, de manera que vale 1 en ese punto y se anula en los puntos asociados a las otras funciones de forma. Estos puntos especiales dentro de cada elemento reciben el nombre de nodos. 2. Las funciones de forma correspondientes a cada elemento están definidas de manera que hay continuidad del campo independiente al pasar de un elemento a otro, es decir, al agregar las aproximaciones locales para construir la aproximación extendida al dominio completo se obtiene una función continua. 3. Las funciones de forma son capaces de representar de manera exacta campos con un mínimo orden de variación, el cual depende del tipo de problema. 4. Las funciones de forma dentro de cada elemento son suficientemente lisas o “suaves” como para permitir el cálculo de las derivadas que sean necesarias para el planteamiento de la forma débil del problema. Es muy ventajoso trabajar siempre, independientemente del elemento concreto, con el mismo dominio de interpolación y con las mismas funciones de forma. De esta manera se simplificaría la evaluación de estas funciones y de sus derivadas y se conseguiría realizar las integrales siempre sobre el mismo dominio de integración, operación programable mucho más fácilmente. Existen dominio en que los elementos finitos como cuadrado o triángulos no son capaces de enmallarlo, por tanto se necesitan elemento que puedan amoldarse a casi cualquier geometría. Esto es posible por medio de la transformación isoparamétrica, la cual es aplicable tanto para triángulos como para cuadriláteros. Para este estudio se transformaran elementos cuadrados de cuatro nodos, ya que estos resultan ser más precisos que los triángulos de tres nodos. 3.4.1 Transformación isoparamétrica. Se puede llevarse a cabo por medio de las siguientes operaciones: 1. Se proyecta el dominio elemental o elemento finito en un dominio estándar Figura 13. Elemento Transformado y η (-1,1) (1,1) 3 4 3 4 PROYECCION x ξ 1 2 1 2 (-1,-1) (1,-1) 2. La proyección puede interpretarse también como un cambio de variables dentro de cada dominio local: x = x(ξ ,η) Cambio de Variable  y = y(ξ ,η) es decir, x = N1x1 + N2x2 + N3x3 + N4x4 3.18 y = N1y1 + N2 y2 + N3 y3 + N4 y4 3. Las funciones de interpolación o funciones de forma se definen en coordenadas locales: N 11 = (1 )(1 1 −ξ −η ) N3 = (1+ ξ )(1+η )4 4 3.19 N 12 = (1+ξ )(1−η ) N 1 4 = (1−ξ )(1+η )4 4 o N 1i = (1+ ξ ξi )(1+ηηi ) 3.20 4 donde ξi y ηi son las coordenadas del nodo i. 4. Se expresa las derivadas de la función Ni(x,y) en términos de las coordenadas locales ξ,η. Usando la regla de la cadena de diferenciación, se obtiene ∂Ni ∂Ni ∂x ∂N= + i ∂y ∂ξ ∂x ∂ξ ∂y ∂ξ ∂Ni ∂N= i ∂x ∂N ∂y+ i ∂η ∂x ∂η ∂y ∂η o bien ∂Ni  ∂Ni   ∂ξ     = J  ∂x  3.21 ∂N   i  ∂N  i   ∂η   ∂y  donde J es la matriz Jacobiana  ∂x ∂y   ∂ξ ∂ξ  J =   3.22  ∂x ∂y  ∂η ∂η  La ecuación 3.21 puede invertirse, ∂Ni  ∂Ni   ∂x   ∂ξ   N  = J −1     3.23 ∂ i  ∂Ni   ∂y   ∂η  Esto requiere que la matriz jacobiana sea no-singular. Estas expresiones se usarán para la obtención de la matriz de rigidez del elemento. El diferencial de área esta dado por dΩ = dx dy = detJ dξ dη 3.24 Para así obtener K e = t ∫ BTe D B det J dξ dη 3.25 Ω Es importante tener en cuenta, que si el elemento está muy distorsionado (tiene ángulos muy distintos del recto), la interpolación de Ni(x,y) en el espacio real puede ser muy distinta de la que tiene lugar en coordenadas locales ξ,η. Esto quiere decir que, las funciones de forma reales pueden llegar a ser incluso cualitativamente muy diferentes de los polinomios definidos en coordenadas locales. Una interpolación polinómica lineal en el dominio estándar no tiene por qué ser también lineal en el espacio real. Es posible que ni siquiera sea completamente polinómica. La distorsión de los elementos da lugar así en la práctica a una degradación de la calidad de la solución aproximada. 3.4.2 Elementos conformes. Para este análisis la formulación por elementos finitos requiere la continuidad de w (ecuación 3.9) y de su primera derivada son llamadas continuidad C1 . Un elemento rectangular con cuatro nodos, con (w, ∂w ∂x, ∂w ∂y) en cada nodo, requiere 12 términos (n=12) en el polinomio de interpolación de w para obtener la expresión: w = a1 + a2x + a3 y + a4xy + a5x 2 + a6 y 2 + a x27 y + a8xy 2 + a9x 3 + a10 y 3 + a11x 3 y + a xy312 3.26 Esta expresión no es un polinomio completo de cuarto orden; es un polinomio de tercer orden. Para un elemento rectangular de cuatro nodos (n=12), no es posible conseguir la continuidad a través de las fronteras de los elementos de w y ∂w ∂η , si se utilizan tres grados de libertad por nodo. Los elementos que violen cualquier condición de continuidad son conocidos como elementos no-conformes [7]. De esa manera el elemento rectangular de cuatro nodos con un polinomio de 12 términos es un elemento no-conforme. ________________ [7] REDDY, J.N. An Introduction to the Finite Element Method. 2 ed. Mc. Graw Hill, 1993, 672 p. Un elemento de cuatro nodos con cuatro grados de libertad (w, ∂w ∂x, ∂w ∂y, ∂2w ∂x ∂y) requiere un polinomio completo de cuarto orden. Y a este se le conoce como un elemento es conforme. Las funciones de interpolación que cumplen con la anterior condición son las funciones de Hermite, y son: Tabla 1. Funciones de interpolación de Hermite φ [8] w 1 (ξ + ξ )2i (ξ ξi − 2)(η +ηi ) 2 (ηη − 2) 16 i ∂w / ∂ξ 1− ξ (ξ + ξ )2 (ξ ξ −1)(η +η )2i i i i (ηηi − 2) 16 Para un nodo i 1 (i =1, 2,3, 4)∂w / ∂η − (ξ +ξ 2i ) (ξ ξi − 2)ηi (η +η ) 2 i (ηη16 i −1) 1 ∂2w / ∂ξ ∂η ξi (ξ +ξ ) 2 i (ξ ξ −1)η (η +η 2 i i i ) (ηηi −1) 16 Para el cálculo de B en la matriz de rigidez en placas con elementos de cuatro grados de libertad se calcula de la siguiente manera: _____________ [8] REDDY, J.N. An Introduction to the Finite Element Method. 2 ed. Mc. Graw Hill, 1993, 672 p. 2 2 −1 ∂2  φ e   ∂xe   ∂ye   ∂x   e ∂ye  i    2  ∂x2   ∂ξ     ∂ξ      ∂ξ ∂ξ     2 e  2 2 ∂ φ i    ∂xe   ∂ye   ∂xe ∂y  x2  = e  ∂   ∂η   ∂η  2          ∂η ∂η    2 e  ∂ φ  i   ∂xe   ∂xe   ∂ye  ∂ye   ∂xe  ∂ye   ∂xe  ∂ye  ∂x ∂y           +     ∂ξ   ∂η   ∂ξ  ∂η   ∂η  ∂ξ   ∂ξ  ∂η  3.27   ∂2φ e   ∂2i x 2 e ∂ y e    ∂ξ 2        ∂ξ 2 ∂ξ 2  ∂φ e i    ∂2φ e× i    ∂2 x ∂2 y   ∂x   2  −  e e       ∂η 2 2 e   ∂η ∂η  ∂φ i    ∂2φ e   ∂2 x ∂2 y   ∂y i e e         ∂ξ ∂η  ∂ξ ∂η ∂ξ ∂η    donde φi son las funciones de forma de Hermite ∂φ e  ∂φ ei i  ∂x   [ ]−1  ∂ξ     = J   3.28 ∂φ e e i  ∂φ i   ∂y   ∂η  finalmente la matriz B es: ∂2φ e 2 e 21 ∂ φ 2 ∂ φ e i   ∂x2 L ∂x2 ∂x2   ∂2φ e ∂2φ e ∂2φ eB = 1 2 L i  ∂x2 ∂x2 ∂x2  3.29   ∂2φ e ∂2φ e ∂21 2 φ e i   L ∂x ∂y ∂x ∂y ∂x ∂y  3.5 TRANSFORMACIÓN DE COORDENADAS GLOBALES La matriz de rigidez deducida en la sección anterior utiliza un sistema de coordenadas locales, ya que los efectos de flexión y de las acciones “en el plano” se obtenían originalmente en dicho sistema. Para ensamblar los elementos y para escribir las ecuaciones de equilibrio apropiadas será necesario un transformación de coordenadas locales (designadas por x’ y’ z’) a un sistema global común (coordenadas x, y, z). Figura 14. Coordenadas locales y globales Y´ Z Z´ X´ Y X´X X Además será mas conveniente al principio definir los nodos del elemento por sus coordenadas globales y a partir de ellas calcular las coordenadas locales, para lo cual necesitamos una transformación inversa. La figura 14 se muestra los dos sistemas coordenados. Las fuerzas y los desplazamientos de un nodo se transforman de sistema global a sistema local por una matriz L de tal forma que ∆ 'i = L ∆i f 'i = L fi 3.30 donde λ 0L =   3.31 0 λ y λ es una matriz de dimensiones 3 X 3 de los cosenos directores de los ángulos que forman entre si los dos sistemas de ejes, o sea λx ' x λx ' y λx ' z  λ = λ  y ' x λy ' y λy ' z  3.32 λz 'x λz ' y λ z ' z  donde λx ' x es el coseno del ángulo formado por los ejes “x” y “x’”, etc. Por tanto, para el conjunto de fuerzas que actúan sobre los nodos de un elemento podemos escribir [∆ ']e = T [ ]e∆ 3.33 por ende la matriz de rigidez de un elemento en coordenadas globales se obtiene por K e = T T [K ']e T 3.34 en las dos ecuaciones anteriores T viene dado por L 0 0 L 0 L 0  T =   3.35 0 0 L   M   que es una matriz diagonal formada por un número de matrices L igual al numero de nodos del elemento. 3.5.1 Elementos orientados en forma arbitraria Figura 15. Coordenadas locales y globales para un elemento cuadrilátero transformado z n m y’ z’ j y i x’ x Una de las direcciones de los ejes locales puede elegirse arbitrariamente, se tomara la del eje x’ orientada en la dirección del lado ij del elemento, tal como se muestra en la figura 14. El lado esta definido por el vector Vij y en función de las coordenadas locales que se tienen.  x j − xi  V  i j =  y j − yi  3.36  z j − zi  los cosenos directores se obtienen dividiendo las componentes de este vector por su longitud, o sea, definiendo un vector unitario λx ' x  x ji  v  1x ' = λ    x ' y  = yl  ji  3.37  ijλ   z x ' z ji  donde lij = x 2 ji + y 2 ji + z 2 ji y xij = x j − xi , etc. Es preciso establecer la dirección de z’, que debe ser normal al plano del elemento. Según las propiedades del producto vectorial de los dos vectores, puede obtenerse dicha dirección a partir de un vector que sea el producto vectorial de los dos lados del elemento. Vz ' = Vij ×Vim 3.38 representa un vector normal al plano del elemento. Los cosenos directores del eje z’ se identifican sencillamente con los cosenos directores Vz’, obteniéndose un vector unitario λz 'x  v = λ z '  z ' y  3.39 λz ' z  finalmente, los cosenos directores del eje y’ se establecen de una manera similar a partir de un vector normal a las direcciones x’ y z’ λy ' x  v  y ' = λy ' y  3.40 λ y ' z  3.6 INTEGRACIÓN NUMÉRICA Considere el problema de evaluar numéricamente una integral unidimensional de la forma 1 I = ∫ f (ξ ) dξ 3.41 −1 El método de cuadratura Gaussiana [9] para evaluar I se da a continuación. Este ________________ [9] CHANDRUPATLA, Tirupathi. Introducción al Estudio de Elemento Finito en Ingeniería. 2 ed. Mexico: Pearson,1999, 462p método ha demostrado ser muy útil en el trabajo con elementos finitos. La extensión a integrales en dos y tres dimensiones puede obtenerse fácilmente. Consideremos la aproximación de n puntos 1 I = ∫ f (ξ ) dξ ≈ w−1 1 f (ξ1 ) + w2 f (ξ2 ) + L + wn f (ξn ) 3.42 donde w1, w2 ,K wn son los pesos y ξ1, ξ2 ,Kξn son los puntos de muestreo o puntos de Gauss. La idea detrás e la cuadratura gaussiana es seleccionar los n puntos de Gauss y los n pesos de manera que la ecuación 3.42 proporcione una respuesta para polinomios f (ξ ) de un grado tan grande como se quiera. En otras palabras, la idea es que si la formula de integración de los n puntos es exacta para todos los polinomios hasta un grado tan alto como sea posible, entonces la formula funcionará bien aun si f no es un polinomio. Tabla 2. Puntos de Gauss y Pesos para cuadratura Gaussiana [10] 1 n ∫ f (ξ ) dξ ≈ w−1 ∑ i f (ξi ) i=1 Número de puntos, n Localización ξi Pesos, wi 1 0.0 2.0 2 ±1 3 = ± 0.5773502692 1.0 3 ± 0.7745966692 0.5555555556 0.0 0.8888888889 4 ± 0.8611363116 0.3478548451 ± 0.3399810436 0.6521451549 ± 0.9061798459 0.2369268851 5 ± 0.5384693101 0.4786286705 0.0 0.5688888889 ± 0.9324695142 0.1713244924 6 ± 0.6612093865 0.3607615730 ± 0.2386191861 0.4679139346 _________________ [10] CHANDRUPATLA, Tirupathi. Introducción al Estudio de Elemento Finito en Ingeniería. 2 ed. Mexico:Pearson,1999, 462 p. 3.6.1 Integrales bidimensionales La extensión de una cuadratura gaussiana a integrales bidimensionales de la forma 1 1 I = ∫ ∫ f (ξ ,η ) dξ dη 3.43 −1 −1 es consecuencia natural ya que 1 n I  ≈ ∫−1 ∑wi f (ξi ,ηi ) dη  i=1  n n I ≈ ∑w j ∑wi f (ξi ,  ηi ) j=1  i=1  n n I ≈ ∑ ∑wi w j f (ξi ,ηi ) 3.44 i=1 j=1 Para ilustrar el uso de la ecuación 3.44 considere la rigidez del elemento para un elemento cuadrilátero e 1 1k = t∫ ∫ BT D B det J dξ dη 3.45 −1 −1 donde B y det J son funciones de ξ y η . Sea φ el elemento ij-ésimo en el integrado. Es decir, sea φ (ξ ,η ) = t (BT D B det J ) ij Entonces, si usamos una regla de 2 x 2, obtenemos k ≈ w2φ (ξ ,η 2i j 1 1 1 ) + w1w2φ (ξ1,η2 ) + w2w1φ (ξ2 ,η1 ) + w2φ (ξ2 ,η2 ) 3.46 donde w1 = w2 = 1.0 , ξ1 =η1 = −0.5773502692 , y ξ2 =η2 = 0.5773502692 . Los puntos de Gauss para la regla de dos puntos usada anteriormente se muestran en la Figura 16. Cuadratura Gaussiana en dos dimensiones η 1 η2 = 3 4 3 ξ 1 1 2 η1 = − 3 1 ξ1 = − 1 ξ 3 2 = 3 1 1 ∫ ∫ f (ξ ,η ) ≈ w21φ (ξ1,η1 ) + w1w2φ (ξ ,η ) + w wφ (ξ ,η ) + w2φ (ξ ,η−1 −1 1 2 2 1 2 1 2 2 2 ) 3.7 AJUSTE POR MINIMOS CUADRADOS PARA UN CUADRILATERO DE CUATRO NODOS Figura 17. Ajuste por mínimos cuadrados para un cuadrilátero η 4 3 q4 q3 4 3 s4 s3 ξ 1 2 s s q1 1 2 q2 1 2 Si q = [ q1, q2, q3, q4 ]T es el vector de valores nodales del elemento que son determinados por ajuste de mínimos cuadrados, entonces se define usando el error en cuatro puntos interiores. Si s = [ s1, s2, s3, s4 ]T representa el vector de valores interpolados en los cuatro puntos interiores y a = [ a1, a2, a T3, a4 ] representa los valores reales de las variables (figura 17), el error puede ser definido como ε = ∑(s − a)T (s − a) e 3.47 ε = ∑(sT s − 2sT a − aT a) e Los cuatro puntos interiores se toman generalmente como los cuatros puntos de integración de Gauss. Los valores del esfuerzo concuerdan bien en esos puntos. Si N denota  N1 N1 1 11 2 N3 N4   2 2 2 2  N =  N1 N2 N3 N4  3 3 3 3  3.48 N   1 N2 N3 N4 N 4 N 4 4 4   1 2 N3 N4  donde N ij es la función de forma N j evaluada en el punto interior I, entonces s puede escribirse como s = Nq 3.49 Insertando esto en la ecuación 3.47, el error ε toma la forma ε = ∑qT N T Nq − 2qT NT a + aT a e notando que NTN es similar a la rigidez del elemento K e y que NT a lo es al vector de fuerza del elemento, puede efectuarse el ensamble de la rigidez y del vector fuerza. Las ecuaciones de la matriz ensamblada ponerse como KQ = F 3.50 La solución de este conjunto de ecuaciones da Q, que es el vector de los valores nodales de la variable considerada para el ajuste por mínimos cuadrados de los valores del elemento. 3.8 ANÁLISIS DINÁMICO 3.8.1 Comportamiento Dinámico de las estructuras elásticas con amortiguamiento lineal. Cuando los desplazamientos de un cuerpo elástico varían en función del tiempo, entran en juego dos tipos de fuerzas adicionales. Las primeras son las fuerzas de inercia, que para una aceleración definida por ∆&& se pueden reemplazar por sus equivalentes estáticas −ρ ∆&& Las componentes de esta fuerza se definen en las mismas direcciones que las del desplazamiento ∆; el coeficiente ρ es simplemente la densidad. La segunda fuerza es la producida por las resistencias (rozamientos) opuestas al movimiento. Estas últimas pueden ser debidas a movimientos microscópicos, a la resistencia del aire, etc., estando en general relacionadas no linealmente con la velocidad ∆& . Sin embargo para simplificar los cálculos siguientes sólo consideraremos las resistencias lineales de tipo viscoso, que en un problema estático equivalente dan lugar a fuerzas por unidad de volumen de intensidad. −µ ∆& En la expresión anterior, µ es una propiedad que se supone puede ser evaluada numéricamente. Las fuerzas (nodales) aplicadas a un elemento están definidas por: f e = − ∫ N Tb dV = − ∫ NTb dV + ∫ NT ρ ∆&& dV + ∫ NT µ ∆& 3.51 V e V e V e V e Donde la primera fuerza es precisamente la debida a las cargas repartidas exteriores, no siendo por tanto necesarias mayores consideraciones acerca del particular. Sustituyendo la ecuación 3.51 en las ecuaciones generales de equilibrio se obtiene finalmente, tras el ensamble, la ecuación diferencial matricial [11] siguiente: M ∆&& + C ∆& + K∆ + f = 0 3.52 En la cual K y f son las matrices globales de rigidez y de fuerzas. Las nuevas matrices C y M son las matrices globales de masa y amortiguamiento y se ensamblan a partir de las submatrices elementales dadas por: C e = ∫ NT µ N dV 3.53 V e M e = ∫ N T ρ N dV 3.54 V e La determinación de la matriz de amortiguamiento C es difícil en la práctica, ya que se desconoce la viscosidad µ. A menudo se hace la hipótesis de que la matriz de amortiguamiento sea una combinación lineal de las matrices de rigidez y de masa, es decir C = α M + β C 3.55 donde α y β se determinan experimentalmente. Esta forma de expresar el amortiguamiento se conoce como “amortiguamiento de Rayleigh” 3.8.2 Método de integración del tiempo de Newmark Tenemos la siguiente ecuación: M ∆&& + C ∆& + K∆ = f 3.56 ___________ [11] REDDY, J.N. An Introduction to the Finite Element Method. 2 ed. Mc. Graw Hill, 1993, 672 p. donde [ M ] = Matriz de Masa [C ] = Matriz de Amortiguamiento [ K ] = Matriz de Rigidez { ∆&& } = Vector de Aceleración { ∆& } = Vector de Velocidad { ∆ } = Vector de Desplazamiento El método de Newmark [12] usa expansión de diferencia finita en el intervalo ∆t : ∆& & && &&n+1 = ∆n + (1−δ ) ∆n + δ ∆n+1  ∆t 3.57 ∆n+1 = ∆n + ∆& &&n ∆t + (1 2 −α ) ∆n +α ∆&& 2 n+1  ∆t 3.58 donde α , δ = Parámetros de integración de Newmark ∆t = tn+1 − tn { ∆n } = Vector de desplazamiento en el tiempo tn { ∆& n } = Vector de velocidad en el tiempo tn { ∆&&n } = Vector de aceleración en el tiempo tn { ∆n+1 } = Vector de desplazamiento en el tiempo tn+1 { ∆& n+1 } = Vector de velocidad en el tiempo tn+1 { ∆&&n+1 } = Vector de aceleración en el tiempo tn+1 Las ecuaciones de gobierno son evaluadas en un tiempo tn+1 M ∆&&n+1 + C ∆& n+1 + K∆n+1 = f 3.59 _____________ [12] ANSYS Theory Referente. 001099. Ninth Edition. SAS IP, Inc. La solución para los desplazamientos en un tiempo tn+1 se obtiene ∆&&n+1 = a0 [∆ − ∆ ] − a ∆& − a ∆&&n+1 n 2 n 3 n 3.60 ∆& & && &&n+1 = ∆n + a6K∆n + a7 ∆n+1 3.61 donde a 10 = 2 a δ 4 = −1α ∆t α a δ a ∆t  δ1 = 5 =  − 2  α ∆t 2  α  3.62 a 12 = a6 = ∆t (1−δ )α ∆t a 13 = −1 a = δ ∆t2α Notando que { u&&n+1 } en la ecuación 3.60 puede ser sustituida dentro de 3.61, las ecuaciones para { u&&n+1 } y {u&n+1 } pueden ser solo expresados en términos desconocidos { un+1 } . Las ecuaciones para { u&&n+1 } y {u&n+1 } son combinadas con 3.59: (a0 M + a1 C + K ) ∆n+1 = f + 3.63 M (a ∆ + a ∆& + a ∆&&0 n 2 n 3 n ) + C (a ∆ + a ∆& + a ∆&&1 n 4 n 5 n ) Una vez obtenido la solución para { un+1 }, las velocidades y aceleraciones son actualizadas con las ecuaciones 3.60 y 3.61. Los parámetros de Newmark son relacionados así: 1 α = (1+ γ )2 , 1δ = + γ 3.64 4 2 donde γ = amplitud del facto de decaimiento La ecuación 3.64 es estable siempre y cuando γ ≥ 0 4 VALIDACION DEL MODELO 4.1 DESCRIPCION DEL MODELO COMPUTACIONAL La herramienta utilizada para desarrollar el modelo computacional por elementos finitos fue MATLAB. Matlab es el nombre abreviado de “MATrix LABoratory”, es un programa para realizar cálculos numéricos con vectores y matrices. Como caso particular puede también trabajar con números escalares, tanto reales como complejos. Una de las capacidades más atractivas es la de realizar una amplia variedad de gráficos en dos y tres dimensiones. Y también posee un lenguaje de programación propio. El modelo computacional fue desarrollado para realizar cálculos de desplazamientos y esfuerzos en placas cuando estas están sometidas cargas en el plano, a flexión o combinadas, de modo estático o dinámico. Fabricadas en materiales isotropicos, ortotropicos o laminados. Aspectos suficientes para el análisis en cuestión. Por medio la visualización grafica de matlab se puede analizar la distribución de desplazamiento y esfuerzos a lo largo de toda la placa. Con esto es posible encontrar con mayor facilidad cuales son las áreas en las cuales se dan los mayores o menores desplazamientos y esfuerzos. 4.2 VALIDACION ANALITICA A continuación se hará una comparación de los resultados obtenidos por una solución analítica vs solución por elementos finitos, para comprobar que la formulación del modelo arroja resultados coherentes. 4.2.1 Viga en voladizo cargada puntualmente en el extremo libre. Figura 18. Viga en voladizo cargada puntualmente en el extremo libre. P Datos del Problema L (m) 2 P (N) 2000 h h (m) 0.2 b (m) 0.20 A B b E (N/m2) 200e9 Calcular Desplazamiento en el punto B mm L Solución Analítica y P L 3 1 3 B = ; I = b ⋅h ; P τ = ; A = b ⋅ h; 3EI 12 xy A I 1= (0.2)(0.2)3 = 1.334×10−4 m4 12 y (2000).(2) 3 B = 9 = 0.0002 m → y = 0.2 mm(3).(200×10 ).(1.334×10−4 ) B Tabla 3. Resultados del análisis por elementos finitos N° Nº a:b Div Div Nodo Deflexión Error Elementos Nodos x y Analisis (mm) (%) 30 48 1:1.3 15 2 32 0.1644 17.8 34 54 1:1.17 17 2 36 0.1713 14.35 40 63 1:1 20 2 42 0.1787 10.70 50 78 1:1.25 25 2 52 0.1861 6.95 75 104 1.33:1 25 3 52 0.1862 6.90 90 124 1:1 30 3 62 0.1905 4.76 105 144 1:1.16 35 3 72 0.1932 3.42 160 205 1:1 40 4 82 0.1950 2.50 200 255 1:1.25 40 5 153 0.1951 2.45 225 276 1.1:1 45 5 138 0.1963 1.80 250 306 1:1 50 5 153 0.1972 1.40 Figura 19. Número de Elementos vs % de Error 20 18 16 14 12 10 8 6 4 2 0 30 34 40 50 75 90 105 160 200 225 250 Nº Elementos Figura 20. Número de Elementos vs Deflexión 0,207 0,198 0,189 0,18 0,171 0,162 0,153 0,144 30 34 40 50 75 90 105 160 200 225 250 Nº Elementos Deflexion (mm) % Error 4.2.2 Viga en voladizo con carga distribuida uniforme Figura 21. Viga en voladizo con carga distribuida uniforme w Datos del Problema L (m) 2 Solu w (N/m) 5000 h h (m) 0.20 ción A B b b (m) 0.20 Anal E (N/m2) 200 x 109 L Calcular ítica Desplazamiento en B mm For w L4mulas y ; I 1 b h3 ; PB = = ⋅ τ = ; A = b.h; 8EI 12 xy A I 1= (0.2)(0.2)3 = 1.334×10−4 m4 12 y (5000).(2) 4 = −4B 9 −4 = 3.75×10 m → yB = 0.375 mm(8).(200×10 ).(1.334 ×10 ) Tabla 4. Resultados del análisis por elementos finitos N° Nº a:b Div Div Nodos Deflexión Error Elementos Nodos x y Cargados (mm) (%) 30 48 1:1.3 15 2 34 - 48 0.3367 10.20 34 54 1:1.17 17 2 38 - 54 0.3474 7.36 40 63 1:1 20 2 44 - 63 0.3582 4.48 50 78 1:1.25 25 2 54 - 78 0.3683 1.78 75 104 1.33:1 25 3 80 - 104 0.3685 1.73 90 124 1:1 30 3 95 - 124 0.3738 0.32 105 144 1:1.16 35 3 110 - 144 0.3768 0.48 160 205 1:1 40 4 166 - 205 0.3786 0.96 200 255 1:1.25 40 5 207 - 246 0.3787 0.98 225 276 1.1:1 45 5 232 - 276 0.3797 1.25 250 306 1:1 50 5 257 - 306 0.3804 1.44 Figura 22. Número de Elementos vs % de Error 15 10 5 0 30 34 40 50 75 90 105 160 200 225 250 Nº Elementos Figura 23. Número de Elementos vs Deflexión 0,4 0,37 0,34 0,31 0,28 0,25 0,22 0,19 0,16 0,13 30 34 40 50 75 90 105 160 200 225 250 Nº Elementos Deflexion (mm) % Error 4.2.3 Viga en voladizo sometida a carga axial Figura 24. Viga en voladizo sometida a carga axial Datos del Problema L (m) 2 a (m) 0.20 P h P (N) 50000 b (m) 0.20 E (N/m2) 200 x 109 A B b Calcular Desplazamiento en el punto B mm L Esfuerzo en el punto B MPa P 50000 σ = → σ = = 1250000 → σ = 1.25MPa A 0.04 PL (50000)(2) δ = = ( = 1.25×10 −5 m → δ = 0.0125mm EA 200×109 )(0.04) Tabla 5. Resultados del análisis por elementos finitos N° Elementos Nº Nodos a:b Div x Div y δ Error Esfuerzos Error (mm) (%) MPa (%) 30 48 1:1.3 15 2 0.012250 2 1.25 0 34 54 1:1.17 17 2 0.012260 1.92 1.25 0 40 63 1:1 20 2 0.012240 2.08 1.25 0 50 78 1:1.25 25 2 0.012210 2.3 1.25 0 75 104 1.33:1 25 3 0.012340 1.28 1.25 0 90 124 1:1 30 3 0.012342 1.26 1.25 0 105 144 1:1.16 35 3 0.012340 1.28 1.25 0 160 205 1:1 40 4 0.012354 1.16 1.25 0 200 255 1:1.25 40 5 0.012388 0.89 1.25 0 225 276 1.1:1 45 5 0.012386 0.91 1.25 0 250 306 1:1 50 5 0.012384 0.92 1.25 0 Figura 25. Número de Elementos vs % de Error 5 4,5 4 3,5 3 2,5 2 1,5 1 0,5 0 30 34 40 50 75 90 105 160 200 225 250 Nº Elementos Figura 26. Número de Elementos vs Deflexión 0,0125 0,01247 0,01244 0,01241 0,01238 0,01235 0,01232 0,01229 0,01226 0,01223 0,0122 0,01217 30 34 40 50 75 90 105 160 200 225 250 Nº Elementos Deflexion (mm) % Error 4.2.4 Cilindro recto sometido a compresión Considere un cilindro sujeto a dos cargas a compresión en la mitad de su longitud, las cuales tienen un valor de 100 lb. Las propiedades y dimensiones del cilindro son E = 10.5 x 106 psi, v = 0.3125, y L = 10.35 pulg, r = 4.953 pulg y t = 0.094 pulg. El cilindro esta sujeto a lo largo de su longitud a 90º de la fuerza aplicada. Figura 27. Cilindro recto F l l 2 2 Y r Z X Z t F Tabla 6. Comparación de resultados Elemento Sol. Analítica Sol. Computacional Error (%) Shell 90 0.1139 0.1138 0.087 Shell 160 0.1139 0.1149 0.877 Solución por Matlab 512 Elementos 0.1139 0.1145 0.527 Figura 28. Desplazamientos en dirección Y con enmallado de 512 elementos Figura 29. Desplazamientos totales con enmallado de 512 elementos 4.2.5 Placa sujeta a una carga uniforme Considere una placa cuadrada (a=b) cuya longitud es de 0.5 m y tiene un espesor de 3 mm. La placa esta soportada simplemente en todos sus bordes y se le aplica una presión uniforme de 2500 Pa en su superficie. Las propiedades de la placa son: E1 = E2 = 70 GPa ν12 = 0.33 Tomando la ecuación B9 del anexo B tenemos: w p0a 4 max = λ D donde m=∞ n=∞ 16 (−1)(m+n) / 2−1 λ = ∑ ∑ 1 1 π 6m n(m2 + n2 )2 Tabla 7. Aproximación de la ecuación de la deflexión máxima m n λ 1 1 0.0041606 1 1,3 0.0041606 1,3 1,3 0.0041606 1,3 1,3,5 0.0041606 1,3,5 1,3,5 0.0041606 La solución de para la deflexión máxima para m=5 y n=5 es 4 wmax = 0.0041606 p0a D Introduciendo los datos del problema se tiene: D Et 3 (2500)(0.5)4 = 12(1−ν 2 ) = 176.7478 → wmax = 0.0041606 176.7478 Finalmente la deflexión máxima es: wmax = 3.678 mm Tabla 8. Resultados utilizados diferentes mallas Elementos Deflexión Error (%) 64 3.716 1.03 81 3.720 1.14 Figura 30. Desplazamientos en dirección Z Figura 31. Giro totales 4.2.6 Placa laminada sujeta a una carga uniforme Figura 32. Placa laminada sujeta a una carga uniforme p =500 Pa 1 mat 1 mm 0 cloth 2 cloth 1.4 mm mat 2 mm b =0.7 m 3 mat cloth 4 cloth 1.4 mm a =2.8 m 5 mat 1 mm (a) (b) Considere el caso de una placa rectangular con largo de a = 2.8 m y ancho de b = 0.78 m que esta sujeta a una presión uniforme de 500 Pa. La placa esta hecha de un laminado simétrico de cinco capas. La capas 1,3 y 5 son capas de mat reforzadas de 450 g/m2, y sus propiedades son: E1 = E2 = 7.72 GPa ν12 = 0.33 G12 = 2.91GPa La capa 3 es doble. Las capas 2 y 4 son dobles con reforzados con cloth de 500 g/m2 y sus propiedades son: E1 = E2 = 13.8 GPa ν12 = 0.12 G12 = 1.87 GPa Calculando la matriz elástica (E) para cada material, se obtiene: Para las capas de mat 1, 3, 5: E11 = 8.66 GPa E12 = 2.86 GPa E22 = E11 = 8.66 GPa E66 = 2.91GPa Para las capas cloth 2, 4: E11 =14 GPa E12 = 1.68 GPa E22 = E11 = 14 GPa E66 = 1.87 GPa La matriz elástica para un laminado se obtiene con la tercera ecuación de 2.62: C11 = 272.64 Nm C12 = 64.834 Nm C22 = C11 = 272.64 Nm C66 = 67.358 Nm De la ecuación C7, el coeficiente D2m−1,2n−1 esta expresado como: D 4 2 22m−1,2n−1 = 272.64 (2m −1) + 6385.64 (2m −1) (2n −1) + 69796.68(2n −1) 4 Tabla 9. Valores de los coeficientes α como una función de n y m m n α x 10-5 (Pa-1) 1 1 1.308 2 2 1.0807 3 3 1.1302 4 4 1.1169 5 5 1.1213 6 6 1.1196 7 8 1.1203 8 9 1.1199 10 20 1.1200 15 15 1.1201 20 20 1.1201 La máxima deflexión wmax es deducida de la ecuación C18 Los valores de α en la serie como una función de m y n están dados en la Tabla 9, en la cual α = 1.201×10−5 Pa−1 Donde la máxima deflexión con a=2.8 m y p0=500 Pa: wmax = 5.728 mm Tabla 10. Valores de los coeficientes β y γ como función de m y n m n β x 10-6 (Pa-1) γ x 10-5 (Pa-1) 1 1 13.08 1.308 2 2 -6.904 1.0484 3 3 5.383 1.1000 4 4 -1.1328 1.0865 5 5 2.4327 1.0906 6 6 3.5091 1.0890 7 7 1.6454 1.0897 8 8 0.7943 1.0893 10 10 0.9604 1.0894 15 15 1.1910 1.0895 30 30 1.1313 1.0894 60 60 1.1371 1.0894 100 100 1.1378 1.0894 150 150 1.1379 1.0894 El calculo del esfuerzo (C24) necesita la determinación de las series β (C21) y γ (C22), los cuales son determinados en la tabla 10 β =1.1379 ×10−6 Pa−1 γ = 1.894 ×10−5 Pa−1 Los esfuerzos en el plano en la capa mat , en el centro de la placa, son: σ m mxx = Axx z σ m myy = Ayy z σ mxy = 0 donde, Am 16 p 2 0a m xx = 4 (E11 β + Em R2γ )π 12 Am 16 p a 2 yy = 0 4 (Em12 β + Em R222 γ )π Introduciendo todos los valores, se obtiene: Amxx = 3.2721×10 8 Nm−1 Amyy = 9.774 ×10 8 Nm−1 Similarmente los esfuerzos en el cloth son: σ cxx = A c xx z σ c cyy = Ayy z σ cxy = 0 y, 2 Ac 16 p axx = 0 (E c4 11 β + Ec 2π 12 R γ ) Ac 16 p a 2 yy = 0 (E c4 12 β + Ec R2γ ) π 22 Introduciendo los valores en A se tiene: Amxx = 1.9883×10 8 Nm−1 Amyy = 1.5727 ×10 8 Nm−1 Tabla 11. Variación de esfuerzos en el centro de la placa. 1.112 3.313 3 mat 3 mat 0.477 3.774 0.785 2.3339 2 2 z (mm) cloth z (mm) cloth 0.199 1 0.327 0.974 1 1.573 mat mat σ x x (MPa) σ y y (MPa) Los esfuerzos transversales cortantes son determinados por: ∂ ∂ ∂ σ xx + σ xy + σ = 0∂x ∂y ∂z xz ∂ ∂ ∂ σ xy + σ + σ = 0∂x ∂y yy ∂z yz El esfuerzo cortante en un punto (x,y) de la placa son escrito de la forma: σ kx z = B k x z ( x, y) z2 + const σ k ky z = By z ( x, y) z 2 + const k = m,c Las constantes son determinadas considerando la continuidad de los esfuerzos y ellas desaparecen en las caras inferior y superior. Capa 1 (mat): σ 1i z = B m i z ( x, y)( z2 − h20 ) , i = x, y h0 = 3.4mm Capa 2 (cloth) σ 2 = Bci z i z ( x, y)( z2 − h21 ) + Bm 2 2i z ( x, y)(h1 − h0 ) i = x, y h0 = 2.4mm Capa 3 (mat) σ 2 = Bm ( x, y)( z2 − h2 2 2i z i z 2 + h1 − h0 ) + Bci z ( x, y)(h22 − h21 ) i = x, y h0 = 1mm La distribución de esfuerzos en las capas 4 y 5 son simétricas a la distribución en las capas 1 y 2. Tabla 12. Resultados del análisis por elementos finitos N° Elementos Nº Nodos a:b Divisiones Divisiones δ Error x y (mm) (%) 324 370 1:1 9 36 5.6 1.75 Figura 33. Desplazamiento en Z Figura 34. Giros Totales Figura 35. Esfuerzos en dirección X en la capa 1 Figura 36. Esfuerzos en dirección X en la capa 5 Figura 37. Esfuerzos en dirección Y en la capa 1 Figura 38. Esfuerzos en dirección Y en la capa 5 Figura 39. Esfuerzo cortante XY en la capa 1 Figura 40. Esfuerzo cortante XY en la capa 5 4.2.7 Análisis de resultados. Antes de analizar los resultados es necesario aclarar que para los ejercicio de vigas sometidos a carga a flexión es mejor trabajar con elementos tipo barra, los cuales dan mejores resultado que los elementos cuadriláteros que se utilizaron en el modelo. Sin embargo, para efectos de validación son útiles estos elementos (cuadriláteros) ya que permiten comparar la veracidad de los resultados de estos. En los ejercicios de vigas desarrollados anteriormente que estaban sometidos a cargas de flexiones se noto que los resultados de desplazamientos más cercanos al valor real, se dieron en mallas mas refinadas. Sin embargo un refinamiento demasiado grande de la malla ocasiona un aumento en el error y además un gasto innecesario de tiempo de cómputo. En la sección 4.2.3 cualquiera de las malla podría considerarse como optima siempre y cuando no sea de interés conocer de manera gráfica la distribución de esfuerzo, debido a que con pocos elementos no se aprecia una correcta distribución de esfuerzos. A medida que se fue aumentado el número de elementos en la malla, la distribución de esfuerzos se hizo más clara. En los ejercicios de placa simple y laminada en las graficas se ven deformadas por efecto de giro y desplazamiento, las cuales han sido exageradas para efecto de visualización, sin embargo en la realidad estos valores no son observables a simple vista. En las graficas de giros da la impresión que los nodos de los bordes presentaran giros, sin embargo los que giran son los nodos proximos a a estos. 4.3 VALIDACION DINAMICA CON UN MODELO REAL Los datos reales se obtuvieron por medio de una prueba realizada a una placa en voladizo la cual se deflectó hacia abajo y luego se suelta. Figura 41. Estado inicial de la placa o o Posición de Equilibrio Deflexión Inicial Una vez se suelta el extremo libre este empieza a vibrar hasta quedar en reposo. Los desplazamientos que se dan durante la vibración son registrados mediante un dispositivo electrónico. El dispositivo electrónico utilizado para registrar los desplazamientos a través del tiempo en el extremo libre de la placa consta básicamente de dos etapas. Figura 42. Comportamiento dinámico de la placa Máx. Elongación o o Posición de Equilibrio Min Elongación La primera etapa se encarga de crear una señal que define los desplazamientos, esta señal resulta de la variación en la cantidad de luz recibida por una fotorresistencia que se encuentra en un costado del extremo libre de la placa, y la luz que esta recibe proviene de un haz de luz emitido por una pequeña lámpara localizada al frente de la fotorresistencia, es decir al otro costado del extremo libre de la placa en la posición de equilibrio de esta. Al vibrar la placa permite que la cantidad de luz directa que recibe la fotorresistencia varié y la resistencia de esta también, este cambio de resistencia luego es enviado a un “Puente de Winstone” el cual se encarga de convertir la variación de resistencia en una variación de voltaje y estos voltajes son utilizados para medir los desplazamientos Figura 43. Adquisición de datos o o Fotorresistencia Emisor de Luz Salida de señal Puente de Winstone . La segunda etapa se encarga de recibir la variación de voltaje y visualizarlos, dicha visualización se lleva a cabo mediante un osciloscopio que recibe la señal emitida en la primera etapa. Previamente se toman valores de voltajes correspondientes a determinados desplazamientos y con estos hallamos una ecuación por regresión lineal para poder determinar el desplazamiento correspondiente a cada valor de voltaje mostrado por la grafica que saca el osciloscopio. 4.3.1 Validación dinámica en una lámina de PRFV. La lámina en voladizo esta fabricada por un laminado de tres capas fibra (mat+roving+mat). Cada una de las láminas es considerada isotrópica y por tanto el ángulo de orientación respecto a los ejes de referencia es cero para cada lámina. Tabla 13. Datos de la lámina Dimensiones Material L (m) t (m) b ( m) ρ (kg/m3) PRFV 0.25 0.0032 0.04 1593.75 Propiedades Capas E (GPa) v t (mm) Mat 450 7.72 0.33 1.1 Roving 13.34 0.12 1.0 Tabla 14. Resultados de la prueba Tiempo Desplazamiento (ms) (cm) 0 3.0712 40 -2.637 80 2.5702 130 -2.136 160 1.9022 210 -1.7686 250 1.4346 300 -1.3344 330 1.0672 380 -0.8334 420 0.633 460 -0.466 500 0.2656 El rango de tiempo considerando en la tabla 13 (500 milisegundos) es un tiempo suficiente para analizar la rigidez del material en cuestión, además permite hallar los coeficientes de amortiguamiento (coeficientes de Rayleigh), los cuales son indispensables para el análisis dinámico de cualquier material. Después de los 500 milisegundos se siguieron dando oscilaciones mucho más pequeñas, que no son tenidas en cuenta debido a que no son indispensables para el calculo de los coeficientes de amortiguamiento, además el tenerlas en cuentas implica más tiempo de análisis y a su vez un aumento en el tiempo de computo en el momento de modelar el efecto. Figura 44. Variación del desplazamiento en el tiempo. 3.1 2.6 2.1 1.6 1.1 0.6 0.1 -0.4 0 40 80 130 160 210 250 300 330 380 420 460 500 -0.9 -1.4 -1.9 -2.4 -2.9 -3.4 4.3.2 Calibración de los coeficientes de Rayleigh. En la sección anterior se explicó como obtener los desplazamientos respecto al tiempo de una placa empotrada en un extremo. Haciendo uso de esto es posible obtener los coeficientes α y β, y así calcular el amortiguamiento de Rayleigh descrito en la sección 3.8.1 Introduciendo todos los datos anteriores en la formulación por elementos finitos, queda por último calcular cuales son los coeficientes α y β que dan la mejor aproximación a los desplazamientos reales. Tras de pruebas de ensayo y error los mejores resultados se obtuvieron con los valores de 0.01 y 0.0008 para α y β. Figura 45. Comparación de resultados 3.1 3.1 2.6 Simulación en Matlab Prueba Real 2.6 2.1 2.1 1.6 1.6 1.1 1.1 0.6 0.6 0.1 0.1 -0.4 0 40 80 130 160 210 250 300 330 380 420 460 500 -0.4 0 42 83 125 167 209 250 292 334 376 417 458 500 -0.9 -0.9 -1.4 -1.4 -1.9 -1.9 -2.4 -2.4 -2.9 -2.9 -3.4 -3.4 Tabla 15. Resultados de pruebas Prueba Real Prueba Matlab Tiempo Desplazamiento Tiempo Desplazamiento (ms) (cm) (ms) (cm) 0 3.0712 0 3.0 40 -2.637 42 -2.65 80 2.5702 83 2.45 130 -2.136 125 -2.02 160 1.9022 167 1.98 210 -1.7686 209 -1.83 250 1.4346 250 1.65 300 -1.3344 292 -1.4 330 1.0672 334 1.1 380 -0.8334 376 -0.98 420 0.633 417 0.84 460 -0.466 458 -0.62 500 0.2656 500 0.38 5 APLICACIÓN DEL MODELO 5.1 DESCRIPCION DE LA PRUEBA REAL Una de las aplicaciones del modelo desarrollado en los capítulos anteriores consistirá en analizar el comportamiento de tanques fabricados en PRFV, para lo cual se utilizaran las pruebas de una empresa que trabaje con PRFV. En este caso se trabajará con Rotofibra Ltda. quienes se dedican a la fabricación de recipientes y tuberías en PRFV. Dentro de sus productos de mayor demanda se encuentra los tanques para almacenamiento de líquidos. En esta entidad cada nuevo diseño de tanque es llenado de agua y soltado a una altura promedio de 1,2 metros (altura promedio de un camión de carga) para evaluar su resistencia mecánica. Esta prueba es aplicada a tanques hechos en plásticos, sin embargo es de interés para la empresa conocer el comportamiento de tanques hechos en PRFV cuando se ven sometidos a una situación similar. El estudio consistirá en simular el comportamiento de tanques fabricados en PRFV cuando son utilizados para almacenar un líquido (en este caso agua) a presión atmosférica, lo cual permitirá observar la distribución de desplazamiento y esfuerzos en el tanque. También se simulará el comportamiento de estos mismos tanques (llenos de agua), pero sometidos a una carga de impacto provocada por la cada de este a una determinada altura. 5.2 MATERIALES DEL TANQUE Rotofibra Ltda. utiliza como materiales principales para la fabricación de tanques en PRFV los siguientes tipos de fibra VFG-Woven Roving 800 y Mat 450, y como resina la Cristalán 805. 5.2.1 VFG-Woven Roving 800. Este producto es un tejido hecho de filamentos continuos colocados en forma vertical y horizontal, sin amarras, para laminación manual. Por lo cual puede tomar distintas formas y curvas. Siendo un tejido pesado podrá transferir su configuración a través del Gel Coat si es colocado cerca de la superficie. Entre sus aplicaciones se encuentran la producción de lanchas, cavas, piezas automotrices, gabinetes y artículos del hogar, implementos deportivos, láminas translúcidas, tanques resistentes a la corrosión y un sin número de artículos de diversos tamaños y áreas. 5.2.2 Mat 450. Es una colchoneta de hilos cortados, distribuidos multidireccionalmente en un solo plano. Es ideal para moldeo abierto a baja presión, como soporte de Geal Coat y como un refuerzo para el laminado estándar. Se usa en sistemas de resinas termofijas, en moldeo manual y como refuerzo para laminado. Algunos usos concretos son: producción de lanchas, cavas, tanques resistentes a la corrosión, piscinas, componentes para camiones, paneles para la construcción, ductos, tuberías y muchas otras partes. 5.2.3 Cristalán 805. Es una resina poliéster ortoftálica muy versátil, para aplicaciones de propósito general, con propiedades mecánicas excelentes y buen desempeño en aplicaciones eléctricas. Es apropiada para laminación manual o moldeo a máquina. Los laminados tienen alta resistencia mecánica, buena rigidez y excelente resistencia al impacto. Se desempeñan bien en más de 200 ambientes químicos. El Cristalán 805 es apropiado para la fabricación de embarcaciones, cabinas, partes para carros, camarotes, cubiertas para máquinas, recipientes para líquidos, etc. 5.6 SIMULACION ESTATICA La simulación se hará a un tanque para almacenar agua a presión atmosférica, el cual tiene una capacidad de almacenamiento de 10000 litros. Tabla 16. Datos del tanque Dimensiones Material L (m) t (m) D ( m) ρ (kg/m3) PRFV 2.45 0.0032 2.3 1593.75 Propiedades Capas E (GPa) v t (mm) Mat 450 7.72 0.33 1.1 Roving 13.34 0.12 1.0 Figura 46. Distribución de fuerzas en tanque de 10000 Litros 0 m t 4234.16 N 0.245 m 8468.32 N 0.49 m Mat 12702.48 N 0.74 m P Rovin R g O Mat 16936.63 N 0.98 m F U 21170.79 N N 1.22 m D I Φ = 2.30 m 25404.95 N 1.47 m D A D 29639.11 N 1.72 m 33873.27 N 1.96 m 38072.86 N 2.20 m 141681.64 N 2.45 m 5.3.1. Resultados del análisis. Figura 47. Desplazamientos en dirección X y Y Dirección x Dirección y Figura 48. Desplazamientos totales Figura 49. Esfuerzos normales Figura 50. Esfuerzos cortantes 5.7 SIMULACION DINAMICA Para la simulación dinámica se utilizará como modelo el mismo tanque de la sección 5.3, el cual se dejará caer de una altura de 1,2 metros. Se analizará el momento de impacto donde la velocidad en ese instante es de 4.8522 m/s. El tiempo a analizar después del impacto será alrededor de 0.2 segundo, tiempo suficiente para observar el comportamiento dinámico de la estructura. Se estarán utilizando los mismo coeficientes de Rayleigh calculados en la sección 4.32 y con diferencial de tiempo de 4 x10-4, espacio de tiempo que permite una resolución que refleja todo el comportamiento en el tiempo dado (0.2 seg.) Para este tanque se analizaran tres zonas; inferior, media y superior. Figura 51. Desplazamiento vs tiempo en la zona inferior Figura 52. Desplazamiento vs tiempo en la zona media Figura 53. Desplazamiento vs tiempo en la zona superior 5.4.1 Análisis de resultados. En primera medida los desplazamiento en direcciones “x” y “y” son iguales gracias a la simetría de las cargas desarrolladas a lo largo del tanque, también se puede observar que en dirección “z” se desarrollaron desplazamientos, aunque esto fueron pequeños comparado con los de “x” y “y”. Los mayores desplazamientos se dieron en la parte inferior del tanque, ya que en esta se localizan las presiones más altas. Por tanto, esta zona suele ser crítica, no olvidando que los desplazamientos en esta zona contribuyen que hallan desplazamiento en dirección z a lo largo de todo el tanque. Los mayores desplazamientos dados en el tanque están alrededor de los 1.5 cm. lo cual es un desplazamiento alto y crítico, por esto la empresa Rotofibra Ltda. utiliza rigidizadores (anillos de refuerzos) a lo largo del tanque. En este caso para un tanque de 10.000 litros utilizan dos rigidizadores. Se puede observar que los esfuerzos normales tienen un comportamiento similar en cuanto a simetría a los desplazamientos en dirección “x” y “y” De los tres puntos analizados en la sección 5.4, se nota que el punto ubicado en la parte superior presenta mayores oscilaciones, ya que está más lejos del soporte del tanque. Los puntos ubicados en la zona inferior y media del tanque sigue vibrando después de los 0.2 segundos, pero sus oscilaciones son considerablemente pequeñas. Observando las pocas oscilaciones en cada punto analizado podemos darnos cuenta que el material de fabricación del tanque es bastante rígido. 5.8 ESTUDIO COMPARATIVO El propósito de esta sección es dar una idea del comportamiento de algunos tanques cuando se varía la configuración del laminado y geometría que fabrican en la actualidad Rotofibra. El estado de fuerzas al que estarán sometidos los diferentes tanques será similar al mostrado en la figura 46 pero de acuerdo a la capacidad de cada uno de ellos. Tabla 17. Comparación de resultados Capacidad de Dimensiones Tanque con Tanque sin Almacenamiento Configuración del tapa tapa (Litros) D L laminado (Desp max) (Desp max) (m) (m) (mm) (mm) Mat 3.00 3.00 1000 0.9 1.57 Mat+Mat 2.10 2.10 Mat+Roving+Mat 1.70 1.70 Mat+Mat 15.00 14.90 5000 1.7 2.2 Mat+Roving+Mat 8.30 8.20 Mat +Roving+Mat+Roving+Mat 4.80 4.70 Mat+Mat 52.60 52.60 20.000 2.8 3.0 Mat+Roving+Mat 29.00 29.00 Mat +Roving+Mat+Roving+Mat 16.50 16.70 5.5.1 Análisis de resultados. En base a los resultados mostrados en la tabla 17 se puede notar claramente la disminución de desplazamiento a medida que se aumenta el número de capas del laminado. También se puede obtener variación en el desplazamiento si la dirección de las mantas de fibra en el laminado tuvieran diferentes ángulos de orientación, pero para este caso las mantas de fibra Mat 450 y Roving 850 tienen propiedades como módulos de elasticidad y Poisson iguales a nivel transversal y longitudinal lo que hace que se que el laminado se comporte igual pesar de variar el ángulo de una fibra con respecto a la otra, por lo cual asumimos ángulo de 0º en el momento de hacer los análisis, llegando a la conclusión que la manta de fibra Mat tiene un comportamiento isotrópico al igual que la manta de fibra Roving siempre que esta se someta a fuerzas en direcciones ortogonales. Comparando los datos de las dos últimas columnas de la tabla 17 se observan los resultados de un tanque de un mismo laminado que se someten al mismo estado de fuerza, pero uno con tapa en la parte superior y el otro no, y a pesar de esta diferencia en su estructura la respuesta de los tanque en la mayoría de los casos es igual y donde no, la diferencia es mínima. Esto nos deja ver que la tapa en un tanque de almacenamiento de líquido, de agua en este caso, no afecta de forma considerable el comportamiento del mismo cuando se somete a una carga estática. Sin embargo por lo general se le fabrican con tapa para evitar que el líquido almacenado se vea contaminado. Y al saber que la tapa no tiene mayor incidencia en el tanque, estas se diseña a partir del criterio del diseñador el cual tendrá en cuenta otros factores, como es el caso de la estética y la corrosión, de hecho el que la tapa sea mas alta en el centro, impide la acumulación de agua o cualquier sustancia y así se evita que el material se descomponga con el transcurrir del tiempo. La zonas mas afectadas serán las mas próximas al borde inferior donde reposa el tanque, de hecho es en este sitio donde se dan los máximos desplazamientos, los cuales se desarrollan en las dirección “x” y “y”. Este comportamiento se da de igual forma como se mostró en las figuras 47 a 50 CONCLUSIONES ¸ A través de todo el estudio se puedo evidenciar que el método de elementos finitos sigue siendo una herramienta para la optimización de los diseños, tal como lo es para el caso de análisis de este trabajo, donde uno de los objetivos era conocer el comportamiento estructural de tanques fabricados en PRFV cuando se sometieran a cargas estáticas y dinámicas. ¸ Otros programas especializados en esta área utilizan elementos Shell para analizar este mismo tipo de problemas, sin embargo con este estudio se pudo demostrar que un elemento de placa cuadrilátero de cuatro nodos funciona bien en geometrías curvas, además que su formulación resulta mas sencilla que la del elemento shell. ¸ Al analizar los resultados del comportamiento con diferentes laminados se noto con claridad que ha medida que se aumentaba el espesor los desplazamiento eran menores, lo cual a primera vista deja entrever que lo mejor es aumentar el espesor del laminado para evitar grandes deformaciones, sin embargo en realidad esto resulta antieconómico ya que esto requiere mas material. Por lo cual, se hace necesario una solución para no tener que aumentar la cantidad de material sin necesidad tampoco de recurrir a otro material que resultaría mas costoso a mediano y largo plazo por su mantenimiento. ¸ El análisis estático en los tanques dejo claro que la zona critica de estos se localiza en los bordes de la parte inferior, por lo cual requieren atención a la hora de fabricarlo y también al momento de la instalarlo, asegurándose que su base quede bien sujeta y además la superficie donde se instala debe estar debidamente nivelada y libre de cualquier objetos para que no se creen cargas en el fondo del tanque, las cuales son bastante perjudiciales para este. ¸ Para controlar desplazamientos es oportuno hacer uso de rigidizadores a lo largo de los tanques. La ubicación de estos rigidizadores dependerá de la geometría y las condiciones a que este sometido el tanque. Hay que tener en cuenta que el uso de estos rigidizadores se debe hacer con una buena elección del espesor del laminado, teniendo cuidado de no saturar el tanque con estos refuerzos de tal forma que se este sobrediseñando o encareciendo le producto. ¸ El análisis del dinámico es de vital importancia, ya que este nos permite observar como responde un tanque cuando se somete a una carga de impacto. Con esta información se pueden determinar los puntos que mas sufren daño un instante después de recibir la carga de impacto, lo cual permite revaluar el diseño de toda la estructura o de las partes mas afectadas para así evitar problemas de rotura. ¸ Los resultados de análisis dinámico dejan ver que la fibra de vidrio es un material compuesto bastante rígido comparado con otros materiales de uso común como el acero. El estudio de este comportamiento es de mucha utilidad para el diseñador ya que le permite cerciorarse que el laminado y la disposición de este sea la adecuada para conseguir determinado comportamiento en su diseño. RECOMENDACIONES ¸ Para la programación de elementos finitos es recomendable eliminar ciclos innecesarios así como variables excesivas o repetidas, ya que todo esto provoca en la solución de los problemas lleve mas tiempo de computo. ¸ Es aconsejable la utilización de un software que facilite la programación como lo es Matlab, el cual trae consigo gran cantidad de subrutinas que lo hacen fácil de trabajar. ¸ Una manera segura de depurar el algoritmo es: primero revisar que los conceptos estén bien aplicados, segundo comenzar a revisar del algoritmo más simple al más complejo. BIBLIOGRAFIAS BERTHELOT, Jean Marie. Composite Materials Mechanical Behavior and structural Analysis. Springer, 1998. 645 p. CHANDRUPATLA, Tirupathi. Introducción al Estudio de Elemento Finito en Ingeniería. 2 ed. Mexico: Pearson, 1999, 462 p. DEREK, Hull. Materiales Compuestos. Barcelona: Reverté SA, 1987. 254 p. DONALDSON, Bruce K. Análisis of Aircraft Structure an Introduction, New York: Mc-Graw Hill, 1993, 935 p ORTIZ BERROCAL, Luis. Elasticidad. 3 ed. Madrid: Mac Graw Hill, 1998. 549 p. REDDY, J.N. An Introduction to the Finite Element Method. 2 ed. Mc. Graw Hill, 1993, 672 p. UGURAL, Ansel C. Stresses in Plates and Shells. 2 ed. New Jersey: Mac Graw Hill, 1981. 502 p. VINSON, J.R. The Behavior of Structure Composed of Composite Materials. 1 ed. Boston: Martinus Nijhoff Publishers, 1987. 323 p. VOLTERRA, Enrico. Advanced Strength of Materials, New Jersey: Prentile-Hall, Inc., 1971, 522 p. Zienkiewics, O. El Método de los Elementos Finitos, Mecánica de sólidos y fluidos Volumen 1. 4 ed. Barcelona: Mc-Graw Hill, 1989, 678 p. Zienkiewics, O. El Método de los Elementos Finitos, Mecánica de sólidos y fluidos Volumen 2. 4 ed. Barcelona: Mc-Graw Hill, 1989 , 510 p. ANEXO A. TRANSFORMACION DEL DIFERENCIAL DE AREA (dA) Se necesita expresar dxdy en las coordenadas locales. Por tanto, considere un mapeo de variables de “x” ,y “y” a u1 y u2 , dado como x = x(u1,u2 ) y = y(u1,u2 ) A1 Suponemos que las ecuaciones anteriores pueden intervine para expresar u1, u2 en términos de “x” y “y “que la correspondencia es única. Si una partícula se mueve desde un punto P, de tal manera que u2 se mantiene constante y u1 varía, entonces se genera una curva en el plano. Llamamos a ésta la curva u2 es generada manteniendo u1 constante y permitiendo que u2 varíe. Sea r el vector de posición de un punto P, Figura A1. Vector posición en un punto P u2 curva T2 dA y S2 T1 ds2 P S u1 curva 1 ds1 j x i r = xi + yj A2 donde i, j son vectores unitarios a lo largo de x, y respectivamente. Considere los vectores T ∂r ∂r1 = T∂u 2 = A3 1 ∂u2 o, en vista de que existe la ecuación A2, T ∂x i + ∂y j T ∂x i + ∂y1 = = j ∂u1 ∂u 2 1 ∂u2 ∂u A4 2 Podemos demostrar que T1 es un vector tangente a la curva u1 y que T2 es tangente a la curva u2 (figura A1). Para analizar estos usamos la definición ∂r lím ∆r= A5 ∂u ∆u→01 ∆u1 donde ∆r = r (u1 + ∆u1 ) − r (u1 ) . En el límite, la cuerda ∆r llega a ser la tangente a la curva u1 (figura A2). Sin embargo, ∂r ∂u1 (∂r ∂u2 ) no es un vector unitario. Para determinar su magnitud, escribimos Figura A2. Vector tangente. P’ P ∆r r r+∆r y x ∂r ∂r ds = 1 A6 ∂u1 ∂s1 du1 donde s1 es la longitud del arco a lo largo de la curva u1 y ds1 es el diferencial de la longitud de arco. La magnitud del vector ∂r = lím ∆r ∂s ∆u→01 ∆s1 es la razón límite de la cuerda a la longitud del arco que es igual a la unidad. Concluyendo que la magnitud del vector ∂r ∂u1 es ds1 du1 . Tenemos  T ds1   t T ds  1 =  = 2 t A7  du  1 2  1   du  2 2  donde t1 y t2 son vectores unitarios tangentes a las curvas u1 y u2 , respectivamente. Usando la ecuación A7. tenemos la representación de los vectores ds1 y ds2 cuyas longitudes son ds1 y ds2 ds1 = t1ds1 = T1du1 A8 ds2 = t2ds2 = T2du2 El diferencial de área dA es un vector con magnitud dA y dirección normal al área elemental, que en este caso es k. El vector dA, una vez observadas las ecuaciones A4 y A8 está dado por la regla del determinante dA = ds1 x ds2 = T1 x T2 du1 du2 i j k ∂x ∂y = 0 du du ∂u1 ∂u 1 2 1 ∂x ∂y 0 ∂u2 ∂u2  ∂x ∂y ∂x ∂y  =  − du du k  ∂u1 ∂u2 ∂u2 ∂u  1 2 1  Marcamos la matriz jacobiana mediante  ∂x ∂y   ∂u ∂u  J =  1 1  A9  ∂x ∂y   ∂u  2 ∂u2  La magnitud puede escribirse ahora como dA = det J du1 du2 A10 que es el resultado deseado. Note que si se trabaja con coordenadas ξ ,η en vez de u1 y u2 , entonces dA = det J dξ dη A11 ANEXO B. SOLUCION DE NAVIER PARA UNA PLACA RECTANGULAR SOPORTADA SIMPLEMENTE Y CARGADA UNIFORMEMENTE Considérese una placa rectangular de lados a y b, simplemente soportada en todo sus bordes y cargada uniformemente p(x,y). En general, la solución del problema de flexión es empleando las series de Fourier para carga y deflexión m=∞ n=∞ p ( x, y) ∑ ∑ p mπ x nπ y= mnsen sen B1 1 1 a b m=∞ n=∞ w( x, y) = ∑ ∑ amnsen mπ x sen nπ y B2 1 1 a b donde pm n y amn representan los coeficientes a determinar. La deflexión debe satisfacer la ecuación diferencial (2.20) con las siguientes condiciones de contorno 2 w = 0 ∂ w2 = 0 ( x = 0, x = a)∂x 2 w ∂ w= 0 2 = 0 ( y = 0, y = b)∂y La solución correspondiente a la carga p(x,y) requiere la determinación de los coeficientes de pm n y amn. Para determinar el coeficiente de Fourier, cada lado de la ecuación B1 es multiplicado por sen m 'π x sen n 'π y dx dy B3 a b e integrando entre los límites de 0,a y 0,b: b a ∫ ∫ p ( x, y)sen m 'π x sen n 'π y dx dy = 0 0 a b m=∞ n=∞ b a ∑ ∑ p p ( x, y)sen mπ x sen nπ y mn ∫ dx dy0 ∫0 1 1 a b Por integración directa esta es: a sen mπ x m 'π x 0 (m ≠ m ') ∫ sen dx =0 a a a / 2 (m = m ') B4 a sen nπ ysen n 'π y  dy =  0 (n ≠ n ') ∫0 b b b / 2 (n = n ') Los coeficientes de la doble expansión de Fourier son: b a p ( x, y)sen mπ xsen m 'π x∫ ∫ dx dy B5 0 0 a a La evaluación de amn en B2 requiere la sustitución de las ecuaciones B1 y B2 en 2.20, dando como resultado: m=∝ n=∝   m 4∑ ∑ a  π  2 m 2 2 4 π   nπ   nπ   p    + +  − mn  mπ x mn         sen sen nπ y = 0 1 1   a   a   b   b   D  a b Esta ecuación debe aplicarse para todas las x y y. Por tanto,  m2 n2 2  a 4 pmnπ  2 + 2  − mn = 0  a b  D donde, a 1 pmn = mn4 2 π D  2 2  B6 (m / a) + (n / b)  Sustituyendo la ecuación B6, se encuentra que la deflexión en la placa es: 1 m=∞ n=∞w ∑ ∑ pmn sen mπ x sen nπ y= 4 2 π D 1 1 (m / a)2 + (n / b)2  a b B7   en la cual pmn esta dado pro la ecuación B5. La ecuación B7 es una solución valida para placas rectangulares son soporte simples sometidas a flexión bajo varias clases de cargas. Cuando una placa rectangular esta sujeta a una carga uniformemente distribuida p(x,y)=p0, el resultado anterior se simplifica considerablemente. Ahora la ecuación B5 después de la integración es p 4 p0mn = 2 (1− cos mπ )(1− cos nπ ) 4 p = 0 1− (1)m  1− (−1)n  π m n π 2mn     o p 4 p0mn = 2 (m,n = 1,3...) m n π Se observa que pmn=0 para incluso valores de m y n, las integrales asumen solo los valores impares. Introduciendo pmn en B7 16 p m=∞ n=∞0 sen (mπ x / a)sen (nπ y / bw )= 6 (m, n = 1,3...) π D ∑ ∑ 2 2( ) ( B8 1 1 mn  m / a + n / b) 2   Claramente, para cargas uniformemente en placas debe deflectarse en forma simétrica. La máxima deflexión ocurre en el centro de la placa (x=a/2, y=b/2) entonces: m=∞ n=∞ (−1)(m+n) / 2−1w 16 p0max = π 6 D ∑ ∑ 2 2 2 B9 1 1 mn  (m / a) + (n / b)  en la ecuación B8 sen (mπ / 2) y sen (nπ / 2) son reemplazados por (−1)(m−1) / 2 y (−1)(n−1) / 2 respectivamente. ANEXO C. SOLUCION DE UNA PLACA LAMINADA SOPORTADAS SIMPLEMENTE Y SUJETA A UNA CARGA UNIFORME Considérese una placa rectangular sujeta a una carga distribuida p=p(x,y), en el caso de una placa ortotrópica la ecuación es ∂4w ∂4D w ∂ 4w 11 ∂x4 + 2( D12 + D66 ) ∂x2∂2 + Dy 22 = p ( x, y) C1 ∂y4 En el caso general de carga transversal puede expandido como una doble serie de senos de Fourier: m=∞ n=∞ p ( x, y) ∑ ∑ p sen mπ x= mn sen nπ y C2 1 1 a b donde el coeficiente qmn esta dado por: 4 b apmn = p ( x, y)sen mπ x sen nπ y dx dy a b ∫0 ∫0 a b C3 La solución del problema de flexión puede ser investigado expresando los desplazamientos en forma una doble serie de Fourier : m=∞ n=∞ w( x, y) = ∑ ∑ a sen mπ x sen nπ ymn C4 1 1 a b La expresión del coeficiente amn es obtenida sustituyendo la ecuación C4 en la ecuación C1 y escribiendo p(x,y) igual a la ecuación C2, es obtiene p π 4 am n = mn m 4 2 2 4 D   + 2( D + D ) m   n  + D  n  C5 11  a  12 66  a    22       b   b  La deflexión en un punto (x,y) es: a4 m=∞ n=∞ aw( x, y) = ∑ ∑ mn4 sen mπ x sen nπ y C6 π 1 1 Dmn a b donde D 4mn = D11m + 2( D12 + D66 ) (mnR) 2 + D 422 (nR) C7 y R = a/b. Los esfuerzos en el plano son deducidos: a 2  m=∞ n=∞k pσ = z ∑ ∑ mn (E k 2 k 2 2x x   11m + E12n R )sen mπ x sen nπ y C8  π  1 1 Dmn a b 2 σ k =  a  m=∞ n=∞ pz ∑ ∑ mn (E k 2 k 2 2y y   12m + E22n R )sen mπ x sen nπ y C9  π  1 1 Dmn a b  a 2  m=∞ n=∞k k p2 z R E ∑ ∑ mn mn cos mπ x cos nπ yσ x y = −   66 C10  π  1 1 Dmn a b Las constante de integración son determinadas considerando la continuidad de los esfuerzos cortantes entre las capas y estas son desaparecen en las caras inferior y superior. En el caso de una carga uniforme, donde p(x,y)=const=p0 la ecuación C3 es: p 16 p0mn = 2 si m y n son imparesπ m n C11 pmn = 0 si m y n son pares Esas dos expresiones puede ser reagrupadas en la siguiente forma: p 16 p0 12m−1,2n−1 = 2 , m, n = 1, 2,K π (2m −1)(2n −1) C12 La deflexión C6 en el punto (x,y) es entonces: sen (2m 1) π x16 p a4 m=∞ n=∞ − sen (2n π y −1) w( x, y) = 0 ∑ ∑ a b C13 π 4 1 1 (2m −1) (2n −1) D2m−1, 2n−1 o 4 m=∞ n=∞ w( x, y) 16 p= 0a4 ∑ ∑cmn sen (2m −1) π x sen (2n −1) π y C14 π 1 1 a b con c 1mn = (2m −1)(2n −1) D C15 2m−1,2n−1 Los esfuerzos en el plano en la capa k son:  ∂2w   − ∂x2   σ  E E  x 11 12 0      2   σ y  = z E E 0  ∂ w 12 22   − 2  C16  ∂y  σ x y  k  0 0 E66 k  ∂2w  −2   ∂x ∂y  donde ∂2w 16 p a2 m=∞ n=∞0 ∑ ∑ (2m 1)2 c sen (2m 1) π x sen (2n 1) π y− ∂x2 = π 4 − mn − −1 1 a b ∂2w 16 p 2 2 m=∞ n=∞ − = 0 a R 2 4 ∑ ∑ (2n −1) 2 cm n sen (2m π x −1) sen (2n 1) π y− ∂y π 1 1 a b C17 2 2 ∂ w− = ∂x∂y 32 p a2R m=∞ n=∞0 ∑ ∑ (2m 1)(2n 1)c sen (2m 1) π x π y− 4 − − mn − sen (2n −1)π 1 1 a b La deflexión es máxima en el centro de la placa, x=a/2, y=b/2. los numeradores en la ecuación C13 son entonces: sen (2m 1) π− sen (2n 1) π− = (−1)m+n−2 2 2 y la máxima deflexión es: w 16 p0a 4 max = 6 α C18 π donde m=∞ n=∞ α = ∑ ∑αmn 1 1 m+n−2 C19 α = (−1)m+n−2 ( −1 mn c ) mn = (2m −1)(2n −1) D2m+1,2n+1 De la ecuación C16 los esfuerzos en el centro de la placa son expresados como:  σ x  E11 E12 0   β   σ  z E E 0  16 p0a 2 =  2  y   12 22  π 4  R γ  C20 σ x y   0 0 E66   0 k k donde m=∞ n=∞ β = ∑ ∑ βmn β 2m n = (2m −1) αmn C21 1 1 m=∞ n=∞ γ = ∑ ∑γ 2mn γ mn = (2n −1) αmn C22 1 1 con 16 p 2 σ k = 0 a k x x 4 (E11β + Ek R2γ ) zπ 12 k 16 p0a 2 σ y y = (Ek β + Ek R2γ ) z C23 π 4 12 22 σ kx y = 0 o también σ k = Akx x x x z σ ky y = A k y y z C24 σ kx y = 0 y 16 p a2Ak = 0 (Ek β + E k R2x x π 4 11 12 γ ) 16 p a2 C25 Ak = 0 k k 2y y π 4 (E12β + E22R γ ) ANEXO D. CODIGO DE PROGRAMACION EN MATLAB function [Bb] = bend(eo, no, x1,x2,x3,x4,y1,y2,y3,y4,J,J11,J12,J21,J22); %Calculo de la transformacion invJ=inv(J); Jo=[J11^2, J12^2, 2*J11*J12; J21^2, J22^2, 2*J21*J22; J11*J21, J12*J22, J21*J12+J11*J22]; JJ=inv(Jo); JM=[0 0;0 0;0.25*(x1-x2)+0.25*(x3-x4), 0.25*(y1-y2)+0.25*(y3-y4)]; Bb=zeros(3,16); ee=[-1 1 1 -1];, nn=[-1 -1 1 1]; for ii = 1 : 4 ei=ee(ii);, ni=nn(ii); %Segundas derivadas dN/de2 ddNie =[1/8*(eo*ei-2)*(no+ni)^2*(no*ni- 2)+1/4*(eo+ei)*ei*(no+ni)^2*(no*ni-2),... -1/8*ei*(eo*ei-1)*(no+ni)^2*(no*ni-2)- 1/4*ei^2*(eo+ei)*(no+ni)^2*(no*ni-2),... -1/8*(eo*ei-2)*ni*(no+ni)^2*(no*ni-1)- 1/4*(eo+ei)*ei*ni*(no+ni)^2*(no*ni-1),... 1/8*ei*(eo*ei-1)*ni*(no+ni)^2*(no*ni- 1)+1/4*ei^2*(eo+ei)*ni*(no+ni)^2*(no*ni-1)]; %Segundas derivadas dN/dn2 ddNin =[1/8*(eo+ei)^2*(eo*ei-2)*(no*ni-2)+1/4*(eo+ei)^2*(eo*ei- 2)*(no+ni)*ni,... -1/8*ei*(eo+ei)^2*(eo*ei-1)*(no*ni-2)- 1/4*ei*(eo+ei)^2*(eo*ei-1)*(no+ni)*ni,... -1/8*(eo+ei)^2*(eo*ei-2)*ni*(no*ni-1)- 1/4*(eo+ei)^2*(eo*ei-2)*ni^2*(no+ni),... 1/8*ei*(eo+ei)^2*(eo*ei-1)*ni*(no*ni- 1)+1/4*ei*(eo+ei)^2*(eo*ei-1)*ni^2*(no+ni)]; %Segundas derivadas dN/dnde ddNien =[1/4*(eo+ei)*(eo*ei-2)*(no+ni)*(no*ni- 2)+1/8*(eo+ei)*(eo*ei-2)*(no+ni)^2*ni+1/8*(eo+ei)^2*ei*(no+ni)*(no*ni- 2)+1/16*(eo+ei)^2*ei*(no+ni)^2*ni,... -1/4*ei*(eo+ei)*(eo*ei-1)*(no+ni)*(no*ni-2)- 1/8*ei*(eo+ei)*(eo*ei-1)*(no+ni)^2*ni-1/8*ei^2*(eo+ei)^2*(no+ni)*(no*ni- 2)-1/16*ei^2*(eo+ei)^2*(no+ni)^2*ni,... -1/4*(eo+ei)*(eo*ei-2)*ni*(no+ni)*(no*ni-1)- 1/8*(eo+ei)*(eo*ei-2)*ni^2*(no+ni)^2-1/8*(eo+ei)^2*ei*ni*(no+ni)*(no*ni- 1)-1/16*(eo+ei)^2*ei*ni^2*(no+ni)^2,... 1/4*ei*(eo+ei)*(eo*ei-1)*ni*(no+ni)*(no*ni- 1)+1/8*ei*(eo+ei)*(eo*ei- 1)*ni^2*(no+ni)^2+1/8*ei^2*(eo+ei)^2*ni*(no+ni)*(no*ni- 1)+1/16*ei^2*(eo+ei)^2*ni^2*(no+ni)^2]; %Primera derivadas dNi/de dNie = [1/8*(eo+ei)*(eo*ei-2)*(no+ni)^2*(no*ni- 2)+1/16*(eo+ei)^2*ei*(no+ni)^2*(no*ni-2),... -1/8*ei*(eo+ei)*(eo*ei-1)*(no+ni)^2*(no*ni-2)- 1/16*ei^2*(eo+ei)^2*(no+ni)^2*(no*ni-2),... -1/8*(eo+ei)*(eo*ei-2)*ni*(no+ni)^2*(no*ni-1)- 1/16*(eo+ei)^2*ei*ni*(no+ni)^2*(no*ni-1),... 1/8*ei*(eo+ei)*(eo*ei-1)*ni*(no+ni)^2*(no*ni- 1)+1/16*ei^2*(eo+ei)^2*ni*(no+ni)^2*(no*ni-1)]; %Primera derivadas dNi/dn dNin = [1/8*(eo+ei)^2*(eo*ei-2)*(no+ni)*(no*ni- 2)+1/16*(eo+ei)^2*(eo*ei-2)*ni*(no+ni)^2,... -1/8*ei*(eo+ei)^2*(eo*ei-1)*(no+ni)*(no*ni-2)- 1/16*ei*(eo+ei)^2*(eo*ei-1)*ni*(no+ni)^2,... -1/8*(eo+ei)^2*(eo*ei-2)*ni*(no+ni)*(no*ni-1)- 1/16*(eo+ei)^2*(eo*ei-2)*ni^2*(no+ni)^2,... 1/8*ei*(eo+ei)^2*(eo*ei-1)*ni*(no+ni)*(no*ni- 1)+1/16*ei*(eo+ei)^2*(eo*ei-1)*ni^2*(no+ni)^2]; %Calculo de las primeras derivadas de Ni respecto a x, y dNi_1=invJ*[dNie(1); dNin(1)];, dNi_2=invJ*[dNie(2); dNin(2)]; dNi_3=invJ*[dNie(3); dNin(3)];, dNi_4=invJ*[dNie(4); dNin(4)]; %Calculo de las segundas derivadas de Ni respecto a e y n ddNi2=[ddNie;ddNin;ddNien]; %Calculo de la matrix Bb Bb(:,ii)=JJ*(ddNi2(:,1)-(JM*dNi_1)); Bb(:,ii+4)=JJ*(ddNi2(:,2)-(JM*dNi_2)); Bb(:,ii+8)=JJ*(ddNi2(:,3)-(JM*dNi_3)); Bb(:,ii+12)=JJ*(ddNi2(:,4)-(JM*dNi_4)); end function [D] = elas(E1x, E2x, v1x, v2x, Gx, tx, anl) DL=zeros(6);, E=zeros(6);, D=DL; nc=length(E1x); if length(tx)==2 % No laminado for k=1:nc al=anl(k); E1=E1x(k);, E2=E2x(k);, v1=v1x(k);, v2=v2x(k);, G=Gx(k);, t=tx; D11=(E1/(1-v1*v2)); D12=(v1*E2/(1-v1*v2)); D22=(E2/(1-v1*v2)); D66=G; Dx=[D11,D12,0;D12,D22,0;0,0,D66]; %Calculo de la matriz elastica ortotropica D11o=D11*(cos(al))^4+2*(D12+2*D66)*(sin(al))^2*(cos(al))^2+D22*(sin(al))^ 4; D12o=(D11+D22- 4*D66)*(sin(al))^2*(cos(al))^2+D12*((sin(al))^4+(cos(al))^4); D13o=(D11-D12-2*D66)*(sin(al))*(cos(al))^3+(D12- D22+2*D66)*(sin(al))^3*(cos(al)); D22o=D11*(sin(al))^4+2*(D12+2*D66)*(sin(al))^2*(cos(al))^2+D22*(cos(al))^ 4; D23o=(D11-D12-2*D66)*(sin(al))^3*(cos(al))+(D12- D22+2*D66)*(sin(al))*(cos(al))^3; D33o=(D11+D22-2*D12- 2*D66)*(sin(al))^2*(cos(al))^2+D66*((sin(al))^4+(cos(al))^4); Do=[D11o,D12o,D13o;D12o,D22o,D23o;D13o,D23o,D33o]; Ad=Do*(t(k)-t(k+1));, Bd=zeros(3);, Dd=(1/12)*Do*((t(k))- (t(k+1)))^3; DL=[Ad Bd; Bd Dd]; D=D+DL; end else for k=1:nc % Laminado al=anl(k); E1=E1x(k);, E2=E2x(k);, v1=v1x(k);, v2=v2x(k);, G=Gx(k);, t=tx; D11=(E1/(1-v1*v2)); D12=(v1*E2/(1-v1*v2)); D22=(E2/(1-v1*v2)); D66=G; Dx=[D11,D12,0;D12,D22,0;0,0,D66]; %Calculo de la matriz elastica ortotropica D11o=D11*(cos(al))^4+2*(D12+2*D66)*(sin(al))^2*(cos(al))^2+D22*(sin(al))^ 4; D12o=(D11+D22- 4*D66)*(sin(al))^2*(cos(al))^2+D12*((sin(al))^4+(cos(al))^4); D13o=(D11-D12-2*D66)*(sin(al))*(cos(al))^3+(D12- D22+2*D66)*(sin(al))^3*(cos(al)); D22o=D11*(sin(al))^4+2*(D12+2*D66)*(sin(al))^2*(cos(al))^2+D22*(cos(al))^ 4; D23o=(D11-D12-2*D66)*(sin(al))^3*(cos(al))+(D12- D22+2*D66)*(sin(al))*(cos(al))^3; D33o=(D11+D22-2*D12- 2*D66)*(sin(al))^2*(cos(al))^2+D66*((sin(al))^4+(cos(al))^4); Do=[D11o,D12o,D13o;D12o,D22o,D23o;D13o,D23o,D33o]; Ad=Do*(t(k)-t(k+1));, Bd=0.5*Do*((t(k))^2-(t(k+1))^2);, Dd=(1/3)*Do*((t(k))^3-(t(k+1))^3); DL=[Ad Bd; Bd Dd]; D=D+DL; end end function [E, z] = elas_E(E1x, E2x, v1x, v2x, Gx, tx, anl) k=1; %Capa 1 z = tx(k); %Distancia del plano medio a una capa DL=zeros(6);, E=zeros(6);, D=DL; al=anl(k); E1=E1x(k);, E2=E2x(k);, v1=v1x(k);, v2=v2x(k);, G=Gx(k);, t=tx; D11=(E1/(1-v1*v2)); D12=(v1*E2/(1-v1*v2)); D22=(E2/(1-v1*v2)); D66=G; Dx=[D11,D12,0;D12,D22,0;0,0,D66]; %Calculo de la matriz elastica ortotropica D11o=D11*(cos(al))^4+2*(D12+2*D66)*(sin(al))^2*(cos(al))^2+D22*(sin(al))^ 4; D12o=(D11+D22- 4*D66)*(sin(al))^2*(cos(al))^2+D12*((sin(al))^4+(cos(al))^4); D13o=(D11-D12-2*D66)*(sin(al))*(cos(al))^3+(D12- D22+2*D66)*(sin(al))^3*(cos(al)); D22o=D11*(sin(al))^4+2*(D12+2*D66)*(sin(al))^2*(cos(al))^2+D22*(cos(al))^ 4; D23o=(D11-D12-2*D66)*(sin(al))^3*(cos(al))+(D12- D22+2*D66)*(sin(al))*(cos(al))^3; D33o=(D11+D22-2*D12- 2*D66)*(sin(al))^2*(cos(al))^2+D66*((sin(al))^4+(cos(al))^4); Do=[D11o,D12o,D13o;D12o,D22o,D23o;D13o,D23o,D33o]; Ad=Do*(t(k)-t(k+1));, Bd=zeros(3);, Dd=(1/12)*Do*((t(k))-(t(k+1)))^3; E=[Do Bd; Bd Do]; function [Ke] = Kelem(xyz_nodos, NNT, x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4, E1x, E2x, v1x, v2x, tx, Gx, anl) Ke = zeros(24);, B = zeros(6,24); T = trans_matrix(xyz_nodos,NNT,x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4); xy_2 = trans_coord(xyz_nodos,NNT,x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4); %Transformacion de coordenadas x1=xy_2(1,1);, x2=xy_2(1,2);, x3=xy_2(1,3);, x4=xy_2(1,4); y1=xy_2(2,1);, y2=xy_2(2,2);, y3=xy_2(2,3);, y4=xy_2(2,4); [D] = elas(E1x, E2x, v1x, v2x, Gx, tx, anl); %Matriz de elasticidad eg = [-0.5773502692; 0.5773502692]; ng = [-0.5773502692; 0.5773502692]; for i = 1:2 for j = 1:2 eo = eg(i); no = ng(j); J11 = -0.25*(1-no)*x1+0.25*(1-no)*x2+0.25*(1+no)*x3- 0.25*(1+no)*x4; %dx/de J12 = -0.25*(1-no)*y1+0.25*(1-no)*y2+0.25*(1+no)*y3- 0.25*(1+no)*y4; %dy/de J21 = -0.25*(1-eo)*x1-0.25*(1+eo)*x2+0.25*(1+eo)*x3+0.25*(1- eo)*x4; %dx/dn J22 = -0.25*(1-eo)*y1-0.25*(1+eo)*y2+0.25*(1+eo)*y3+0.25*(1- eo)*y4; %dy/dn J = [J11 J12; J21 J22]; Bp = plane(eo, no, x1,x2,x3,x4,y1,y2,y3,y4,J); Bb = bend(eo, no, x1,x2,x3,x4,y1,y2,y3,y4,J,J11,J12,J21,J22); B(1:3,1:8) = Bp; B(4:6,9:24) = Bb; Ke = Ke+B'*D*B*det(J); end end Ke=T'*Ke*T; function [M]=mass(xyz_nodos, conx_data, NNT, Nel, tx, d, conex) load Uknow M=zeros(6*NNT); t=2*tx(1); % Espesor total de la lamina eg=[-0.5773502692; 0.5773502692]; ng=[-0.5773502692; 0.5773502692]; for h=1:Nel x1=xyz_nodos(conx_data(h,2),2);, y1=xyz_nodos(conx_data(h,2),3);, z1=xyz_nodos(conx_data(h,2),4); x2=xyz_nodos(conx_data(h,3),2);, y2=xyz_nodos(conx_data(h,3),3);, z2=xyz_nodos(conx_data(h,3),4); x3=xyz_nodos(conx_data(h,4),2);, y3=xyz_nodos(conx_data(h,4),3);, z3=xyz_nodos(conx_data(h,4),4); x4=xyz_nodos(conx_data(h,5),2);, y4=xyz_nodos(conx_data(h,5),3);, z4=xyz_nodos(conx_data(h,5),4); T = trans_matrix(xyz_nodos,NNT,x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4); % Transformacion de matrices xy_2 = trans_coord(xyz_nodos,NNT,x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4); %Transformacion de coordenadas x1=xy_2(1,1);, x2=xy_2(1,2);, x3=xy_2(1,3);, x4=xy_2(1,4); y1=xy_2(2,1);, y2=xy_2(2,2);, y3=xy_2(2,3);, y4=xy_2(2,4); Me = zeros(24); for i=1:2 for j=1:2 eo = eg(i); no = ng(j); Mp = massp(x1, x2, x3, x4, y1, y2, y3, y4, t, d, eo, no); Mb = massb(x1, x2, x3, x4, y1, y2, y3, y4, t, d, eo, no); Meo=zeros(24); Meo(1:8,1:8)=Mp; Meo(9:24,9:24)=Mb; Me = Me + Meo; end end Me=T'*Me*T; for m=1:24 for n=1:24 M(conex(h,m),conex(h,n))=M(conex(h,m),conex(h,n))+Me(m,n); end end end % Reducción de la matriz de masa load Uknow c = 1; cmax = 6*NNT; while c <= cmax if Uknow(c) == 0 M(:,c) = []; M(c,:) = []; Uknow(c) = []; c = 0; cmax = cmax - 1; end; c = c + 1; end function [Mb]=massb(x1, x2, x3, x4, y1, y2, y3, y4, t, d, eo, no) Mb=zeros(16); N1=zeros(1,4);, N2=N1;, N3=N1;, N4=N1; eh=[-1 1 1 -1];, nh=[-1 -1 1 1]; for j=1:4 ei=eh(j); ni=nh(j); Nw=(1/16)*((eo+ei)^2)*(eo*ei-2)*((no+ni)^2)*(no*ni-2); Nx=-(1/16)*ei*((eo+ei)^2)*(eo*ei-1)*((no+ni)^2)*(no*ni-2); Ny=-(1/16)*((eo+ei)^2)*(eo*ei-2)*ni*((no+ni)^2)*(no*ni-1); Nxy=(1/16)*ei*((eo+ei)^2)*(eo*ei-1)*ni*((no+ni)^2)*(no*ni-1); N1(j)=Nw; N2(j)=Nx; N3(j)=Ny; N4(j)=Nxy; end N = [N1 N2 N3 N4]; J11 = -0.25*(1-no)*x1+0.25*(1-no)*x2+0.25*(1+no)*x3-0.25*(1+no)*x4; J12 = -0.25*(1-no)*y1+0.25*(1-no)*y2+0.25*(1+no)*y3-0.25*(1+no)*y4; J21 = -0.25*(1-eo)*x1-0.25*(1+eo)*x2+0.25*(1+eo)*x3+0.25*(1-eo)*x4; J22 = -0.25*(1-eo)*y1-0.25*(1+eo)*y2+0.25*(1+eo)*y3+0.25*(1-eo)*y4; J=[J11 J12; J21 J22]; Mb = t*d*N'*N*det(J); function [Mp] = massp(x1, x2, x3, x4, y1, y2, y3, y4, t, d, eo, no) Mp=zeros(8); N1=(0.25)*(1-eo)*(1-no); N2=(0.25)*(1+eo)*(1-no); N3=(0.25)*(1+eo)*(1+no); N4=(0.25)*(1-eo)*(1+no); N=[N1 N2 N3 N4 0 0 0 0; 0 0 0 0 N1 N2 N3 N4]; J11 = -0.25*(1-no)*x1+0.25*(1-no)*x2+0.25*(1+no)*x3-0.25*(1+no)*x4; J12 = -0.25*(1-no)*y1+0.25*(1-no)*y2+0.25*(1+no)*y3-0.25*(1+no)*y4; J21 = -0.25*(1-eo)*x1-0.25*(1+eo)*x2+0.25*(1+eo)*x3+0.25*(1-eo)*x4; J22 = -0.25*(1-eo)*y1-0.25*(1+eo)*y2+0.25*(1+eo)*y3+0.25*(1-eo)*y4; J=[J11 J12; J21 J22]; Mp = t*d*N'*N*det(J); function [K]=matrix(conx_data, xyz_nodos, Nel, NNT, E1x, E2x, v1x, v2x, tx, Gx, anl, conex) K=zeros(6*NNT); for h=1:Nel x1=xyz_nodos(conx_data(h,2),2);, y1=xyz_nodos(conx_data(h,2),3);, z1=xyz_nodos(conx_data(h,2),4); x2=xyz_nodos(conx_data(h,3),2);, y2=xyz_nodos(conx_data(h,3),3);, z2=xyz_nodos(conx_data(h,3),4); x3=xyz_nodos(conx_data(h,4),2);, y3=xyz_nodos(conx_data(h,4),3);, z3=xyz_nodos(conx_data(h,4),4); x4=xyz_nodos(conx_data(h,5),2);, y4=xyz_nodos(conx_data(h,5),3);, z4=xyz_nodos(conx_data(h,5),4); Ke = Kelem(xyz_nodos, NNT, x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4, E1x, E2x, v1x, v2x, tx, Gx, anl); for m=1:24 for n=1:24 K(conex(h,m),conex(h,n))=K(conex(h,m),conex(h,n))+Ke(m,n); end end end function [Bp] = plane(eo, no, x1,x2,x3,x4,y1,y2,y3,y4,J) Bp=zeros(3,8); %calculo de las derivadas de Ni/e,n DN1c = [-0.25*(1-no);-0.25*(1-eo)]; DN2c = [0.25*(1-no);-0.25*(1+eo)]; DN3c = [0.25*(1+no);0.25*(1+eo)]; DN4c = [-0.25*(1+no);0.25*(1-eo)]; DN1x = inv(J)*DN1c; DN2x = inv(J)*DN2c; DN3x = inv(J)*DN3c; DN4x = inv(J)*DN4c; %Calculo de la matriz Bp Bp=[DN1x(1) DN2x(1) DN3x(1) DN4x(1) 0 0 0 0;... 0 0 0 0 DN1x(2) DN2x(2) DN3x(2) DN4x(2);... DN1x(2) DN2x(2) DN3x(2) DN4x(2) DN1x(1) DN2x(1) DN3x(1) DN4x(1)]; function [S]=stress(conx_data, NNT, xyz_nodos, conex, E1x, E2x, v1x, v2x, Gx, tx, anl, Nel) load U %Transformacion de U a coordenadas locales Tg=zeros(6*NNT); for h=1:Nel x1=xyz_nodos(conx_data(h,2),2);, y1=xyz_nodos(conx_data(h,2),3);, z1=xyz_nodos(conx_data(h,2),4); x2=xyz_nodos(conx_data(h,3),2);, y2=xyz_nodos(conx_data(h,3),3);, z2=xyz_nodos(conx_data(h,3),4); x3=xyz_nodos(conx_data(h,4),2);, y3=xyz_nodos(conx_data(h,4),3);, z3=xyz_nodos(conx_data(h,4),4); x4=xyz_nodos(conx_data(h,5),2);, y4=xyz_nodos(conx_data(h,5),3);, z4=xyz_nodos(conx_data(h,5),4); T = trans_matrix(xyz_nodos,NNT,x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4); for m = 1:24 for n = 1:24 Tg(conex(h,m),conex(h,n)) = T(m,n); end end end U=inv(Tg)*U; %Matriz de resultandos sols=zeros(NNT,6); sol(:,1)=U(1:NNT'); sols(:,2)=U(NNT+1:2*NNT'); sols(:,3)=U(2*NNT+1:3*NNT'); sols(:,4)=U(3*NNT+1:4*NNT'); sols(:,5)=U(4*NNT+1:5*NNT'); sols(:,6)=U(5*NNT+1:6*NNT'); %Resultados para cada elemento F1 = sol(:,1);, F2 = sol(:,2);, F3 = sol(:,3);, F4 = sol(:,4);, F5 = sol(:,5);, F6 = sol(:,6); Ux(1:4,conx_data(:,1)) = F1(conx_data(:,2:5)'); Uy(1:4,conx_data(:,1)) = F2(conx_data(:,2:5)'); Uz(1:4,conx_data(:,1)) = F3(conx_data(:,2:5)'); Ox(1:4,conx_data(:,1)) = F4(conx_data(:,2:5)'); Oy(1:4,conx_data(:,1)) = F5(conx_data(:,2:5)'); Oxy(1:4,conx_data(:,1)) = F6(conx_data(:,2:5)'); Ue=[Ux; Uy; Uz; Ox; Oy; Oxy]; %Calculo de esfuerzos [E, z] = elas_E(E1x, E2x, v1x, v2x, Gx, tx, anl); ee = [-0.5773502692; 0.5773502692; 0.5773502692; -0.5773502692]; nn = [-0.5773502692; -0.5773502692; 0.5773502692; 0.5773502692]; S = zeros(3,Nel); for h=1:Nel r=1;, s=3; xy_2 = trans_coord(xyz_nodos,NNT,x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4); x1=xy_2(1,1);, x2=xy_2(1,2);, x3=xy_2(1,3);, x4=xy_2(1,4); y1=xy_2(2,1);, y2=xy_2(2,2);, y3=xy_2(2,3);, y4=xy_2(2,4); for i=1:4 eo=ee(i);, no=nn(i); J11 = -0.25*(1-no)*x1+0.25*(1-no)*x2+0.25*(1+no)*x3- 0.25*(1+no)*x4; %dx/de J12 = -0.25*(1-no)*y1+0.25*(1-no)*y2+0.25*(1+no)*y3- 0.25*(1+no)*y4; %dy/de J21 = -0.25*(1-eo)*x1-0.25*(1+eo)*x2+0.25*(1+eo)*x3+0.25*(1- eo)*x4; %dx/dn J22 = -0.25*(1-eo)*y1-0.25*(1+eo)*y2+0.25*(1+eo)*y3+0.25*(1- eo)*y4; %dy/dn J = [J11 J12; J21 J22]; B = zeros(6,24); Bp = plane(eo, no, x1,x2,x3,x4,y1,y2,y3,y4,J); Bb = bend(eo, no, x1,x2,x3,x4,y1,y2,y3,y4,J,J11,J12,J21,J22); B(1:3,1:8) = Bp; B(4:6,9:24) = Bb; Sn = E*B*Ue(:,h); So = Sn(1:3)+z*Sn(4:6); %Esfuerzos en el nodo i S(r:s,h) = So; r = r+3;, s = s+3; end end function [St]=stress_t(conx_data, NNT, xyz_nodos, conex, E1x, E2x, v1x, v2x, Gx, tx, anl, Nel) load Um load Uknow %Construyendo los vectores de solución [m,n,p]=size(Um); j = 0; for k=1:p for i = 1 : 6*NNT if Uknow(i) == 1 j = j + 1; U(i,1,k) = Um(j,1,k); end end j=0; end %Transformacion de U a coordenadas locales Tg=zeros(6*NNT); for h=1:Nel x1=xyz_nodos(conx_data(h,2),2);, y1=xyz_nodos(conx_data(h,2),3);, z1=xyz_nodos(conx_data(h,2),4); x2=xyz_nodos(conx_data(h,3),2);, y2=xyz_nodos(conx_data(h,3),3);, z2=xyz_nodos(conx_data(h,3),4); x3=xyz_nodos(conx_data(h,4),2);, y3=xyz_nodos(conx_data(h,4),3);, z3=xyz_nodos(conx_data(h,4),4); x4=xyz_nodos(conx_data(h,5),2);, y4=xyz_nodos(conx_data(h,5),3);, z4=xyz_nodos(conx_data(h,5),4); T = trans_matrix(xyz_nodos,NNT,x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4); for m = 1:24 for n = 1:24 Tg(conex(h,m),conex(h,n)) = T(m,n); end end end for i=1:p U(:,:,i)=inv(Tg)*U(:,:,i); end %Matriz de resultandos U_x=U(1:NNT,1,:); U_y=U(NNT+1:2*NNT,1,:); U_z=U(2*NNT+1:3*NNT,1,:); U_ox=U(3*NNT+1:4*NNT,1,:); U_oy=U(4*NNT+1:5*NNT,1,:); U_oxy=U(5*NNT+1:6*NNT,1,:); St=zeros(12,Nel,p); for k=1:p F1 = U_x(:,:,1);, F2 = U_y(:,:,i);, F3 = U_z(:,:,i);, F4 = U_ox(:,:,i);, F5 = U_oy(:,:,i);, F6 = U_oxy(:,:,i); Ux(1:4,conx_data(:,1)) = F1(conx_data(:,2:5)'); Uy(1:4,conx_data(:,1)) = F2(conx_data(:,2:5)'); Uz(1:4,conx_data(:,1)) = F3(conx_data(:,2:5)'); Ox(1:4,conx_data(:,1)) = F4(conx_data(:,2:5)'); Oy(1:4,conx_data(:,1)) = F5(conx_data(:,2:5)'); Oxy(1:4,conx_data(:,1)) = F6(conx_data(:,2:5)'); Ue=[Ux; Uy; Uz; Ox; Oy; Oxy]; %Calculo de esfuerzos [E, z] = elas_E(E1x, E2x, v1x, v2x, Gx, tx, anl); ee = [-0.5773502692; 0.5773502692; 0.5773502692; -0.5773502692]; nn = [-0.5773502692; -0.5773502692; 0.5773502692; 0.5773502692]; S = zeros(3,Nel); for h=1:Nel r=1;, s=3; xy_2 = trans_coord(xyz_nodos,NNT,x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4); x1=xy_2(1,1);, x2=xy_2(1,2);, x3=xy_2(1,3);, x4=xy_2(1,4); y1=xy_2(2,1);, y2=xy_2(2,2);, y3=xy_2(2,3);, y4=xy_2(2,4); for i=1:4 eo=ee(i);, no=nn(i); J11 = -0.25*(1-no)*x1+0.25*(1-no)*x2+0.25*(1+no)*x3- 0.25*(1+no)*x4; %dx/de J12 = -0.25*(1-no)*y1+0.25*(1-no)*y2+0.25*(1+no)*y3- 0.25*(1+no)*y4; %dy/de J21 = -0.25*(1-eo)*x1-0.25*(1+eo)*x2+0.25*(1+eo)*x3+0.25*(1- eo)*x4; %dx/dn J22 = -0.25*(1-eo)*y1-0.25*(1+eo)*y2+0.25*(1+eo)*y3+0.25*(1- eo)*y4; %dy/dn J = [J11 J12; J21 J22]; B = zeros(6,24); Bp = plane(eo, no, x1,x2,x3,x4,y1,y2,y3,y4,J); Bb = bend(eo, no, x1,x2,x3,x4,y1,y2,y3,y4,J,J11,J12,J21,J22); B(1:3,1:8) = Bp; B(4:6,9:24) = Bb; Sn = E*B*Ue(:,h); So = Sn(1:3)+z*Sn(4:6); %Esfuerzos en el nodo i S(r:s,h) = So; r = r+3;, s = s+3; end end k St(:,:,k)=S; end function time(xyz_nodos, conx_data, NNT, Nel, tx, d, conex) load Kr load desp U = desp; % Vector de Desplazamiento disp('Calculando la matriz de masa ...') M=mass(xyz_nodos, conx_data, NNT, Nel, tx, d, conex); % Matriz de Masa Ace = 9.81*ones(size(U)); % Vector de Aceleracion Vel = zeros(size(U));, Vel(1:651)=4.8522; F= zeros(size(U)); a = 0.01;, b = 0.0008; % Coeficientes de Rayleigh C = a*M + b*Kr; % Matriz de amortiguamiento t=0.2; DT=4e-4; p1=0.25;, p2=0.5; % Coeficientes de Newmark a0=1/(p1*DT^2);, a1=p2/(p1*DT);, a2=1/(p1*DT);, a3=(1/(2*p1))-1; a4=(p2/p1) -1;, a5=(DT/2)*((p2/p1)-2);, a6=DT*(1-p2);, a7=p2*DT; %Archivo de pelicula p=length(0:DT:t);, [m,n]=size(U); Um=zeros(m,n,p); disp('Resolviendos los sistemas de ecuaciones ...') Ko = inv(a0*M+a1*C+Kr); pd=145; j=1; for i=0:DT:t Um(:,:,j)=U; ((t/DT)+1)-j U2(j)=U(pd);, Vel2(j)=Vel(pd);, Ace2(j)=Ace(pd); Mo = M*((a0*U)+(a2*Vel)+(a3*Ace));, Co = C*(a1*U+a4*Vel+a5*Ace); Fo = F+Mo+Co; Ut=Ko*Fo; Ace1=a0*(Ut-U)-a2*Vel-a3*Ace; Vel1=Vel+a6*Ace+a7*Ace1; Vel=Vel1;, Ace=Ace1;, U=Ut; j=j+1; end disp('Guardando datos ...') save t save DT save U2 save Vel2 save Ace2 save Um disp('Calculo finalizado ...') function [xy_2] = trans_coord(xyz_nodos,NNT,x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4) Xo = [0.25*(x1+x2+x3+x4) 0.25*(y1+y2+y3+y4) 0.25*(z1+z2+z3+z4)]; %Vector de posicion Xo Vij= [(x2-x1) (y2-y1) (z2-z1)]; %Vector Vij (direccion ij) Vim= [(x4-x1) (y4-y1) (z4-z1)]; %Vector im (direccion im) Vzo= cross(Vij,Vim); %Vector unitario Vz' Vz = Vzo/norm(Vzo); Vxo = [(0.5*(x3+x2)-Xo(1)) (0.5*(y3+y2)-Xo(2)) (0.5*(z3+z2)-Xo(3))]; %Vector C-nj Vx = Vxo/norm(Vxo); %Vector unitario Vx' Vy = cross(Vz,Vx); %Vector unitario Vy' A = [Vx;Vy;Vz]; %Se transforman a coordenadas nodales xy_2(:,1)=A*[x1 y1 z1]'; xy_2(:,2)=A*[x2 y2 z2]'; xy_2(:,3)=A*[x3 y3 z3]'; xy_2(:,4)=A*[x4 y4 z4]'; function [T]=trans_matrix(xyz_nodos,NNT,x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4) Tsup=zeros(12); T=zeros(24); Xo = [0.25*(x1+x2+x3+x4) 0.25*(y1+y2+y3+y4) 0.25*(z1+z2+z3+z4)]; %Vector de posicion Xo Vij= [(x2-x1) (y2-y1) (z2-z1)]; %Vector Vij (direccion ij) Vim= [(x4-x1) (y4-y1) (z4-z1)]; %Vector im (direccion im) Vzo= cross(Vij,Vim); %Vector unitario Vz' Vz = Vzo/norm(Vzo); Vxo = [(0.5*(x3+x2)-Xo(1)) (0.5*(y3+y2)-Xo(2)) (0.5*(z3+z2)-Xo(3))]; %Vector C-nj Vx = Vxo/norm(Vxo); %Vector unitario Vx' Vy = cross(Vz,Vx); %Vector unitario Vy' A = [Vx;Vy;Vz]; for i = 1:3 %T es la matriz de transformacion global k=4*i-3; for m =0:3 for j = 1:4:9 Tsup(k+m,j+m)=A(i,(j+3)/4); end end end T(1:12,1:12)=Tsup; T(13:24,13:24)=Tsup;