Repositorio Institucional de Documentos

Universidad de Zaragoza Escuela de Ingenier´ ıa y Arquitectura Trabajo Fin de Master Master de Mec´ anica Aplicada Programa oficial de Posgrado en Ingenier´ ıa Mec´ anica y materiales ESTIMACI´ ON DEL ERROR A POSTERIORI DEL M´ ETODO DE ELEMENTOS FINITOS CON LA TEOR´ IA VARIACIONAL DE LAS MULTIESCALAS APLICADO A ELASTICIDAD LINEAL Autor: Diego Irisarri Jim´enez Director: Guillermo Hauke Bernardos ´ Area de mec´anica de Fluidos Curso 2010/2011 Diciembre 2011 Estimaci´on del error a posteriori del m´etodo de elementos finitos con la teor´ıa variacional de las multiescalas aplicado a elasticidad lineal Resumen El presente trabajo investiga t´ecnicas de estimadores de error a posteriori para el m´etodo de elementos finitos aplicado al problema de elasticidad lineal. Para hallar el error que comete el m´etodo de elementos finitos se emplea el m´etodo variacional multiescala. Con la ayuda de este m´etodo y partiendo de los resultados obtenidos de la simulaci´on de elementos finitos, se obtiene una estimaci´on del error cometido. En concreto, en este proyecto se ha estudiado la estimaci´on de error para problemas de tensi´on plana. A la hora de analizar el acierto en la estimaci´on de error, es necesario conocer previamente el error real que comete el m´etodo de elementos finitos. Por ello, se ha calculado la soluci´on anal´ıtica de los casos estudiados de tensi´on plana, empleando las ecuaciones b´asicas de elasticidad. Los problemas estudiados de tensi´on plana son aplicados sobre dominios rectangulares. Sobre estos dominios se han impuesto distintos tipos de solicitaciones en el contorno para estudiar c´omo afectan los distintos tipos de esfuerzos en los estimadores de error. Tambi´en se analiza su influencia con los principales par´ametros que caracterizan el problema de elasticidad plana. Con el objetivo de cuantificar la calidad de los estimadores de error, se establece la denominada eficiencia del estimador de error, el cual relaciona el error estimado con el error realmente cometido. ´ Indice general 1. Introducci´on y objetivos 1 1.1. Objetivo del trabajo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2. Motivaci´on del trabajo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.3. Estimaci´on del error . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.4. Metodolog´ıa................................... 4 2. Elasticidad bidimensional 5 2.1. Estado tensional plano . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.2. Problemas de tensi´on plana . . . . . . . . . . . . . . . . . . . . . . . . . . 9 3. Problema el´astico lineal 11 3.1. Formulaci´on fuerte del problema . . . . . . . . . . . . . . . . . . . . . . . . 11 3.2. Formulaci´on d´ebil del problema . . . . . . . . . . . . . . . . . . . . . . . . 12 3.3. M´etodo de aproximaci´on de Galerkin . . . . . . . . . . . . . . . . . . . . . 14 3.4. Comparativa ecuaci´on de calor y elasticidad . . . . . . . . . . . . . . . . . 15 4. Estimaci´on de error a posteriori 16 4.1. M´etodo variacional de las multiescalas . . . . . . . . . . . . . . . . . . . . 16 4.2. Estimaci´on de error . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 4.3. Funci´onbubble ................................. 20 4.4. Expresi´on del estimador de error . . . . . . . . . . . . . . . . . . . . . . . . 22 4.5. Eficiencia global y eficiencia local . . . . . . . . . . . . . . . . . . . . . . . 23 5. Resultados 24 5.1. Problema de la ecuaci´on de calor . . . . . . . . . . . . . . . . . . . . . . . 24 5.2. Problemas sometidos a flexi´on en el contorno . . . . . . . . . . . . . . . . . 26 5.2.1. Soluci´on anal´ıtica de los problemas . . . . . . . . . . . . . . . . . . 27 5.2.2. Soluci´on de elementos finitos . . . . . . . . . . . . . . . . . . . . . . 28 5.2.3. An´alisis de la eficiencia del error . . . . . . . . . . . . . . . . . . . . 30 5.2.4. An´alisis de sensibilidad del estimador . . . . . . . . . . . . . . . . . 34 5.3. Viga en voladizo y viga biapoyada . . . . . . . . . . . . . . . . . . . . . . . 37 5.3.1. Soluci´on anal´ıtica de los problemas . . . . . . . . . . . . . . . . . . 37 5.3.2. Soluci´on por elementos finitos . . . . . . . . . . . . . . . . . . . . . 38 5.3.3. An´alisis de la eficiencia del error . . . . . . . . . . . . . . . . . . . . 38 6. Conclusiones 42 i ´ INDICE GENERAL ´ INDICE GENERAL Bibliograf´ıa 44 Anexos 45 A. Soluci´on anal´ıtica: Problema flexi´on Pura 47 B. Soluci´on anal´ıtica: Problema flexi´on compuesta 50 C. Soluci´on anal´ıtica: Problema de flexi´on en los dos ejes 52 D. Soluci´on anal´ıtica: Problema viga empotrada 55 E. Soluci´on anal´ıtica: Problema viga apoyada 59 ii Cap´ıtulo 1 Introducci´on y objetivos Los m´etodos computacionales son una de las ´ultimas herramientas a disposici´on del ingeniero a la hora de resolver modelos matem´aticos de un fen´omeno f´ısico. Se empezaron a desarrollar a mediados del siglo XX tras la irrupci´on de los computadores. En un principio se limitaron a apoyar o a extender resultados emp´ıricos. Hoy en d´ıa son capaces de ofrecer soluciones confiables a modelos f´ısicos plasmados en ecuaciones diferenciales y a pesar de su relativa novedad, se han integrado de una manera muy notable en los m´as variados procesos industriales. El r´apido desarrollo de los ordenadores ha supuesto un aumento de su capacidad de c´alculo al mismo tiempo que un abaratamiento de su coste. La mec´anica de medios continuos es un campo donde los m´etodos computacionales son de gran utilidad y al mismo tiempo suponen un gran reto. Ambas situaciones se deben a la gran complejidad de los modelos matem´aticos que describen su comportamiento. Hay que recordar que las soluciones anal´ıticas que existen en elasticidad son limitadas y sirven para geometr´ıas sencillas, o bien asumen hip´otesis simplificativas que desvirt´uan la soluci´on obtenida. Sin embargo, los m´etodos de elementos finitos se adaptan a cualquier geometr´ıa y permiten abordar el problema con toda su complejidad a cambio de tiempo de computaci´on. 1.1. Objetivo del trabajo El objetivo de este Trabajo Fin de Master es investigar la estimaci´on de error que se comete en la resoluci´on num´erica por el m´etodo de elementos finitos de las ecuaciones de elasticidad lineal. Para ello, se caracterizan y eval´uan estimadores de error para el m´etodo de elementos finitos basado en el M´etodo Variacional de las Multiescalas (VMS). La evaluaci´on consiste en la comparaci´on del error estimado con el error exacto. En el estudio de este trabajo se persiguen los siguientes objetivos concretos: Hallar un estimador de error adecuado para elasticidad lineal, en el que se consiga 1 CAP´ ITULO 1. INTRODUCCI´ ON Y OBJETIVOS 1.2. MOTIVACI´ ON DEL TRABAJO una precisi´on y robustez razonables Evaluar la influencia en la estimaci´on del error para los par´ametros caracter´ısticos de la elasticidad lineal. Partiendo de la formulaci´on d´ebil del problema y aplicando el VMS, se emplean estos estimadores de error, los cuales involucran directamente el residuo de la soluci´on de elementos finitos. De esta manera, se pretende proporcionar al usuario de programas de FEM una medida de la calidad de la soluci´on obtenida. 1.2. Motivaci´on del trabajo El prop´osito de realizar este trabajo, parti´o de la investigaci´on previamente realizada en el departamento de Mec´anica de Fluidos sobre estimaci´on de error ([10], [11], [12]) en la ecuaci´on del transporte. En particular, para ecuaciones diferenciales el´ıpticas, como la ecuaci´on del calor, los estimadores de error expl´ıcitos empleados proporcionan unos buenos resultados. Las ecuaciones diferenciales que gobiernan el problema el´astico lineal son el´ıpticas. Sin embargo, en relaci´on con la ecuaci´on de calor, la ecuaci´on de elasticidad presenta algunas diferencias significativas que a˜naden complejidad al problema y que se detallar´an m´as adelante. Desde este punto de partida y debido a la elevada utilizaci´on de programas de m´etodo de elementos finitos que existe hoy en d´ıa en el dise˜no mec´anico y estructural, mediante este trabajo, se estudia la estimaci´on de error en un campo en el que existe un gran inter´es en investigar el error que se comete en la soluci´on obtenida por el m´etodo de elementos finitos. 1.3. Estimaci´on del error La estimaci´on de error es uno de los campos de investigaci´on de mayor inter´es en la mec´anica computacional. Hoy en d´ıa, los m´etodos computacionales son fiables, ´utiles y cada vez m´as r´apidos, de tal manera que ya est´an plenamente integrados en el proceso productivo. La estimaci´on de error surge como una mejora natural de los mismos. Como en otros campos de la ciencia o de la ingenier´ıa, es deseable conocer con qu´e margen de error se ajustan a la realidad los c´alculos o incluso los resultados de un experimento. En los m´etodos experimentales, una medida no se considera completa sin una estimaci´on del error cometido. La tendencia en los m´etodos computacionales deber´ıa ser la misma, ya no s´olo por una cuesti´on de rigor cient´ıfico, sino porque una correcta evaluaci´on del error puede suponer un ahorro considerable en el coste computacional (y por tanto econ´omico) de las simulaciones. Las soluciones obtenidas por los m´etodos computacionales en general y el M´etodo de los elementos finitos, en particular, dependen fuertemente de la discretizaci´on. La elecci´on 2 CAP´ ITULO 1. INTRODUCCI´ ON Y OBJETIVOS 1.3. ESTIMACI´ ON DEL ERROR de la malla surge por tanto del compromiso entre la precisi´on de la soluci´on requerida y el coste computacional que se est´a dispuesto a asumir. La estrategia de mallado m´as interesante es aquella en la que se concentran mayor n´umero de elementos en la zona donde ´estos son necesarios. As´ı, se impone mayor resoluci´on en las zonas donde los gradientes son intensos. Cabe destacar que obtener el mallado adecuado ser´ıa trivial si tuvi´esemos ya una soluci´on donde se observen todos los fen´omenos relevantes. Las t´ecnicas de estimaci´on de error pueden desempe˜nar un papel crucial en esta situaci´on. Para esto es necesario que el estimador proporcione el error en cada elemento. Si se conoce el error de la simulaci´on en cada punto con fiabilidad, y se relaciona dicho error con el tama˜no de los elementos de la malla, se podr´ıa elegir un umbral de error que se considere aceptable y elaborar una malla que obtenga ese resultado. Actuando as´ı se obtendr´ıan simulaciones con un error controlado aumentando la densidad de la malla s´olo donde sea necesario. A esta estrategia se le llama mallado adaptativo. Los principales estimadores de error empleados hoy en d´ıa pueden ser clasificados en los siguientes grupos: M´etodos expl´ıcitos: Se llaman m´etodos expl´ıcitos porque emplean la soluci´on obtenida mediante elementos finitos de forma directa, sin resolver ecuaciones diferenciales. Estas t´ecnicas estiman el error a partir del residuo de la soluci´on aproximada de elementos finitos. Para ello, emplean el residuo en el interior del elemento, en las fronteras entre elementos y en el contorno. M´etodos impl´ıcitos: Estos m´etodos emplean la soluci´on aproximada de elementos finitos de forma indirecta. Para la estimaci´on del error es necesario resolver un problema de ecuaciones diferenciales aplicados a elementos individuales o a un conjunto de ellos. M´etodos de recuperaci´on: Estos m´etodos aprovechan las propiedades superconvergentes de las soluciones. Se centran en el posproceso de la soluci´on num´erica aplicado a distintos grupos de elementos. Es sabido que conforme una malla se hace m´as fina, la soluci´on num´erica converge a la soluci´on exacta. Es decir, si la malla fuera infinitamente densa la soluci´on num´erica coincidir´ıa con la soluci´on anal´ıtica. Aprovechando este concepto, el m´etodo variacional de las multiescalas, consiste en realizar una separaci´on de escalas, separando escalas resueltas y no resueltas. Este m´etodo es un caso particular de los m´etodos expl´ıcitos. La escala resuelta es la soluci´on obtenida de la malla empleada, y la escala no resuelta proporciona el error producido por la discretizaci´on. Esta descomposici´on en escala resuelta y no resuelta se introduce en la formulaci´on variacional del problema a estudiar para realizar una estimaci´on del error. 3 CAP´ ITULO 1. INTRODUCCI´ ON Y OBJETIVOS 1.4. METODOLOG´ IA 1.4. Metodolog´ıa El planteamiento general que se ha seguido a la hora de realizar este proyecto ha sido en primer lugar, establecer unos problemas concretos de tensi´on plana y hallar anal´ıticamente su soluci´on exacta. Seguidamente, se ha empleado el VMS para obtener una formulaci´on de los estimadores de error a posteriori y pronosticar el error. Por ´ultimo, se analizan las eficiencias obtenidas con los estimadores de error para los distintos problemas tratados. Los pasos seguidos en el desarrollo del trabajo se esquematizan en la figura 1.1. Figura 1.1: Esquema de trabajo Para la realizaci´on del proyecto se han empleado diversos software. La programaci´on de la soluci´on exacta, del estimador de error y del c´alculo de eficiencias se ha implementado en Fortran. Posteriormente, con la herramienta Tecplot, se han realizado las gr´aficas de resultados. Mediante la presente memoria se expone de manera descriptiva el trabajo desarrollado a lo largo de este proyecto. Primeramente, se presentan las ecuaciones que gobiernan la tensi´on plana, para asentar los conceptos b´asicos. Despu´es, se plantea la formulaci´on fuerte y la formulaci´on d´ebil del problema de elasticidad. Seguidamente, se muestra la formulaci´on del VMS y los estimadores de error empleados. Por ´ultimo, se analizan los resultados obtenidos con los estimadores de error y se obtienen unas conclusiones. 4 Cap´ıtulo 2 Elasticidad bidimensional En un cuerpo s´olido sometido a un estado el´astico plano, el rasgo m´as relevante es que tanto la matriz de tensiones como la matriz de deformaci´on permanecen constantes en una direcci´on perpendicular a un plano, llamado plano director. De esta manera, es posible tratar el problema el´astico en dos dimensiones, estudiando el problema en el plano director. Seg´un el estado tensional y de deformaci´on que se presente, se distinguen dos tipos de estados tensionales: Deformaci´on plana y tensi´on plana. En deformaci´on plana, la deformaci´on en una recta perpendicular al plano director es nula. Sin embargo, para el caso de tensi´on plana, la tensi´on es nula en la direcci´on perpendicular al plano director, existiendo deformaci´on en la misma. En este trabajo, se han estudiado problemas de tensi´on plana, si bien la formulaci´on de ambos problemas es similar. En ambos casos, el campo de desplazamientos en el plano director queda perfectamente determinado si se conocen los desplazamientos en xey. u(x, y) v(x, y) donde u(x, y) y v(x, y) son los desplazamientos en xeyrespectivamente. En el siguiente apartado, se muestra m´as concretamente, el estado de tensiones y deformaciones que se producen en tensi´on plana. 2.1. Estado tensional plano En la figura 2.1 se muestran los ejes adoptados en el plano director para un s´olido cualesquiera sometido a un estado tensional plano. Se cumple que las tensiones solamente dependen de xey, de manera que τxz = τyz=σz=0. Se denomina u,vywa los desplazamientos en los ejes x,yyzrespectivamente. 5 CAP´ ITULO 3. PROBLEMA EL´ ASTICO LINEAL 3.2. FORMULACI´ ON D´ EBIL DEL PROBLEMA los desplazamientos, u, se obtiene por medio de las ecuaciones de comportamiento y de compatibilidad. Ecuaciones de comportamiento: σx= (λ∗+ 2G)εx+λ∗εy σx= (λ∗+ 2G)εy+λ∗εx τxy =Gγxy De manera que, las componentes del tensor de Cauchy pueden expresarse en funci´on de las componentes de la matriz de deformaci´on, εx,εyyγxy. Relaci´on desplazamiento-deformaciones: εx=∂u ∂x εy=∂v ∂y γxy =∂u ∂y +∂v ∂x Por tanto, queda determinada la relaci´on entre las tensiones y los desplazamientos. Denotando como Da la matriz que relaciona tensi´on y deformaci´on, y H, al operador diferencial que relaciona las deformaciones con los desplazamientos, se pueden expresar las tensiones como σ=DHu. 3.2. Formulaci´on d´ebil del problema La formulaci´on d´ebil se obtiene a partir de la formulaci´on fuerte, multiplicando por una funci´on de peso e integrando por partes. ZΩ w· ∇ · σdΩ = ZΩ −w·fdΩ (3.1) ZΩ ∇ · (w·σ)dΩ−ZΩ ∇w:σdΩ = ZΩ −w·fdΩ (3.2) Aplicando el teorema de la divergencia y ordenando t´erminos ZΩ σ:∇wdΩ = ZΩ w·fdΩ + ZΓh w·σndΓh(3.3) 12 CAP´ ITULO 3. PROBLEMA EL´ ASTICO LINEAL 3.2. FORMULACI´ ON D´ EBIL DEL PROBLEMA ZΩ ∇w:σdΩ = ZΩ w·fdΩ + ZΓh w·hdΓh(3.4) Expresando σen funci´on de los desplazamientos, u. ZΩ ∇w: (DHu)dΩ = ZΩ w·fdΩ + ZΓh w·hdΓh(3.5) La ecuaci´on (3.5) se puede expresar de forma compacta mediante los siguientes operadores: a(·,·) es una forma bilineal (·,·) es el producto escalar en L2(Ω) (·,·)wes el producto escalar en L2(w) De este modo la ecuaci´on (3.5) queda: a(w,u) = (w,f) + (w,h)Γh(3.6) Para asegurar que las integrales que han surgido en la formulaci´on d´ebil se puedan calcular, hay que imponer como requisito a las funciones de peso, w, y a las funciones soluci´on, u, que su derivada primera al cuadrado est´e acotada. Las funciones que satisfacen esta condici´on son llamadas funciones H1. Por tanto, los espacios de funciones donde se buscan la funci´on soluci´on y las funciones de peso son los siguientes: δ={u|u∈H1,u=gen Γg} υ={w|w∈H1,w= 0 en Γg} La ´unica diferencia entre los dos espacios de funciones es el valor que deben tener en la parte del contorno donde se aplica la condici´on de contorno Dirichlet. De esta manera, las funciones de δson las que satisfacen la condici´on de contorno esencial y las funciones de υsatisfacen las condiciones de contorno esenciales homog´eneas. 13 CAP´ ITULO 3. PROBLEMA EL´ ASTICO LINEAL 3.3. M´ ETODO DE APROXIMACI´ ON DE GALERKIN 3.3. M´etodo de aproximaci´on de Galerkin Para desarrollar el m´etodo de elementos finitos, es necesario crear un espacio de funciones de dimensi´on finita que sea una aproximaci´on de δyυ. Estos espacios de funciones se denotan como δhyυhy su dimensi´on est´a relacionada con la discretizaci´on que se haga del dominio. δh⊂δ υh⊂υ Por lo tanto, la formulaci´on d´ebil discreta del problema para un n´umero de elementos nel se formula: nel X e=1 ZΩe ∇we: (DHue)dΩe= nel X e=1 ZΩe we·fdΩe+ nel X e=1 ZΓe h we·hdΓe(3.7) Considerando un elemento aislado de la discretizaci´on, como se observa en la figura 3.1, se puede expresar las componentes del vector desplazamiento en funci´on de los desplazamientos de sus nodos. u=u1N1+u2N2+u3N3+u4N4 v=v1N1+v2N2+v3N3+v4N4 Figura 3.1: Representaci´on de elemento cuadril´atero y sus variables nodales En la expresi´on anterior uyvson los desplazamientos en direcci´on xey. Por otro lado Nison las funciones de forma del nodo ial que est´a asociado. Las funciones Ni, son funciones polin´omicas que valen uno en el nodo iy cero en el resto. A su vez, las funciones de peso we, son las mismas que las funciones de forma empleadas para aproximar los desplazamientos. Las siguientes expresiones muestran las funciones de 14 CAP´ ITULO 3. PROBLEMA EL´ ASTICO LINEAL 3.4. COMPARATIVA ECUACI´ ON DE CALOR Y ELASTICIDAD forma para cuadril´ateros. we= (N1, N2, N3, N4)TPara cuadrilateros Las inc´ognitas del problema son los desplazamientos nodales ui, que se determinan resolviendo el sistema de ecuaciones lineales que resulta de la discretizaci´on. 3.4. Comparativa ecuaci´on de calor y elasticidad En esta secci´on se compara la formulaci´on de la ecuaci´on de calor y la formulaci´on de elasticidad vista en las anteriores secciones. La ecuaci´on que gobierna la ecuaci´on de calor es la ecuaci´on de Laplace. A continuaci´on se muestra las ecuaciones para un dominio en el que la conductividad en xes igual a la conductividad en y. La condici´on de contorno natural es el flujo en el contorno y la condici´on de contorno esencial es la temperatura en el contorno.    ∇2T= 0 en Ω T=gen Γg k·∂T ∂n n=hen Γh donde Tes la temperatura, gla temperatura impuesta en el contorno, hel flujo de calor impuestas en el contorno qes el vector flujo de calor q= (qx, qy)Tynel vector normal hacia fuera del contorno, n=(nx, ny). Por ´ultimo el flujo de calor se expresa como: q=−k∇Tdonde k=kx0 0ky Las principales diferencias de la formulaci´on del problema de calor y la ecuaci´on de elasticidad plana son: En la ecuaci´on de calor, la inc´ognita es un escalar, la temperatura; mientras que en elasticidad es un vector, los desplazamientos en xy en y. La matriz constitutiva de las ecuaciones de elasticidad es m´as compleja que en el caso de calor. Aunque ambas matrices son sim´etricas, la matriz constitutiva de la ecuaci´on del calor es diagonal, mientras que para elasticidad, presenta t´erminos no nulos fuera de la diagonal. En el cap´ıtulo 5 se muestra la resoluci´on de un problema concreto de la ecuaci´on de calor empleando el mismo m´etodo de estimaci´on de error que en el caso de elasticidad plana. 15 Cap´ıtulo 4 Estimaci´on de error a posteriori 4.1. M´etodo variacional de las multiescalas Sean xeydos puntos cualesquiera sobre el dominio Ω′. El problema de elasticidad es una ecuaci´on diferencial de segundo orden con las condiciones de contorno naturales y esenciales. Se puede representar de forma gen´erica    Lu=fen Ω u=gen Γg Bu=hen Γh donde L, es el operador de segundo orden y B, el operador diferencial aplicado en el contorno para la condici´on de contorno natural. Como se ha visto anteriormente, particularizando los operadores diferenciales y las variables para el caso de elasticidad lineal.    ∇ · σ+f= 0 en Ω u=gen Γg σn =hen Γh El campo de desplazamientos soluci´on del problema el´astico se puede descomponer en suma de los desplazamientos de la escala resuelta, u, y los desplazamientos de la escala no resuelta, u′. El desplazamiento u, representa la soluci´on obtenida mediante FEM y u′es la diferencia entre la soluci´on exacta y la soluci´on FEM, es decir, el error. De igual manera. las funciones de peso se pueden dividir an´alogamente en wyw′. u=u+u′ w=w+w′ Incorporando esta descomposici´on de las funciones de peso y los desplazamientos en la formulaci´on d´ebil. 16 CAP´ ITULO 4. ESTIMACI´ ON DE ERROR A POSTERIORI4.1. M´ ETODO VARIACIONAL DE LAS MULTIESCALAS ZΩ ∇ · (w+w′) : DH ·(u+u′)dΩ = ZΩ (w+w′)·fdΩ + ZΓh (w+w′)·hdΓ (4.1) Gracias a la linealidad de la ecuaci´on anterior, respeto a los desplazamientos, ´esta ecuaci´on puede ser dividida en dos problemas, uno de escala grosera y otro de escala fina. Problema de escala grosera: ZΩ ∇ · w: (DH ·(u+u′))dΩ = ZΩ w·fdΩ + ZΓh w·hdΓ (4.2) Problema de escala fina: ZΩ ∇ · w′: (DH ·(u+u′))dΩ = ZΩ w′·fdΩ + ZΓh w′·hdΓh(4.3) Tanto la formulaci´on de escala fina como la formulaci´on de escala grosera, gracias a la bilinealidad y recordando la f´ormula (3.6) se pueden agrupar en los siguientes t´erminos: a(w,u) = −a(w,u′) + (w,f) + (w,h)Γh(4.4) a(w′,u′) = −a(w′,u) + (w′,f) + (w′,h)Γh(4.5) En el m´etodo de elementos finitos, el dominio del problema, Ω, se subdivide en elementos con dominio Ωey contorno Γe(e=1,2,. . . , nel). La agrupaci´on de los dominios y el contorno de los elementos se denotan respectivamente: Ω′=∪nel e=1Ωe Γ′=∪nel e=1Γe Por ´ultimo, el conjunto del contorno interior de los elementos con los que se ha discretizado el dominio se expresa como Γint = Γ′\Γ. Se denomina J·Kal operador salto de una funci´on a trav´es de una discontinuidad, como sucede en los contorno entre elementos. Seg´un la notaci´on de la figura 4.1, referida al 17 CAP´ ITULO 4. ESTIMACI´ ON DE ERROR A POSTERIORI4.1. M´ ETODO VARIACIONAL DE LAS MULTIESCALAS vector normal de los elementos que comparten el mismo contorno, el salto de vse expresa como: Jv·nK=v+·n++v− ·n−(4.6) Figura 4.1: Notaci´on para el operador salto Las funciones para elementos finitos son suaves dentro del elemento pero bastas en el contorno entre elementos. Por tanto, la integraci´on de a(w′,u) resulta a(w′,u) = Pnel i=1 a(w′,u)Ωe =Pnel i=1[(w′,Lu)Ωe+ (w′,Bu)Γe] = (w′,Lu)Ω′+ (w′,JBuK)Γint + (w′,Bu)Γh De esta manera el problema de la escala fina queda a(w′,u′) = −(w′,Lu)Ω′−(w′,JBuK)Γint −(w′,Bu)Γh(4.7) +(w′,f) + (w′, h)Γh =−(w′,Lu−f)Ω′−(w′,JBuK)Γint −(w′,Bu−h)Γh Agrupando en un mismo t´ermino, los saltos en el contorno interior entre los elementos JBuKy en el contorno exterior Bu−h. a(w′,u′) = −(w′,Lu−f)Ω′−(w′,JBuK)Γe(4.8) donde el t´ermino JBuKrepresenta las saltos que se producen entre las tensiones en las fronteras entre elementos. 18 CAP´ ITULO 4. ESTIMACI´ ON DE ERROR A POSTERIORI 4.2. ESTIMACI´ ON DE ERROR 4.2. Estimaci´on de error La escala fina se puede resolver de forma anal´ıtica empleando las funciones de Green. Llamando G′(x,y) a la matriz de funciones de Green del problema de escalas finas resulta u′(x)≈ − ZΩ′ G′(x,y)(Lu−f)(y)dΩy−ZΓe y G′(x,y)(JBuK)(y)dΓe y(4.9) donde G′(x,y) = ge x(x,y)ge y(x,y) ge y(x,y)ge x(x,y) Las funciones ge x(x,y) y ge y(x,y) son las funciones de green definidas en el eje xy en el eje yrespectivamente. En el presente problema de elasticidad lineal, el primer t´ermino de la ecuaci´on (4.9) es despreciable para elementos lineales y bilineales y, por tanto, la fuente de error del m´etodo de elementos finitos proviene de los saltos entre las acciones exteriores y la tensi´on que presentan los elementos situados en el contorno y los saltos de tensiones que se producen en los contornos de los elementos. u′(x)≈ − ZΓe y G′(x,y)(JBuK)(y)dΓe y(4.11) Aplicando la inecuaci´on de H¨olders u′(x)≤G′(x,y)Lp(Γe y)||JBuK||Lq(Γe)(4.12) donde 1 ≤p,q≤ ∞ , 1/p + 1/q = 1. Tomando la norma Lr u′(x)Lr(Ωe)≤kG′(x,y)kLp(Γe y)Lr(Ωe x)||JBuK||Lq(Γe)(4.13) Seleccionando p= 1 y q=∞ u′(x)Lr(Ωe)≤kG′(x,y)kL1(Γe y)Lr(Ωe x)||JBuK||L∞(Γe)(4.14) 19 CAP´ ITULO 4. ESTIMACI´ ON DE ERROR A POSTERIORI 4.3. FUNCI´ ON BUBBLE Mediante el teorema de la traza, es posible realizar una aproximaci´on de la norma de la funci´on de Green de la escala fina sobre el contorno Γey de la norma de la funci´on de Green en el interior del elemento Ωemediante la siguiente expresi´on: kG′(x,y)kL1(Γe y)Lr(Ωe x)≈1 2 meas(Γe) meas(Ωe)kGe(x,y)kL1(Ωe y)Lr(Ωe x)(4.15) Si la funci´on Ge(x,y) no cambia de signo en el interior del elemento kGe(x,y)kL1(Ωe y)Lr(Ωe x)=kBe(x)kLr(Ωe x)(4.16) donde la matriz Be(x) es la matriz de funciones bubble del residuo. La manera de determinar ´estas funciones bubble se estudia en el siguiente apartado. Be=be 0xbe 0y be 0ybe 0x As´ı pues, la expresi´on (4.14) se puede expresar como u′(x)Lr(Ωe)≤ kBe(x)kLr(Ωe x)||JBuK||L∞(Γe)(4.17) En este caso, ||JBuK||L∞(Γe)representa el salto m´aximo que existe entre las tensiones en el contorno del elemento. La funci´on bubble aporta al estimador de error la informaci´on espec´ıfica del problema f´ısico a resolver. En la formulaci´on del estimador de error, la funci´on bubble se introduce a trav´es de la matriz τerr el,Lr como se ver´a m´as adelante. 4.3. Funci´on bubble Para encontrar la soluci´on de desplazamientos de la escala fina, u′, es necesario determinar la funci´on bubble, que a su vez es la integral de la funci´on Green. La formulaci´on para obtener la funci´on bubble o burbuja es: LBe=Ien Ωe Be= 0 en Γe donde Bees la funci´on bubble y Ies la matriz identidad. 20 CAP´ ITULO 4. ESTIMACI´ ON DE ERROR A POSTERIORI 4.3. FUNCI´ ON BUBBLE El operador diferencial, L, es el mismo que hemos empleado anteriormente para resolver el problema de tensi´on plana. El problema se eval´ua en un elemento Ωe, cuyo contorno es Γe. La funci´on bubble se ha obtenido de forma num´erica. Para ello, mediante el problema de elementos finitos empleado para resolver los problemas de tensi´on plana, se ha resuel- to un problema con carga distribuida unidad sobre un ´area rectangular y condiciones de contorno Dirichlet homog´eneas. El m´odulo el´astico y el coeficiente de Poisson que se han impuesto son los mismos que los del problema a analizar. En la figura 4.2 se muestra el planteamiento del problema. Figura 4.2: Planteamiento del problema para hallar la funci´on bubble La soluci´on en desplazamientos del problema anterior, determina la funci´on bubble. Para realizar la integral de la funci´on bubble se dispuso con una malla muy fina para que la integral fuera lo m´as exacta posible. En la figura 4.3 se muestra tanto la malla como los vectores desplazamiento de cada nodo. Los desplazamientos que se obtienen tanto en el eje xcomo en el eje y, se integran por separado para obtener la matriz τerr el que nos servir´a posteriormente para obtener la estimaci´on del error. Como se aprecia, existen desplazamientos tanto en xcomo en y, que se denominar´an ub= (ubx, uby). Por tanto, si se toma la norma L1yL2de la funci´on bubble con dichos desplazamientos Norma L1:||be 0x||L1(Ω) =Pnel i=1 RΩe|ubx|dΩebe 0yL1(Ω) =Pnel i=1 RΩe|uby|dΩe(4.18) 21 CAP´ ITULO 5. RESULTADOS 5.2. PROBLEMAS SOMETIDOS A FLEXI ´ ON EN EL CONTORNO Problema (c) u=Mx EIx x·y−νMy 2IyEx2−My 2EIy y2+My 2IyEc2 v=−νMx 2EIx y2−Mxx2 EIx +νMx 2EIx c2+My EIy xy donde Mx,MyyNson los esfuerzos de momento en x, momento en yy axil respectivamente. Ix=Rc −cy2·dy es el momento de inercia en el eje xyIy=Rl/2 −l/2x2·dx es el momento de inercia en el eje y. 5.2.2. Soluci´on de elementos finitos Para realizar las simulaciones se ha realizado un mallado de 20 elementos en xy 4 elementos en ycomo se muestra en la figura 5.5. Figura 5.5: Forma del mallado En las figuras 5.6, 5.7 y 5.8 se muestra la deformada de los tres casos distintos analizados. Se muestra las cargas aplicadas y los vectores del desplazamiento de los nodos. 28 CAP´ ITULO 5. RESULTADOS 5.2. PROBLEMAS SOMETIDOS A FLEXI ´ ON EN EL CONTORNO Figura 5.6: Deformada del problema (a) Figura 5.7: Deformada del problema (b) Figura 5.8: Deformada del problema (c) 29 CAP´ ITULO 5. RESULTADOS 5.2. PROBLEMAS SOMETIDOS A FLEXI ´ ON EN EL CONTORNO 5.2.3. An´alisis de la eficiencia del error Dela figura 5.9 a 5.14, se muestra en mapa de colores la eficiencia local tanto en el eje xcomo en el eje ypara los problemas (a), (b) y (c). Se muestra la eficiencia tanto con la norma L1como con la norma L2. Problema (a) En este caso, se ha aplicado un momento en los extremos del dominio de Mx= 1000 N m. En las figuras 5.9 y 5.10 se observan las eficiencias locales. Figura 5.9: Eficiencia en x. Problema (a) Figura 5.10: Eficiencia en y. Problema (a) 30 CAP´ ITULO 5. RESULTADOS 5.2. PROBLEMAS SOMETIDOS A FLEXI ´ ON EN EL CONTORNO Problema (b) El esfuerzo axil aplicado es Nx=−1000N y el el momento en x es Mx= 1000N m Figura 5.11: Eficiencia en x. Problema (b) Figura 5.12: Eficiencia en y. Problema (b) Se observa que la eficiencia es la misma para el problema (a) que para el problema (b). Este resultado es coherente, ya que la ´unica diferencia entre ambos problemas es el esfuerzo axil en el eje xque existe en el problema (b). La soluci´on que proporciona el m´etodo de elementos finitos de una barra traccionada es exacta, ya que en este caso los desplazamientos son lineales. Por esta raz´on, el esfuerzo axil aplicado no introduce error alguno en la soluci´on de elementos finitos del problema (b). Por otro 31 CAP´ ITULO 5. RESULTADOS 5.2. PROBLEMAS SOMETIDOS A FLEXI ´ ON EN EL CONTORNO lado, nuestro estimador de error tampoco se ve afectado. De ah´ı, que la eficiencia entre los dos problemas sea la misma. Problema (c) En este problema, como ya se ha comentado, se aplica momento de flexi´on en ambas caras del dominio rectangular. En la cara superior e inferior el momento es My= 100000N m y en las caras laterales el momento es Mx= 1000N m. La diferencia de la magnitud del momento aplicado en el eje xy en eje yse debe a que la inercia en y es mayor que en xy por tanto se ha aumentado el momento Mypara que la tensi´on en las caras del dominio sean del mismo orden. Figura 5.13: Eficiencia en x. Problema (c) Figura 5.14: Eficiencia en y. Problema (c) 32 CAP´ ITULO 5. RESULTADOS 5.2. PROBLEMAS SOMETIDOS A FLEXI ´ ON EN EL CONTORNO Gr´aficos de eficiencias locales Para evaluar la eficiencia local en los elementos de los problemas, se han agrupado dichas eficiencias seg´un la tabla 4.1. En las figuras 5.15, 5.16 y 5.17 se representan los diagramas de barras que muestran la frecuencia con que se repite una eficiencia local en los elementos tanto en la norma L1como L2. Figura 5.15: Eficiencia en xy en y. Problema (a) Figura 5.16: Eficiencia en xy en y. Problema (b) 33 CAP´ ITULO 5. RESULTADOS 5.2. PROBLEMAS SOMETIDOS A FLEXI ´ ON EN EL CONTORNO Figura 5.17: Eficiencia en xy en y. Problema (c) Eficiencia Global Por ´ultimo, en la siguiente tabla 5.1 se muestran las eficiencias globales obtenidas tanto en el eje xcomo en el y. Eficiencia x Eficiencia y Norma L1 Norma L2 Norma L1 Norma L2 Probl (a) 2.475 2.73 0.308 0.358 Probl (b) 2.475 2.73 0.308 0.358 Probl (c) 1.544 1.816 0.309 0.360 Tabla 5.1: Eficiencia globales en xeyen problemas a, b y c 5.2.4. An´alisis de sensibilidad del estimador Es interesante comprobar c´omo se comporta el estimador de error al variar los principales par´ametros que definen el problema. En este caso, se ha estudiado la eficiencia para diferentes m´odulos el´asticos y coeficientes de Poisson, los cuales definen el comportamien- to de un material el´astico lineal. Tambi´en se ha analizado la influencia del tama˜no del elemento de la malla. Modificaci´on de M´odulo el´astico, E Para el an´alisis de sensibilidad aplicado al m´odulo el´astico, E, se ha tomado ν= 0.3 yE= 2.1e10,2.1e11y2.1e12 N m. Se ha observado que el estimador proporciona la misma eficiencia independientemente su valor, como se aprecia en la figura 5.18. Esto se debe a la proporcionalidad que existe entre el m´odulo el´astico y los desplazamientos, que es la misma que existe entre el m´odulo el´astico y el error estimado. 34 CAP´ ITULO 5. RESULTADOS 5.2. PROBLEMAS SOMETIDOS A FLEXI ´ ON EN EL CONTORNO Figura 5.18: Influencia del m´odulo el´astico en la eficiencia. Norma L1yL2 De esta manera, al doblar el valor del m´odulo el´astico, los desplazamientos que proporciona la soluci´on de elementos finitos y la soluci´on exacta es justamente el doble. Por tanto, el error real que comete el m´etodo de elementos finitos es el doble. Por otro lado, el error estimado tambi´en se dobla ya que la funci´on bubble dobla su valor y el estimador es proporcional a esta funci´on bubble. Por tanto, al realizar el c´alculo de la eficiencia, en la que se divide error estimado y error real, se obtiene el mismo valor. Coeficiente de Poisson, ν Desde un punto de vista f´ısico, los valores que puede tomar el coeficiente de Poisson de un material va de 0 a 0,5 sin llegar a estos extremos. Por ello, se ha tomado como m´odulo el´astico E= 2.1e11 N m y se ha realizado un barrido del par´ametro coeficiente de Poisson desde 0,05 a 0,45. En las figura 5.19 y 5.20 se muestra los valores de eficiencia global obtenidos para los tres problemas analizados y para las normas L1yL2. Los problemas (a) y (b) se han agrupado al tener la misma eficiencia. Tama˜no del elemento Una cuesti´on importante a analizar en un estimador de error es ver si la eficiencia se ve afectada por el tama˜no de los elementos. Para estudiar la sensibilidad de la eficiencia, se ha partido de un tama˜no de elemento hy se ha ido haciendo el elemento m´as peque˜no h/2 yh/4. En este caso, hes el lado del elemento cuadril´atero. 35 CAP´ ITULO 5. RESULTADOS 5.2. PROBLEMAS SOMETIDOS A FLEXI ´ ON EN EL CONTORNO Figura 5.19: Influencia del coeficiente de Poisson en la eficiencia. Norma L1 Figura 5.20: Influencia del coeficiente de Poisson en la eficiencia. Norma L2 En las figuras 5.21 y 5.22 aparece la eficiencia global obtenida con los distintos tama˜nos de malla para las normas L1yL2con ν= 0.3 y E= 2.1e11 N m. Se aprecia que la influencia es m´ınima, exceptuando la eficiencia para el eje xen los problemas (a) y (b). Figura 5.21: Influencia del tama˜no de malla en la eficiencia. Problema a y b 36 CAP´ ITULO 5. RESULTADOS 5.3. VIGA EN VOLADIZO Y VIGA BIAPOYADA Figura 5.22: Influencia del tama˜no de malla en la eficiencia. Problema c 5.3. Viga en voladizo y viga biapoyada En este apartado, se muestran los resultados obtenidos en las simulaciones realizadas para dos casos en los que el estado tensional de los cuerpos es de mayor complejidad. En la figura 5.23 aparecen los problemas que se han estudiado. El primer problema, se trata de una viga rectangular en voladizo con una carga P en el extremo opuesto del empotramiento. El segundo caso, el problema a resolver es una viga rectangular biapoyada en la que se ha aplicado una carga distribuida en la cara superior. Figura 5.23: Problemas de viga empotrada y viga apoyada Al igual que en la secci´on anterior, se ha seleccionado un dominio con una longitud de 10 m y una anchura de 2 m. 5.3.1. Soluci´on anal´ıtica de los problemas Como en los problemas anteriormente analizados, hay que conocer la soluci´on exacta del problema a tratar. De esta manera, se puede hallar el error que se comete con el m´etodo de elementos finitos. A continuaci´on, se muestran las soluciones anal´ıticas de los 37 Bibliograf´ıa LIBROS [1] Hauke, G.;Introducci´on al m´etodo de elementos finitos en fluidodin´amica computacional, Primera edici´on, 2010. [2] Timoshenko, S.P., Goodier J.N. Theory of elasticity, McGraw-Hill Book Company, 1970. [3] Ortiz Berrocal, L., Elasticidad, 3aedici´on, McGrawHill, 1998. [4] O˜nate, E. C´alculo de estructuras por el m´etodo de elementos finitos. Analisis el´astico lineal, 2aedici´on, CIMNE, 1995. [5] Hughes, T.J.R. The finite element method: Linear static and dynamic finite element analysis, Prentice-Hall, Engelwoods Cliffs,1987. [6] Ainsworth, M., Oden J.T. A posteriori error estimation in finite element analysis, Wiley-interscience, 2000. [7] Smith, I.M., Griffith, D.V. Programming the finite element method, 4th edition, Wiley, 2004. [8] Roach, G.F. Green’s function, 2nd edition, Cambridge university press, 1982. [9] Stakgold, I., Holst M., Green´s function and boundary value problems, 3rd edition, Wiley, 2011. ART´ ICULOS [10] Hauke G., Doweidar M.H. and Fuster D. A posteriori error estimation for computational fluid dynamics: The variational multiscale approach. (Eds.: E. Ramm, R. Borst), Series: Lecture Notes in Applied and Computational Mechanics, Springer, Vol 55 pp. 1–20, 2010. ISSN 1613-7736, ISBN 978-90-481-9808-5. [11] Hauke G., Doweidar M.H. and Miana M. Proper intrisic scales for a-posteriori multiscale error estimation. Comput. Methods Appl. Mech. Engrg., 2005. [12] Hauke G., Doweidar M.H. and Miana M. The multiscale approach to eror estimation and adaptivity. Comput. Methods Appl. Mech. Engrg., 2005. 44 BIBLIOGRAF´ IA BIBLIOGRAF´ IA [13] Masud A., Truster T.J. and Bergman M. A variational multiscale a posteriori error estimation method for mixed form of nearly incompressible elasticity. Comput. Methods Appl. Mech. Engrg., 2011. 45 Anexos 46 Anexo A Soluci´on anal´ıtica: Problema flexi´on Pura En este caso, el problema es una viga de gran canto sometida a flexi´on pura. Las dimensiones y ejes adoptados para el problema de flexi´on pura se observan en la figura A.1. Seg´un se muestra, las cargas que se aplican en el contorno del dominio, debido al momento aplicado en los extremos, es una carga distribuida de forma triangular, es decir, la magnitud de la carga aplicada es proporcional a la distancia respecto a la fibra neutra. Figura A.1: Problema flexi´on pura. Dimensiones y ejes La distribuci´on de la carga en el contorno es igual a las tensiones que se producen en la viga en una secci´on perpendicular al eje x. Por otro lado, la tensi´on en y,σy, y el cortante, τxy, son nulos. 47 ANEXO A. SOLUCI ´ ON ANAL´ ITICA: PROBLEMA FLEXI´ ON PURA σx=Mx Ix y σy= 0 τxy = 0 (A.1) donde Mxes el momento aplicado e Ix=Rc −cy2·dy es el momento de inercia. Conociendo el estado tensional en el dominio del problema, es inmediato hallar el estado de deformaci´on. εx=∂u ∂x =1 E(σx−νσy) = Mx EIx y(A.2) εy=∂v ∂y =1 E(σy−νσx) = −νMx EIx y(A.3) Integrando e imponiendo las condiciones de contorno apropiadas se obtiene el campo de desplazamientos. u=Mx EIx y·x+f1(y) (A.4) v=−νMx 2EIx y2+f2(x) (A.5) Para determinar f1(x) y f2(x) de las ecuaciones anteriores, se aplica la condici´on de que el esfuerzo cortante es cero. τxy =G(∂u ∂y +∂v ∂x) = 0 →Mx EIx x+∂f1(y) ∂y +∂f2(x) ∂x = 0 (A.6) Por lo tanto, para cumplir la ecuaci´on anterior las funciones f1(y) y f2(x) se definen como se muestra a continuaci´on. f2(x) = −Mx 2EIxx2+A f1(y) = 0 Sustituyendo en (A.4) y (A.5), las funciones f2(x) y f1(y) u=Mx EIx y·x(A.7) v=−νMx 2EIx y2−Mx 2EIx x2+A(A.8) 48 ANEXO A. SOLUCI ´ ON ANAL´ ITICA: PROBLEMA FLEXI´ ON PURA La constante Ase debe obtener aplicando una condici´on de contorno. En esta caso, la condici´on de contorno es que el desplazamiento vertical en el apoyo es cero. x= 0 y=−c} → v=−ν Mx 2EIxc2+A= 0 ⇒A=νMx 2EIxc2 Finalmente, los desplazamientos para el problema de flexi´on pura en una viga de gran canto quedan: u=Mx EIx y·x(A.9) v=−νMx 2EIx y2−Mxx2 2EIx +νMx 2EIx c2(A.10) 49 Anexo B Soluci´on anal´ıtica: Problema flexi´on compuesta Este problema presenta la misma geometr´ıa y condiciones de apoyo que en el caso del anexo A. La ´unica diferencia es la carga a la que est´a sometida la viga. En los laterales, se ha a˜nadido al momento del problema de flexi´on pura, un esfuerzo axil. En la figura B.1 se muestra la disposici´on de las cargas y las tensiones aplicadas en el contorno. Figura B.1: Problema flexi´on compuesta. Dimensiones y ejes Debido a la elasticidad lineal, se puede obtener la soluci´on anal´ıtica mediante el principio de superposici´on, a partir de la soluci´on del anexo A. De esta manera, la soluci´on exacta es la suma de la soluci´on del anexo A y la soluci´on de la viga sometida a esfuerzo axil. 50 ANEXO B. SOLUCI ´ ON ANAL´ ITICA: PROBLEMA FLEXI´ ON COMPUESTA La viga sometida a una carga axil, Nx, el estado tensional y de deformaci´on es Tensiones: σx=Nx 2c σy=τxy = 0 Deformaciones: εx=σx E εy=−νσx E Los desplazamientos, uyv, se determinan integrando las expresiones de las deformaciones. εx=∂u ∂x →u=Nx 2cEx+A(B.1) εy=−ν∂v ∂y →v=−νNx 2cE y+B(B.2) Las constantes de integraci´on se determinan imponiendo las condiciones de contorno. En este caso, la condici´on es que tanto el desplazamiento en xcomo en yes cero en el apoyo. x= 0 y=−c} → u=A= 0 v= 0 = νNx 2cE c+B→B=−νNx 2cE c Por tanto, los desplazamientos generados por el esfuerzo axil son u=Nx 2cE x(B.3) v=−νNx 2cE (y+c) (B.4) Sumando esta soluci´on a la soluci´on del anexo A se obtiene la soluci´on anal´ıtica de este problema u=Mx EIx y·x+Nx 2cE x(B.5) v=−νMx 2EIx y2−Mxx2 EIx +νMx 2EIx c2−νNx 2cE (y+c) (B.6) 51 Anexo C Soluci´on anal´ıtica: Problema de flexi´on en los dos ejes Este problema tiene una geometr´ıa y unas condiciones de contorno iguales al anexo A. La carga a la que est´a sometido es tambi´en la misma que en el anexo A pero a˜nadiendo un momento, My, en la cara superior de la viga, como se aprecia en la figura C.1. Por tanto, la soluci´on anal´ıtica de desplazamientos se puede descomponer como suma de dos soluciones parciales. La primera de ellas ser´a la soluci´on debida al momento Mxaplicado en las caras laterales, obtenida en el anexo A. La otra soluci´on, es la debida al momento Myque solicita la cara superior e inferior, y que se expone a continuaci´on. Figura C.1: Problema flexi´on en los dos ejes. Dimensiones y ejes En primer lugar, para hallar los desplazamientos que produce el momento My, es necesario conocer las tensiones y deformaciones que produce este momento. 52 ANEXO C. SOLUCI ´ ON ANAL´ ITICA: PROBLEMA DE FLEXI ´ ON EN LOS DOS EJES Tensiones: σy=My Iyx σx=τxy = 0 Deformaciones: εx=∂u ∂x =−νMy EIyx εy=∂v ∂y =My EIyx Los desplazamientos uyvse obtienen integrando las ecuaciones de deformaci´on. εx=∂u ∂x →u=−νMy 2EIy x2+f1(y) (C.1) εy=−ν∂v ∂y →v=My EIy xy +f2(x) (C.2) Las ecuaciones f1(y) y f2(x) se pueden determinar a partir de la condici´on de que el cortante es nulo. τxy = 0 →∂u ∂y +∂v ∂x = 0 (C.3) Sustituyendo uyven la ecuaci´on (C.3), Myy EIy +∂f1(y) ∂y +∂f2(x) ∂x = 0 (C.4) por tanto ∂f1(y) ∂y =−Myy EIyy∂f2(x) ∂x = 0 y las funciones f1(y) y f2(x) quedan f2(x) = A f1(y) = −Myy2 2EIy +B(C.5) Los desplazamientos uyvse pueden expresar como u=−νMy 2EIy x2−Myy2 2EIy +B(C.6) 53 ANEXO E. SOLUCI ´ ON ANAL´ ITICA: PROBLEMA VIGA APOYADA - Cara laterales: Zc −c τxydy =±ql, Zc −c σxdy = 0,Zc −c σxydy = 0 (E.2) Las dos ´ultimas ecuaciones de (E.2) establecen, que sobre los extremos no hay esfuerzo longitudinal ni ning´un momento flector. Las condiciones (E.1) y (E.2) se cumplen al combinar varias de las funciones de Airy. φ1=a 6x2y3φ2=b 2x2φ3=c 2x2y(E.3) A partir de estas funciones de Airy se generan las siguientes tensiones σx=a(x2y−2 3y3) σy=1 3ay3+cy +b τxy =−axy2−cx Reemplazando estas componentes de tension en las condiciones (E.1) se determinan las constantes a, b y c. a=−3 4 q c3b=−q 2c=3 4 q c(E.4) Teniendo en cuenta que 2c3/3 es el momento de inercia Ide la secci´on transversal. σx=−3q 4c3(x2y−2 2y3) = −q 2I(x2y−2 2y3) (E.5) σy=−3q 4c3(1 3y3−c2y+2 3c3) = −q 2I(1 3y3−c2y+2 3c3) (E.6) τxy =−3q 4c3(c2−y2)x=−q 2I(c2−y2)x(E.7) Como se puede ver, esta soluci´on de tensiones cumplen con las condiciones de tensiones (E.1) y las dos primeras condiciones de (E.2). La ´ultima condici´on no es satisfecha. ´ Esta condici´on expresa que en las caras laterales, no existe un par aplicado. Por tanto, impondremos una flexi´on simple σx=d·y, que proviene de la funci´on de Airy φ4=d 6y3. El coeficiente dse obtiene al imponer precisamente la condici´on comentada en los extremos x=±l. Zc −c σxydy =Zc −c [−3q 4c3(l2y−2 2y3) + dy]ydy = 0 (E.8) 60 ANEXO E. SOLUCI ´ ON ANAL´ ITICA: PROBLEMA VIGA APOYADA por tanto d=3 4 q c(l2 c2−2 5) (E.9) De manera que la tensi´on σxqueda σx=−q 2I(x2y−2 2y3) + 3 4 q c(l2 c2−2 5)y=−q 2I(x2y−2 2y3) + q 2I(l2−2c2 5)y(E.10) Una vez obtenido el estado tensional, se puede hallar f´acilmente las deformaciones aplicando la ley de Hooke εx=∂u ∂x =σx E=−q 2EI (x2y−2 3y3) + q 2EI (l2−2c2 5)y(E.11) εy=∂v ∂y =σy E=−q 2EI (1 3y3−c2y+2 3c3) (E.12) γxy =∂u ∂y +∂v ∂x =τxy G=−q 2GI (c2−y2)x(E.13) Por tanto los desplazamientos se pueden expresar como u=−q 2EI (x3y 3−2 3y3x) + q 2EI (l2−2c2 5)yx +f1(y) (E.14) v=−q 2EI (1 12y4−1 2c2y2+2 3c3y) + f2(x) (E.15) Las funciones f1(y) y f2(x) se determinan sustituyendo las funciones (E.14) y (E.15) en la ecuaci´on (E.13). De esta manera las funciones f1(y) y f2(x) que resultan f1(y) = A(E.16) f2(x) = −2 5 q EI c2x2+q EI x4 24 −νqc2x2 4EI +B(E.17) donde AyBson constantes de integraci´on. ´ Estas se determinan mediante las condiciones de apoyo. Sustituyendo en (E.14) y (E.15) las funciones f1(y) y f2(x) y reordenando t´erminos las funciones uyvquedan: 61 ANEXO E. SOLUCI ´ ON ANAL´ ITICA: PROBLEMA VIGA APOYADA u=q 2EI (l2x−x3 3)y+x(2 3y3−2 5c2y) + νx(1 3y3−c2y+2 3c3)+A(E.18) v=−q 2EI y4 12 −c2y2 2+2 3c3y+ν[(l2−x2)y2 2+y4 6−1 5c2y2](E.19) −q 2EI [l2x2 2−x4 12 −1 5c2x2+ (1 + 1 2ν)c2x2] + B La constante Ade (E.21) se observa que es cero (A= 0), ya que el desplazamiento horizontal en la secci´on central es nulo. Por otro lado, la constante Bse puede hallar con la condici´on de apoyo en los extremos, en la que el desplazamiento ves nulo. B=5 24 ql4 EI [1 + 12 5 c2 l2(4 5+ν 2)] (E.20) Por tanto, los desplazamientos uyvde la soluci´on anal´ıtica son: u=q 2EI (l2x−x3 3)y+x(2 3y3−2 5c2y) + νx(1 3y3−c2y+2 3c3)(E.21) v=−q 2EI y4 12 −c2y2 2+2 3c3y+ν[(l2−x2)y2 2+y4 6−1 5c2y2]−(E.22) −q 2EI [l2x2 2−x4 12 −1 5c2x2+ (1 + 1 2ν)c2x2]+ +5 24 ql4 EI [1 + 12 5 c2 l2(4 5+ν 2)] 62