Aspectos básicos de la modelización matemática y numérica de los medios continuos

FACULTAD DE MATEMÁTICAS DEPARTAMENTO DE ECUACIONES DIFERENCIALES Y ANÁLISIS NUMÉRICO Aspectos básicos de la Modelización Matemática y Numérica de los Medios Continuos Trabajo de Fin de Grado realizado por Susana Pérez Martínez Tutor Tomás Chacón Rebollo Septiembre de 2019 Índice general Índice general 1 1. Motivación 5 2. Planteamiento del problema abstracto 7 3. Discretización 15 3.1. Planteamiento del problema aproximado .............. 15 3.2. Cálculo efectivo de la solución discreta ............... 18 3.3. Interpretación geométrica del método de Galerkin ......... 18 3.4. Buen planteamiento del problema .................. 19 3.5. Análisis de Error ........................... 21 4. Ecuación de Elasticidad 26 4.1. Modelización del problema ...................... 27 4.1.1. Ecuaciones de Equilibrio ................... 27 4.1.2. Relaciones constitutivas de materiales elásticos . . . . . . 31 4.1.3. Las ecuaciones de Elasticidad no Lineal .......... 35 4.1.4. Las ecuaciones de Elasticidad Lineal ............ 35 4.2. Elasticidad Lineal .......................... 36 4.2.1. Desigualdades de Korn .................... 36 4.2.2. Existencia, Estabilidad y Regularidad de la Solución Débil 40 4.3. Discretización del problema ..................... 45 4.4. Implementación en FreeFem++ ................... 46 4.4.1. Caso elástico en 2D ..................... 47 4.4.2. Caso elástico en 3D ..................... 48 5. Ecuación de Difusión-Convección-Reacción 50 5.1. Introducción al Problema ...................... 50 5.1.1. Modelización del Problema ................. 50 5.1.2. Justificación física de transmisión del calor en sólidos . . . 54 5.2. Planteamiento del Problema ..................... 57 5.3. Discretización del Problema ..................... 61 5.3.1. La Aproximación Conforme H1(Ω) ............. 62 1 ÍNDICE GENERAL 2 5.3.2. La Aproximación no Conforme Crouzeix-Raviart ..... 65 5.4. Implementaciones en FreeFem++ .................. 69 5.4.1. Primer Caso: Comparativa Materiales ........... 69 5.4.2. Segundo Caso: Habitación .................. 75 A. Elementos del Análisis Funcional 79 A.1. Espacios de Hilbert y de Banach .................. 79 A.2. Espacios Lp.............................. 83 A.3. Teoría de Distribuciones ....................... 84 A.3.1. Sobre el concepto de solución de una EDP ......... 84 A.3.2. Espacio de Distribuciones D0(Ω) .............. 86 A.3.3. Derivación débil ....................... 90 A.4. Espacios de Sobolev ......................... 91 A.4.1. Los espacios de Hilbert Hk(Ω) ............... 93 A.5. Convergencias ............................. 95 B. Código FreeFem++ 97 B.1. Simulaciones Elasticidad ....................... 97 B.2. Simulaciones Convección-Difusión ..................101 3 Abstract The main purpose of this paper is the construction and treatment of equations that are useful in Continuum Mechanichs from both, physical and mathematical points of view. In particular, we will be considering the Advection Difussión equation and the Elastisicy equation. We will work in Sobolev Spaces, which are the natural enviroment for PDEs to be developed. Through the theorems of Lax Milgram and Banach-Neças-Babuska, we will study under which conditions a configuration setting is well posed. Once the existence and uniqueness of the analytical solution is proven, we will approximate our problem using Galerkin and Petrov Galerkin methods. The existence and uniqueness of the approximate solution is studied with very similar results to the analytical ones, and its covergence to the exact solution will be a result of having both, consistency and stability properties. Next, the Elasticity and Advection-Diffusion equations will be modeled. We will study the conformal approximation of Elasticity problems with different boundary conditions and the conformal and nonconformal approximation of Advection-Difussion problems, also with different boundary conditions. Finally, we will briefly present some implementations of the solutions of these problems in FreeFem++. 4 Capítulo 1 Motivación En su deseo de conocer la Naturaleza, el ser humano como ente no solo animal, sino también racional ha buscado desde tiempos inmemorables la creación de modelos que permitan anteponerse al comportamiento de éste, nuestro Universo. Lo que como recién nacido fueron pruebas de error-acierto o copias de comportamientos que resultaban efectivos, de la mano del proceso de maduración natural de la mentalidad matemática del hombre, se han ido traduciendo en modelos analíticos. Dichos modelos, fieles a una Filosofía aristotélica, si bien tienen origen en lo empírico, no anhelan, sino y como objetivo último la trascendencia a otro mundo mucho más imperecedero y elegante en cuanto a su forma: la lógica matemática y su perpleja capacidad de adaptación al mundo caduco. Madre de cualquier disciplina natural actual (Medicina, Biología, Física...), padre del antropocéntrico estado de bienestar, e hijo de ese germen insaciable por conocer y comprender que caracteriza a nuestra especie, el Análisis Funcional y Modelado de Ecuaciones en Derivadas Parciales (EDPs) constituye una de las más boyantes ramas de estudio del momento. En particular, este trabajo se centra en la modelización matemática de dos aspectos básicos de los que nace la Arquitectura de nuestros días: la convección difusión tratando los fluídos y la elasticidad en los sólidos. Ocurre que, generalmente, los modelos no tienen soluciones analíticas y tenemos que recurrir a herramientas propias del Análisis Numérico para poder aproximarnos a dicho comportamiento. Entre estas herramientas se encuentra una fundamental y que marca el sendero del trabajo que aqui se presenta: el Método de los Elementos Finitos, cuyo nombre fue primeramente acuñado por Ray. W. Clough en 1952 en el marco del Boeing Summer Faculty Program. Si bien estos métodos surgieron en el ámbito de la ingeniería, poco a poco fueron sirviendo para solventar una familia de problemas cada vez más amplia y variada. Como consecuencia, fue increscendo la necesidad de crear una teoría matemática que sustentara de algún modo la validez de estos métodos; produciéndose su salto 5 CAPÍTULO 1. MOTIVACIÓN 6 intrahistórico al mundo académico. Todo esto, unido a la necesidad de dar una base formal fuerte a la nueva ”función impulso o Delta de Dirac” tan necesaria en las ramas de la Física (desde Teoría de Circuitos hasta la entonces recién nacida Mecánica Cuántica, pasando por conceptos de estadística de partículas y principios de incertidumbre), Laurent Schwartz en su deseo subversivo por revolucionar el mundo de sus días, desarrolla una teoría por la que fue galardonado con la Medalla Fields: la Teoría de Distribuciones. El concepto de función se relaja. Las nociones de derivabilidad también. Las distribuciones han nacido, y con ellas la tierra prometida donde poder coexistir desde 1960 hasta nuestros días: Los espacios de Sobolev. Es en estos espacios donde plantearemos nuestros problemas, haciéndolo con un nuevo tipo de formulación, la formulación variacional o débil. Es aquí también donde buscaremos las soluciones; chocándonos contra el muro a priori infranqueable de encontrar una solución analítica y cerrada a nuestra inquietud vestida de problema variacional. Veremos que gracias a herramientas muy potentes como son los teoremas de Lax-Milgram y de Banach-Neças-Babuska podremos asegurar la existencia y unicidad de dicha solución aunque nuestra maquinaria analítica no esté capacitada para ”localizarla”. Es aquí donde nace el mundo de la discretización, y su ambición de preparar nuestro problema para encontrar una solución numérica que, de algún modo aproxime a la solución analítica que no podemos hallar. Aparecen los Métodos de Galerkin y de Petrov Galerkin y querremos ver que de la misma forma que podemos asegurar la existencia y unicidad de la solución analítica, existe y es única la mejor aproximación bajo ciertas condiciones. Estudiaremos las distintas formas que existen de discretizar un problema abstracto cuyo buen planteamiento está asegurado. Veremos que no es necesario aproximar los espacios por subespacios suyos: basta con otro tipo de espacio de Banach que de alguna forma aproxime al original. Continuando con nuestro deseo de crear modelos útiles en la Mecánica del Medio Contínuo y, en especial, en la Arquitectura, se tratarán dos ecuaciones de vital importancia en este mundo: la Ecuación de la Elasticidad y la Ecuación de Convección- Difusión. No sólo estudiaremos éstas desde el punto de vista numérico, si no que partiendo de ciertas leyes físicas basadas en el empirismo propio que a esta ciencia caracteriza, modelaremos las ecuaciones que rigen tanto el transporte de temperaturas, como la elasticidad de los cuerpos. Finalmente, se analizarán de forma muy breve varios ejemplos prácticos cuya simulación se ha realizado con el potente software de EF de FreeFem++. Como comentario final, he de decir que dado que para conseguir el objetivo último que motivó la elección de éste trabajo he tenido que realizar un largo viaje en el sentido de tocar campos muy diversos, no he buscado más que un estudio horizontal y básico (que no vertical y concreto) en este nuevo universo en el que comienzo a sumergirme: El Mundo del Modelado del Medio Continuo. Capítulo 2 Planteamiento del problema abstracto Se presenta a continuación la formulación variacional de las EDPs y resultados esenciales que nos aseguran el buen planteamiento del problema abstracto. Para su redacción, nos hemos centrado especialmente en las lecturas de [6]yde[16]. En concreto, comenzaremos ilustrando la importancia del Teorema de Lax- Milgram clásico para asegurar el buen planteamiento de un probema formulado en su modo variacional. Proposición 2.1. Sean VyWespacios de Banach, con Vespacio de Banach reflexivo. Sea a:W×V→Runa forma bilineal continua. Entonces existe un operador A∈ L(W, V 0)que verifica: hAw, viV0,V =a(w, v)∀v∈B, w ∈H. Demostración. Veamos que efectivamente Aasí definido es un operador lineal, continuo y acotado entre los espacios de Banach V0yV. Consecuencia de la linealidad de a(w, ·)en Vse tiene que, Aw es un operador lineal para cada w∈W. Por otra parte y como consecuencia de la linealidad de a(·, v)en W, se demuestra que Aes lineal. Además, como a(·,·)es continuo se tiene que existe c > 0con: hAw, viV0,V =a(w, v)< c||w||W||v||V∀w∈W, v ∈V. Si ahora dividimos por ||v||V, se obtiene que cada Aw es continuo por estar 7 CAPÍTULO 2. PLANTEAMIENTO DEL PROBLEMA ABSTRACTO 8 acotado por c||w||W: hAw, viV0,V ||v||V ≤c||w||W lo que tomando supremo sobre los v∈Vse obtiene que: ||Aw||V0≤c||w||W Finalmente, se toma supremo sobre los v∈Vy se obtiene que: ||A||L(W,V 0)≤c. Sea Vespacio de Hilbert. Considérese el siguiente problema abstracto: PV1Hallar u∈V tal que a(u, v) = f(v),∀v∈V donde a∈ L(V×V, R)yf∈ L(V, R). Recuérdese que todo espacio de Hilbert es Banach reflexivo por ser espacio normado con norma inducida por el producto escalar y por el Teorema de representación de Riesz (Teorema A.1). Se recomienda consultar Anexo. El siguiente teorema, extraído de [12], nace en la década de 1950 [16] y reúne las condiciones que aseguran el buen planteamiento del problema abstracto. Posteriormente, fue reformulado para extenderse al ámbito de los EF. Teorema 2.1. (Lax-Milgram, LM) Sea a∈ L(V×V, R)una forma bilineal, contíinua y coercitiva y f∈V0lineal y continua. Bajo estas condiciones, el problema PV1 está bien planteado. Demostración. Sea A∈ L(V, V 0)el operador lineal asociado a a(·,·), que sabemos que existe por la Proposición 2.1. Notemos por fuaAu. Por ser aelíptica, se tiene que: α||u||2 V≤a(u, u) = ||fu||V0||u||V=||Au||V0||u||V lo que en particular implica que: A está acotado inferiormente: α||u||V≤ ||Au||V0. A es inyectiva. Para verlo, recuérdese que un operador lineal es inyectivo si y sólo si ker(A) = {0}. Sea v∈ker(A), se tiene que: 0≤α||v||V≤ ||Av||V0= 0 luego, ||v||V= 0 yv= 0. Capítulo 3 Discretización En general nos encontramos con la imposibilidad de resolver problemas EDPs desde un punto de vista analítico, i.e., dar la forma de solución cerrada y exacta del problema, pese a que sea posible demostrar la existencia y unicidad de ésta como ya se ha visto con los resultados anteriores. Es por esto por lo que es necesario recurrir a otros métodos, los métodos numéricos, que den una solución aproximada del problema en cuestión. La idea de éstos reside en la transformación del problema de su forma global con dimensión infinita a una forma finita que sea resoluble mediante técnicas numéricas. Las incógnitas resultantes (llamadas grados de libertad) se trabajarán de una u otra manera según la técnica utilizada para construir una solución aproximada. Este proceso se llama discretización del problema. En el caso concreto de los EF (elementos finitos), los grados de libertad son los valores de la solución en ciertos nodos (puntos) de Ω. Para construir la solución aproximada se usa interpolación de esos valores nodales, para lo que es necesaria la construcción de una malla que se apoye en los nodos. El resultado es que podemos dar una solución a costa de cometer un error, siendo el objetivo último la estimación para el control de este error cometido. Nuestro objetivo en este capítulo es el estudio de los aspectos teóricos de la discretización del Problema abstracto, buscando en último lugar dejar la configuración discreta preparada para poder aplicar posteriormente algún método práctico de aproximación numérica. Este capítulo fue extraído de [16]. 3.1. Planteamiento del problema aproximado La idea es aproximar los espacios WyVpor espacios de dimensión finita Wh yVh, donde el subíndice hse refiere al parámetro de la discretización. Como ya 15 CAPÍTULO 3. DISCRETIZACIÓN 16 se ha dicho, dejaremos el problema preparado la para la aplicación de métodos prácticos de construcción de tales espacios, a los que asociaremos normas ||·||Wh y|| · ||Vh. Considérese, además: W(h) = W+Wh(3.1) en el que asumimos que podemos definir una norma tal que cumpla: 1. ||wh||W(h)=||wh||Wh,∀wh∈Wh. 2. ||w||W(h)≤c||w||W,∀w∈W. Observación 3.1. La construcción del espacio W(h)tiene relevancia cuando el espacio de soluciones Wno se aproxima por un subespacio suyo, si no que se busca otro espacio discreto que lo aproxime bajo ciertas condiciones y que facilite de algún modo la construcción de la solución discreta uh. El método de Petrov-Galerkin (MPG) constituye una aproximación del problema (Pa2). Puede formularse como: MPG Hallar uh∈Whtal que ah(uh, vh) = fh(vh),∀vh∈Vh donde ahyfhson las formas aproximadas de ayfrespectivamente. Un caso más particular sería la construcción de una solución sobre el problema aproximado de (Pa1), que quedaría como: MGHallar uh∈Vhtal que ah(uh, vh) = fh(vh),∀vh∈Vh. Este es el conocido como Método de Galerkin (MG). Definición 3.1. Se dirá que el problema está aproximado de forma conforme si verifica que Wh⊂WyVh⊂V. Observación 3.2. Tener una aproximación conforme permite no tener que aproximar ayfdefinida en W×VyV, por la forma aproximada ahyfh definidas respectivamente en Wh×VhyVh. Definición 3.2. Se dirá que el problema aproximado tiene la propiedad de proximidad en el caso que cumpla que ∀w∈W, ∃wh∈Whtal que l´ım h→0||w−wh||W(h)→0, es decir, el espacio Whaproxima bien a W. Definición 3.3. Se dirá que MPG es consistente si ahpuede ser extendida a W(h)×Vh(es decir, tiene sentido hablar de ah(u, vh)) y si la solución exacta uverifica el problema aproximado, i.e.: ah(u, vh) = fh(vh)∀vh∈Vh. CAPÍTULO 3. DISCRETIZACIÓN 17 Se dirá que MPG es asintóticamente consistente si ahestá acotado de manera uniforme en Wh×Vh, y ∃Πh:W→Whun operador lineal y acotado, tal que: ∀w∈W, ||Πhw−w||W(h)≤c´ınf wh∈Wh ||w−wh||W(h) y l´ım h→0sup vh∈Vh |fh(vh)−ah(Πhu, vh)| ||vh||Vh= 0. Esto implicará que la EDP está bien aproximada de forma asintótica. Definición 3.4. Se define el error de consistencia,Rh, como: Rh(u) = sup vh∈Vh |fh(vh)−ah(Πhu, vh)| ||vh||Vh (3.2) Lema 3.1. El error de consistencia es independiente al operador Πhusado, siempre que se verifique la propiedad de proximidad. Demostración. Sea Π1 hun operador según la definición anterior y Π2 hotro operador tal que satisface que ∀w∈W, ||Π2 hw−w||W(h)≤c´ınfwh∈Wh||w−wh||W(h). Entonces, ∀vh |fh(vh)−ah(Π2 hu, vh)|=|fh(vh)−ah(Π2 hu, vh)±ah(Π1 hu, vh)| ≤ |fh(vh)−ah(Π1 hu, vh)|+|ah(Π1 hu, vh)−ah(Π2 hu, vh)|=|fh(vh)−ah(Π1 hu, vh)|+ |ah(Π1 hu−Π2 hu, vh)|=⇒R2h≤R1h+c||ah||Wh,Vh´ınf wh∈Wh ||u−wh||W(h), de donde tomando el límite h→0y usando la propiedad de proximidad, se tiene que R2h→0y, por tanto, la consistencia asintótica mediante el operador Π2 h. Lema 3.2. (Ortogonalidad de Galerkin) Si el problema aproximado es consistente se tiene que: ∀vh∈Vh, ah(u−uh, vh)=0.(3.3) Demostración. Dado que es consistente, ah(u, vh) = fh(vh). Por su parte la solución del problema aproximado verifica que ah(uh, vh) = fh(vh). Restando ambas expresiones se tiene el resultado. CAPÍTULO 3. DISCRETIZACIÓN 18 3.2. Cálculo efectivo de la solución discreta A continuación, veremos que el problema aproximado es equivalente a un sistema lineal mediante el método de Petrov-Galerkin. Para ello, denotamos M= dim WhyN= dimVh,BW={ψ1, ..., ψM}una base de WhyBV={ϕ1, ..., ϕN} una base de Vh. La solución uhpuede, por tanto, escribirse como combinación de los elementos de BW, i.e.: uh= M X j=1 Ujψj siendo U= (Uj)1≤j≤M∈RMel vector de coordenadas de uhreferidos a la base BW. Por otra parte, considérese la matriz A ∈ RN,M , definida por: Aij =ah(ψj, ϕi) 1 ≤i≤N1≤j≤M, (3.4) conocida como Matriz de Rigidez; y sea, F ∈ RN, tal que: Fi=fh(φi),1≤i≤N. (3.5) Finalmente vemos que, efectivamente, el problema (MPG) es equivalente a un problema lineal: AU=F,(3.6) donde las incógnitas son las coordenadas de uhen BM. Observación 3.3. Basta que asea simétrica para que Atambién lo sea. Si el problema aproximado es coercitivo, consistente y conforme, la matriz de rigidez Aes definida positiva. 3.3. Interpretación geométrica del método de Galerkin Cuando la forma bilineal a(·,·)es continua, simétrica y coercitiva, la solución uh del problema aproximado MG es la proyección ortogonal de uen Vhcon respecto al producto escalar inducido por a(·,·); producto escalar que es equivalente al producto escalar canónico de V. Dado u∈V, su proyección ortogonal uh∈Vhestá caracterizada por las siguientes propiedades equivalentes: Distancia mínima entre uyuh: a(u−uh, u −uh) = m´ın vh∈Vh a(u−vh, u −vh). CAPÍTULO 3. DISCRETIZACIÓN 19 Ortogonalidad: a(u−uh, vh)=0 ∀vh∈Vh. Gracias a esto puede comprenderse mejor la estimación de error a priori que se realizará próximamente como consecuencia del Lema de Lax- Milgram, en el sentido de que uhes el elemento de Vhmás próximo a u en la norma inducida por la forma bilineal a(·,·), y puede utilizarse en la práctica para obtener estimaciones del error que se comete en la aproximación y controlarlo de algún modo, ya que si se es capaz de construir una aproximación cómoda ¯uhtal que no minimice distancia, sabemos que el error cometido al utilizar uhserá menor que el que se comete al usar la aproximación cómoda ¯uh. 3.4. Buen planteamiento del problema El objetivo ahora es estudiar cuándo tenemos un buen planteamiento del problema aproximado. El caso particular de tener una aproximación conforme, coercitiva y consistente es relativamente inmediato ya que teniendo en cuenta que Vh es un subespacio finito de V,Vhes cerrado. Por ser un subespacio cerrado de un espacio de Hilbert, es, a su vez, espacio de Hilbert. Por tanto, consecuencia inmediata del Lema de Lax-Milgram sobre Vhes la siguiente proposición. Proposición 3.1. Sea Vespacio de Hilbert, aforma bilineal y continua. Si Vhes un espacio finito-dimensional y se tiene que aes coercitiva en V Vh⊂V. Entonces el problema aproximado MG está bien planteado en el sentido de Hadamard, i.e. existe una única solución del sistema de ecuaciones y se puede hacer una primera estimación del error cometido en la aproximación: ||u−uh||V≤M α´ınf vh∈Vh ||u−vh||V.(3.7) Estudiamos ahora el caso más general, donde las condiciones del Teorema BNB en el caso del Problema MPG pueden ser traducidas como: BNB1(h) ∃αh>0,´ınfwh∈Whsupvh∈Vh ah(wh, vh) ||wh||Wh||vh||Vh > αh BNB2(h) ∀vh∈Vh,∀wh∈Wh, ah(wh, vh) = 0 =⇒vh= 0. Proposición 3.2. Se tiene que: 1. BNB1(h)⇐⇒ Ker(A) = {0}. 2. BNB2(h)⇐⇒ rg(A) = dimVh. CAPÍTULO 3. DISCRETIZACIÓN 20 3. Si dim Vh= dim Wh,BNB1(h)⇐⇒ BNB2(h). Demostración. 1. Veamos que BNB1(h) =⇒Ker(A) = {0}. Tomemos Para ello sea ξ=PM j=1 Xjψjy comencemos obteniendo una equivalencia a que X∈Ker(A). X∈Ker(A)⇐⇒ (∀i∈ {1, ..., N},PM j=1 Ai,jXj= 0) ⇐⇒ (∀i, ah(ξ, ϕi) = 0) ⇐⇒ supv∈Vhah(ξ, v)=0. Ahora bien, BNB1(h) hace que si se tiene que (supv∈Vhah(ξ, v) = 0) entonces ξ= 0 y, por tanto, Ker(A) = {0}. Por otra parte, veamos que Ker(A) = {0}=⇒BNB1(h). Para realizar la demostración vamos a razonar al absurdo. Para ello, asumamos que Ker(A) = {0}y supongamos que, sin embargo, ∃whn ∈Whsubsucesión, tal que kwhnk= 1,l´ım n→∞ whn =why que cumple: sup v∈Vh ah(whn, v) kvkVh ≤1 n. Dado que la esfera unidad es un subespacio compacto de Wh, existe una subsucesión que seguiremos denotando como whn, que converge a whcuando n→ ∞. Dicho límite whha de cumplir que kwhkWh= 1 ysup v∈Vh = 0. Esto implica que X∈Ker(A)donde Xes el vector de coordenadas de wh. Dado que Ker(A) = {0}se tiene que X= 0, lo que contradice que kwhkWh= 1. 2. La prueba de 2es similar a la de 1sin más que aplicar lo mismo a AT. 3. Consecuencia directa de 1,2y el Teorema del Rango (Teorema 3.1). Finalmente, el siguiente teorema expresa las condiciones necesarias para asegurar el buen planteamiento del problema. Teorema 3.1. (Teorema del Rango) Asúmase que: ahestá acotada en Wh×Vhyfhes contínua en Vh. Se cumple BNB1(h). VhyWhtienen igual dimensión. CAPÍTULO 3. DISCRETIZACIÓN 21 En tal caso, el problema aproximado está bien planteado en el sentido de Hadamard y se verifica la estimación a priori: ||uh||Wh≤1 αh ||fh||V0 h. Siendo la demostración inmediata a partir del tercer punto de la Proposición 3.2 y el Teorema BNB (Teorema 2.5). 3.5. Análisis de Error Queremos ahora hallar el error uh−upara diferentes configuraciones, donde uh resuelve el problema aproximado y uel problema exacto. Caso General A continuación presentamos un teorema para estimar el error de la situación más general. Se exigen las cuatro propiedades que se necesitan para probar la convergencia de la solución aproximada: estabilidad uniforme, acotación uniforme, consistencia asintótica y que se verifique la propiedad de proximidad. Existe un principio, conocido como Principio de Lax que establece que estabilidad y consistencia implica convergencia. No obstante, dicho principio no implica ni continuidad ni que se verifique la propiedad de proximidad. Teorema 3.2. Supóngase que se verifica lo siguiente: 1. La condición BNB1(h)se satisface uniformemente con h, y dim(Wh) = dim(Vh). 2. La forma bilinear ahestá uniformemente acotada en Wh×Vh. 3. El problema aproximado es asintóticamente consistente 4. El problema aproximado tiene la propiedad de proximidad. Entonces, denotando como Rhal error de consistencia, se verifica que: ||u−uh||W(h)≤1 αh Rh(u) + c´ınf wh∈Wh ||u−wh||W(h)(3.8) y l´ım h→0||u−uh||W(h)= 0.(3.9) CAPÍTULO 3. DISCRETIZACIÓN 22 Demostración. Dado que ahestá acotado en Wh×Vh, utilizamos el operador proyector Πh:W→Whpreviamente usado para la consistencia asintótica. Por linealidad tenemos que: ah(uh−Πhu, vh) = fh(vh)−ah(Πh, vh)(3.10) Por otra parte, aplicando la condición BNB1(h)auh−Πhu(obsérvese que uh−Πhu∈Wh), implica que αh||uh−Πhu||Wh≤Rh(u). Y dado que las normas || · ||Why|| · ||W(h)son iguales en Wh, por la desigualdad triangular y la consistencia asintótica, se tiene que: ||uh−u||W(h)≤1 αh Rh(u) + ||u−Πhu||W(h)≤ 1 αh Rh(u) + c´ınf wh∈Wh ||w−wh||W(h), quedando probado la estimación (3.8). Finalmente, haciendo uso de la propiedad de proximidad y de la consistencia asintótica, se llega a (3.9): l´ım h→0||u−uh||W(h)= 0. Casos particulares Entre los casos particulares que vamos a estudiar, comencemos tratando el caso de tener una aproximación no consistente y no conforme. Asumimos que la forma bilineal ah(·,·)puede ser extendida a W(h)×Vh, tal que ah(w, vh) tenga sentido para w∈Wyvh∈Vh. El siguiente resultado se conoce como Segundo Lema de Strang [16]. Lema 3.3. (Strang 2) Asúmase lo siguiente: 1. Se cumple la condición BNB1(h) y dim(Wh) = dim(Vh). 2. La forma bilineal aproximada ahestá acotada en W(h)×Vh. Se tiene entonces la siguiente estimación de error: ku−uhkW(h)≤1 + kahkW(h),Vh αh´ınfwh∈Whku−whkW(h) +1 αhsupvh∈Vh |fh(vh)−ah(u,vh)| kvhkVh . Demostración. Sea wh∈Wh. Entonces, ah(uh−wh, vh) = ah(uh±u−wh, vh) = ah(uh−u, vh) + ah(u−wh, vh) =fh(vh)−ah(u, vh) + ah(u−wh, vh). CAPÍTULO 3. DISCRETIZACIÓN 23 La condición BNB1(h) implica que αhkuh−whkW(h)≤sup vh∈Vh |fh(vh)−ah(u, vh)| kvhkVh +kahkW(h),Vhku−whkW(h). Aplicamos la desigualdad triangular al primer miembro para obtener que: αhkuh−whkW(h)=αhkuh±u−whkW(h)≤αhku−uhkW(h)+αhku−whkW(h). Por último, reagrupamos y obtenemos la estimación de error para una aproximación no conforme y no consistente: ku−uhkW(h)≤1 + kahkW(h),Vh αh´ınfwh∈Whku−whkW(h) +1 αhsupvh∈Vh |fh(vh)−ah(u,vh)| kvhkVh . Observación 3.4. Recuérdese que cuando tenemos una configuración consistente la solución real satisface la versión aproximada de la forma bilineal, i.e. ah(u, vh) = fh(vh). Es por esto por lo que el último término de la estimación anterior se anula, quedando la expresión del error para el caso de la aproximación no conforme y consistente de la forma siguiente: ku−uhkW(h)≤ 1 + kahkW(h),Vh αh!´ınf wh∈Wh ku−whkW(h). Estudiamos ahora el caso de la aproximación no consistente y conforme. Por ser consistente, se tiene que Wh⊂WyVh⊂V. Como consecuencia, es evidente que W=W(h). Esta vez, sin embargo, no asumimos que k·kW(h)y k · kWsean iguales. Como diferencia de la sección anterior, ah(·,·)no puede ser extendida a W×Vh. Tenemos el siguiente resultado [16]: Lema 3.4. (Strang 1) Asúmase lo siguiente: 1. Wh⊂WyVh⊂V. 2. Se cumple la condición BNB1(h) y dim(Wh) = dim(Vh). 3. La forma bilineal ahestá acotada en Wh×Vh, y aestá acotada en W×Vh cuando Westá equipada con la norma k·kWh. Se tiene entonces la siguiente estimación de error: ku−uhkW(h)≤1 αhsupvh∈Vh |f(vh)−fh(vh)| kvhkVh +´ınfwh∈Whh1 + kakW(h),Vh αhku−whkW(h)+1 αhsupvh∈Vh |a(wh,vh)−ah(wh,vh)| kvhkVhi CAPÍTULO 3. DISCRETIZACIÓN 24 Demostración. Sea wh∈W(h). La condición BNB1(h) implica que: αhkuh−whkW(h)≤sup vh∈Vh ah(uh−wh, vh) kvhkVh . Obsérvese que: a(u−wh, vh)−f(vh) + a(wh, vh)=0. Sumando la observación anterior a ah(uh−wh, vh)se llega a que: ah(uh−wh, vh) = fh(vh)−ah(wh, vh) = a(u−wh, vh) + a(wh, vh)−ah(wh, vh) + fh(vh)−f(vh). Obsérvese, además que aplicando la desigualdad triangular a kuh−whkW(h), se tiene que: αhkuh−whkW(h)≤αhku−uhkW(h)+αhku−whkW(h) Finalmente, tenemos la estimación que buscábamos: ku−uhkW(h)≤1 αhsupvh∈Vh |f(vh)−fh(vh)| kvhkVh +´ınfwh∈Whh1 + kakW(h),Vh αhku−whkW(h)+1 αhsupvh∈Vh |a(wh,vh)−ah(wh,vh)| kvhkVhi. El último caso que analizamos es el de la aproximación consistente y conforme, con ah=ayfh=f, recordemos que el problema aproximado es: Halla uh∈Whtal que a(uh, vh) = f(vh),∀vh∈Vh(3.11) El siguiente resultado se conoce como Lema de Cèa y fue extraido de [20]. Lema 3.5. (Lema de Cèa) Asúmanse las hipótesis siguientes: 1. Se cumple la condición BNB1(h) y dim(Wh) = dim(Vh). 2. La forma bilineal aproximada ahestá acotada en W(h)×Vh con Vh⊂VyWh⊂W,ah=ayfh=f. Sea uhuna solución aproximada del problema (3.11). Se tiene la siguiente estimación de error: ku−uhkW≤1 + kakW,V αh´ınf wh∈Wh ku−whkW.(3.12) Demostración. Sea wh∈Wh. La ortogonalidad de Galerkin implica que: ∀vh∈Vh, a (uh−wh, vh) = a(u−wh, vh). CAPÍTULO 4. ECUACIÓN DE ELASTICIDAD 31 implicando que el Primer Tensor de Piola-Kirchoff satisface la siguiente ecuación diferencial: −div T(x) = f(x)para todo x∈Ω.(4.8) Usando ahora la identidad: ∇Φ(x)−1T(x) = ∇Φ(x)−1[det∇Φ(x)˜ T(Φ(x))]∇Φ(x)−T podemos deducir de la simetría de ˜ T(˜x)que la matriz ∇Φ(x)−1T(x)es también simétrica. Es entonces claro ver que las ecuaciones establecidas en el Teorema 4.1.1 son equivalentes a las siguientes ecuaciones definidas sobre la configuración de referencia: −div T(x) = f(x)para todo x∈Ω T(x)n(x) = g(x)para todo x∈Γ1 ∇Φ(x)−1T(x)∈S3para todo x∈Ω (4.9a) (4.9b) (4.9c) donde el subconjunto Γ1de ∂Ωestá definido por ˜ Γ1= Φ(Γ1). Finalmente, definimos como el Segundo Tensor de Piola-Kirchoff,Σ : Ω → S3, a: Σ(x) := ∇Φ(x)−1T(x)(4.10) para todo x∈Ω. Y podemos reescribir por última vez ecuaciones de equilibrio en la configuración de referencia en términos de Σcomo: −div(∇Φ(x)Σ(x)) = f(x)para todo x∈Ω (∇Φ(x)Σ(x))n(x) = g(x)para todo x∈Γ1. (4.11a) (4.11b) Las incógnitas en cualquiera de los sistemas de ecuaciones de equilibrio obtenidos son la deformación del cuerpo y los tensores de Piola-Kirchoff. Con el objetivo de determinar estas incógnitas, ambos sistemas hay que complementarlos con una ecuación que relacione los campos. Este es el objetivo de la siguiente sección. 4.1.2. Relaciones constitutivas de materiales elásticos Está claro que el tensor de esfuerzo tiene que depender de la deformación inducida en el cuerpo por las fuerzas que actúan sobre él. Esta dependencia se establece mediante la ecuación constitutiva del material, por medio de una función de respuesta específica del material considerado. En lo que sigue, consideraremos un tipo de estos materiales, los materiales elásticos. Definición 4.2. Se dice que un material es elástico si existe una función T#:¯ Ω×M3 +→M3tal que: T(x) = T#(x, ∇Φ(x)) ∀x∈Ω. CAPÍTULO 4. ECUACIÓN DE ELASTICIDAD 32 O, de forma equivalente, si existe una función Σ#: Ω ×M3 +→S3tal que : Σ(x)=Σ#(x, ∇Φ(x)) ∀x∈Ω. Tanto Σ#como T#se conocen como función respuesta del material. Una función respuesta no puede ser arbitraria. En Física, existe una axioma que establece que un observable ha de ser invariante con respecto a transformaciones ortogonales del sistema de referencia. En el caso particular de materiales elásticos, el observable es el campo de tensiones ¯ t, que por lo explicado, ha de ser invariante frente a transformaciones ortogonales de la base del sistema de referencia utilizado para representarlo. Este es el conocido Axioma de la Invarianza de marco de referencia para sistemas materiales. En el siguiente teorema damos formalismo matemático en términos de la función de respuesta. Teorema 4.2. Un material elástico satisface el Axioma de la Invarianza de marco de referencia para sistemas materiales si y solamente si T#(x, QF) = QT#(x, F)x∈ΩyQ∈O3 +yF∈M3 + o, equivalentemente satisface que: Σ#(x, QF)=Σ#(x, F)x∈ΩyQ∈O3 +yF∈M3 +. Lo anterior indica que la función de respuesta Σ#depende de Fsólamente vía la matriz simétrica definida positiva U:= (FTF)1 2. Para verlo, realizamos una factorización polar F=RU, con R:= FU−1∈O3 +y vemos que el Axioma de la Invarianza de marco de referencia para Σ#puede expresarse como: Σ#(x, F)=Σ#(x, U)∀x∈ΩyF=RU ∈M3 +. Lo anterior implica que el Segundo Tensor de Piola Σ : ¯ Ω :→S3depende del campo de deformaciones Φ : ¯ Ω→R3exclusivamente a través del tensor métrico asociado C:= ∇ΦT∇Φ, es decir: Σ(x) = ˜ Σ(x, C(x)) ∀x∈¯ Ω,(4.12) y la función ˜ Σ : ¯ Ω×S3 >→S3 >está definido por ˜ Σ(x, C)=Σ#(x, C 1 2)∀x∈¯ ΩyC∈S3 >. Definición 4.3. Se dice que un material elástico es isotrópico en un punto x∈Ωsi la respuesta producida es igual en todas sus direcciones, i.e. el tensor de tensiones de Cauchy ˜ T(x)es invariante frente a rotaciones O3 +del sistema de referencia en torno al punto x. Nos restringimos a materiales isotrópicos, que son aquellos materiales que son isotrópicos (Definición 4.3) en todos sus puntos. El siguiente teorema ofrece una caracterización de cúando la respuesta de un material elástico es isotrópica. CAPÍTULO 4. ECUACIÓN DE ELASTICIDAD 33 Teorema 4.3. Un material elástico ocupando una configuración ¯ Ωes isotrópico si y solamente si: T#(x, FQ) = T#(x, F)Q∀x∈ΩyQ∈O3 +yF∈M3 +, o equivalentemente, Σ#(x, FQ) = QTΣ#(x, F)Qpara todo x∈ΩyQ∈O3 +yF∈M3 +. Otra propiedad que pueden satisfacer los material elásticos es la homogeneidad. Definición 4.4. Un material elástico ocupando una configuración ¯ Ωes homogéneo si su función respuesta es independiente al punto x∈¯ Ωconsiderado. Esto significa que la función respuesta T#oΣ#no depende de su primera variable, i.e., existen aplicaciones a las que se seguirán denotando respectivamente T#yΣ#tales que: T#(x, F) = T#(F)∀x∈ΩyF∈M3 + y Σ#(x, F)=Σ#(F)∀ ∈ ΩyF∈M3 +. Definición 4.5. Una configuración ¯ Ωse dice que está en estado natural o libre de tensiones si verifica que: T#(x, I)=0 ∀x∈Ω o equivalentemente, Σ#(x, I)=0 ∀x∈Ω. Hemos visto en la expresión (4.12) que el Segundo Tensor de Piola-Kirchoff Σ : Ω→S3puede expresarse en función del campo de deformaciones de la forma: Σ(x) = ˜ Σ(x, C(x)),donde C(x) = ∇ΦT(x)∇Φ(x)para todo x∈Ω. El Teorema de Rivlin-Ericksen establece que si el material es isotrópico, podemos simplificar notablemente la dependencia de Σcon respecto a C(x). Veámoslo. Teorema 4.4. (Teorema de Rivlin-Ericksen) Si un material elástico es isotrópico y satisface el principio de invarianza frente a transformaciones del sistema de referencia, entonces existen funciones γ# i: Ω ×R3→Rcon i∈ {1,2,3}tales que: Σ(x) = γ0(x)I+γ1(x)C(x) + γ2(x)C2(x)for all x∈Ω donde γi(x) = γ# i(x, trC, tr(Cof C),det C). Demostración. La idea principal de la prueba reside en que toda matriz verifica su polinomio característico. No incluimos la demostración debido a su extensión, pero para el lector interesado consultar [13] o [14]. CAPÍTULO 4. ECUACIÓN DE ELASTICIDAD 34 Observación 4.1. Nótese que trC(x),tr(cof C(x)),ydetC(x)son los invariantes principales de la matriz C(x). Pese a que el teorema anterior ha reducido notablemente la cantidad de respuestas posibles que puede dar un sistema elástico e isotrópico que satisfazca la invarianza frente a cambios de sistema de referencia, aún las posibles funciones respuesta pertenecen a un conjunto demasiado general y extenso como para que pueda trabajarse con él. Es por esto por lo que ahora nos centramos en un tipo especial de deformaciones, aquellas cuyo estado deformado no dista en exceso de la configuración de referencia. En términos del campo de desplazamiento u: Ω →R3, el tensor métrico C cumple que: C(x) = I+ 2E(x),(4.13) donde E(x)es conocido como el Tensor de tensiones de Green-St Vernant en el punto x∈Ωy en función del campo de desplazamientos, puede escribirse como: E(x) := 1 2∇uT(x) + ∇u(x) + ∇uT(x)∇u(x)(4.14) Hemos supuesto que la deformación de la configuración de referencia es pequeña, por lo que los tensores E(x)han de ser ”pequeños” para cualquier x∈Ω. Es por esto por lo que podemos hacer un desarrollo de Taylor de la función respuesta dada por el Teorema de Rivlin-Ericksen (4.4): trC(x) = 3 + 2trE(x) tr(Cof C(x)) = 3 + 4trE(x) + o(|E(x)|) detC(x) = 1 + 2trE(x) + o(|E(x)|) C2(x) = 1 + 4E(x) + o(|E(x)|) que asumiendo que las funciones γ# ison lo suficientemente regulares, se tiene que: Σ(x)=Σ#(x, I) + {(λ(x)trE(x))I+ 2µ(x)E(x)}+ox(|E(x)|)(4.15) Observación 4.2. Tanto µ(x)como λ(x)son funciones independientes al campo de desplazamientos u(x). Si además el material es homogéneo, entonces µyλson constantes. Si Ωes un estado natural o libre de tensiones del sistema, entonces el candidato obvio para poder establecer una relación constitutiva es: Σ(x) = λ(trE(x))I+ 2µE(x)(4.16) para todo x∈Ω. En este caso, µyλse conocen como las constantes de Lamé del material. CAPÍTULO 4. ECUACIÓN DE ELASTICIDAD 35 Definición 4.6. Se dirá que un material es del tipo St Venant-Kirchhoff si tiene como relación constitutiva la ofrecida en la Ecuación (4.16). Observación 4.3. Las constantes de Lamé se determinan experimentalmente. Dado que nuestro análisis por ahora tiene un fin puramente teórico, no estamos interesados en sus valores exactos, sino en el hecho de saber que son positivas. En numerosas ocasiones las constantes de Lamé aparecen escritas en función de el coeficiente de Poisson νy el módulo de Young Emediante las expresiones: ν=λ 2(λ+µ)yE=µ(3λ+ 2µ) λ+µ. 4.1.3. Las ecuaciones de Elasticidad no Lineal Ahora introduciremos las ecuaciones de la elasticidad no lineal en un problema de condiciones de contorno. El motivo de introducirlas y, sin embargo, hacerlo de una forma tan breve, es porque las necesitaremos en la siguiente sección para deducir las ecuaciones de elasticidad lineal. Como es natural por todo lo expuesto hasta el momento, nacen de la combinación de las ecuaciones de equilibrio previamente obtenidas junto con la relación constitutiva del material considerado. Además le añadiremos condiciones de contorno sobre Γ0:= ∂Ω\Γ1. Asumiendo que el material ofrece una función respuesta conocida, dada por T# o por Σ#y que está fijado en Γ0, concluimos que la deformación inducida por las fuerzas que se aplican con densidades fygviene dada por el problema no lineal de condiciones de contorno siguiente: −div T(x) = f(x), x ∈Ω Φ(x) = x, x ∈Γ0 T(x)n(x) = g(x), x ∈Γ1 donde T(x) = T#(x, ∇Φ(x)) = ∇+ Φ(x)Σ#(x, ∇Φ(x)),para todo x∈Ω.(4.17) 4.1.4. Las ecuaciones de Elasticidad Lineal Las ecuaciones de Elasiticidad lineal se obtienen a partir de las de elasticidad no lineal, bajo la hipótesis de que esta vez el cuerpo sufre una deformación pequeña con respecto a su configuración de referencia, es decir que: Φ(x) = x+u(x)con |∇u(x)|  1para todo x∈Ω. Es por tanto razonable aproximar la parte no lineal del modelo anterior por su parte lineal, obtenida de sustituir la función de respuesta de la Ecuación (4.17) por la parte afín de uen la aproximación por series de Taylor: T#(x, ∇Φ(x)) = T#(x, I) + ∂T# ∂F (x, I)∇u(x) + o(|∇u(x)|). CAPÍTULO 4. ECUACIÓN DE ELASTICIDAD 36 Así, el desplazamiento en respuesta a las fuerzas aplicadas satisface el problema de contorno lineal siguiente: −div T#(·, I) + ∂T # ∂F (·, I)∇u=fen Ω Tn =gsobre Γ1 u= 0 sobre Γ0. Estas son las Ecuaciones de Elasticidad Lineal. Un caso particularemente importante es aquel en el que el cuerpo está constituido de un material elástico, isotrópico y homogéneo tal que su configuración de referencia constituya un estado natural. Su ecuación constitutiva será entonces: Σ(x) = λ(trE(x))I+ 2µE(x) + o(|E(x)|)(4.18) con λ, µ > 0las constantes de Lamé del material. Se tiene entonces que ∀x∈Ω, E(x) = 1 2∇ΦT(x)∇Φ(x)−I=1 2∇uT(x) + ∇u(x)+ox(|∇u(x)|)(4.19) y T(x) = ∇Φ(x)Σ(x)=(I+∇u(x))(λ(trE(x))I+ 2µE(x)) =λ 2tr∇uT(x) + ∇u(x)+µ∇uT(x) + ∇u(x)+ox(|∇u(x)|). Finalmente, Ecuaciones de Elasticidad Lineal puesden escribirse como: −div σ(x) = f(x), x ∈Ω u(x)=0, x ∈Γ0 σ(x)n(x) = g(x), x ∈Γ1 (4.20a) (4.20b) (4.20c) donde σ(x) = λ(tre(x))I+ 2µe(x) e(x) = 1 2∇uT(x) + ∇u(x). (4.21a) (4.21b) En la siguiente sección se probará que este sistema lineal de EDPs tiene solución única en espacios de Sobolev. 4.2. Elasticidad Lineal 4.2.1. Desigualdades de Korn A continuación, vamos a estudiar la existencia y unicidad de solución cuando tratamos con la ecuación ya presentada que modela el comportamiento linealelástico que sufre un material. Mediante el lema de J.L.Lions podrá establecerse CAPÍTULO 4. ECUACIÓN DE ELASTICIDAD 37 la desigualdad fundamental de Korn que implicará que las ecuaciones de elasticidad lineal tienen una sulución lineal, estable y regular siempre que los datos para el contexto del problema cumplan ciertas propiedades de regularidad. El contenido correspondiente a esta sección se ha extraído de [7]. Observación 4.4. Si v∈L2(Ω) entonces v∈H−1(Ω) y∂iv∈H−1(Ω),1≤ i≤N. Para demostrarlo basta ver que dado ψ∈D(Ω): | hv, ψi | =ZΩ vψ≤ ||v||L2(Ω)||ψ||H1(Ω) | h∂iv, ψi | =−ZΩ v∂iψ≤ ||v||L2(Ω)||∂iψ||L2(Ω) ≤ ||v||L2(Ω)||ψ||H1(Ω) La implicación contraria recibe el nombre del primer matemático capaz de probarla, J.L.Lions. Lema 4.1. (J.L.Lions) Sea v∈D0(Ω). Se tiene que: Si v, ∂iv∈H−1(Ω),1≤i≤Nentonces v∈L2(Ω). Teorema 4.5. (Primera Desigualdad de Korn) Sea Γ0⊂Γtal que la medida de Γ0>0. Entonces, existe una constante Ctal que: ||e(v)||L2(Ω;S3)≥C||v||H1(Ω;R3),∀v∈H1 Γ0(Ω;R3), donde H1 Γ0(Ω;R3)se refiere al espacio v∈H1(Ω;R3); v= 0 en Γ0ye(v)a la matriz con eij(v) := 1 2(∂ivj+∂jvi),∀i, j ∈ {1,2,3}. Demostración. Para la demostración procedemos en los pasos siguientes: 1. Teniendo en cuenta que eij(v) := 1 2(∂ivj+∂jvi), se tiene que la desigualdad de Korn es una consecuencia de ∂ijvk=∂iejk(v)+∂jeik(v)−∂keij(v). Si v∈L2(Ω;R3)ye(v)∈L2(Ω;S3)la Observación 4.4 muestra que ∂ijvk∈H−1(Ω). Dado que las funciones ∂jvktambién pertenecen al espacio H−1(Ω) y el Lema de J.L.Lions (Lema 4.1) muestra que ∂jvk∈L2(Ω). Lo que implica que el espacio: E(Ω;R3) := v∈L2(Ω;R3); e(v)∈L2(Ω;S3) es justamente el espacio H1(Ω; R3). 2. El espacio E(Ω; R3)es de Hilbert equipado con la norma: ||v||E(Ω;R3):= kvkL2(Ω;R3+ke(v)kL2(Ω;S3(4.22) de igual manera que lo es H1(Ω) con su norma ||·||H1(Ω). Recordemos que esto es obvio ya que ya vimos anteriormente que los espacios Hk(Ω) son CAPÍTULO 4. ECUACIÓN DE ELASTICIDAD 38 espacios de Hilbert. Por otra parte, la aplicación: id :v∈H1(Ω;R3)7→ v∈E(Ω;R3)) es claramente lineal, biyectiva y continua, por lo que el Teorema de la Aplicación Abierta (Teorema 2.3), muestra que la aplicación id así definida es abierta. Existe, por tanto, C > 0tal que: ||v||H1(Ω;R3)≤C||v||E(Ω;R3),∀v∈E(Ω;R3) lo que es equivalente a escribir que: ||v||L2(Ω;R3)+||e(v)||L2(Ω;R3)≥C−1||v||H1(Ω;R3),∀v∈H1(Ω;R3). 3. Establecemos además, que si v∈H1 Γ0(Ω;R3)satisface que e(v)=0, entonces v= 0. Esto es una consecuencia del Punto 1, que muestra que si v∈H1 Γ0(Ω;R3)tal que satisface e(v)=0entonces ha de satisfacer también: ∂ijvk= 0 en Ω y por tanto vtiene que ser afín [7], es decir, de la forma v(x) = b+Ax con b∈R3yA∈M3y dado que la parte simétrica del gradiente de v (que es justamente e(v)), se anula en Ωla matriz Aha de ser además antisimétrica. Dado que el rango de una matriz no nula y antisimétrica de orden 3, es necesariamente 2, el lugar geométrico de los puntos que satisfacen b+Ax = 0 puede ser o una recta en R3o el vacío, dependiendo para cada caso de si el sistema lineal tiene o no soluciones. Pero b+Ax = 0 en Γ0que tiene medida positiva, entonces b= 0 yA= 0, lo que hace que v= 0 en Ω. 4. Finalmente, obtenemos la desigualdad del enunciado por contradicción. Supongamos que es falsa. Existe una sucesión {vn}n∈N⊂H1 Γ0(Ω;R3)tal que: ||vn||H1(Ω;R3)= 1 ∀n∈N ||e(vn)||L2(Ω;S3)→0cuando n→ ∞. Dado que Ωes un dominio, la inclusión H1(Ω;R3)⊂L2(Ω;R3)es compacta por el teorema de Rellich-Kondrasov (Teorema A.9). La sucesión {vn}n∈Npor estar acotada en H1(Ω;R3)contiene una subsucesión vσ(n)n∈N que converge en L2(Ω;R3)y donde σ:N→Nes una función creciente. Dado que la sucesiones vσ(n)n∈Nye(vσ(n))n∈Nconvergen respectivamente en los espacios L2(Ω; R3)yL2(Ω; S3)son sucesiones de Cauchy. Entonces la sucesión vσ(n)n∈Nes una sucesión con respecto de la norma || · ||E(Ω;R3), así que también lo es con respecto a || · ||H1(Ω;R3)por la desigualdad establecida en el Punto 2. El espacio H1 Γ0(Ω;R3)es completo al CAPÍTULO 4. ECUACIÓN DE ELASTICIDAD 39 ser subespacio cerrado de H1(Ω;R3), entonces existe v∈H1 Γ0(Ω;R3)tal que: vσ(n)→v en H1(Ω;R3) Dado que ese límite satisface e(v) = l´ım n→∞ e(vσ(n))=0 de lo que sigue que v= 0 por el Punto 2. Pero esto contradice la relación ||v||H1(Ω;R3)= l´ımn→∞ ||vσ(n)||H1(Ω;R3)= 1. Luego queda probada la desigualdad de Korn. La desigualdad establecida en la Parte 2 de la prueba se llama Desigualdad de Korn sin condiciones de contorno. El resultado de unicidad correspondiente a la Parte 3 se llama Lema del Desplazamiento Rígido Infinitesimal, es decir, si v∈H1(Ω;R3)es tal que e(v)=0,vtiene que ser necesariamente: v(x) = a+b×x∀x∈Ωya, b ∈R3. Esto hace que el espacio de todos los desplazamientos rígidos infinitesimales tenga dimensión finita y se denota como: Rig Ω;R3:= w: Ω →R3;w(x) = a+b∧x, a, b ∈R3. Observación 4.5. En el caso especial de que Γ0=∂Ω, la Primera Desigualdad de Korn (Teorema 4.5) es una consecuencia trival de la identidad: ZΩ |e(v)|2dx =ZΩ |∇v|2dx ∀v∈H1 0(Ω;R3) sin más que integrar por partes dos veces. Teorema 4.6. (Segunda Desigualdad de Korn) Para todo v∈H1Ω; R3, existe C > 0tal que: ke(v)kL2(Ω;S3)≥C´ınf w∈Rig(Ω;R3)kv+wkH1(Ω;R3). Demostración. La prueba de esta desigualdad se obtiene del Teorema (4.5) con las Partes 2 y 3 adaptadas como sigue. Sea v∈H1Ω;R3tal que satisface e(v)=0necesariamente pertenece al espacio Rig Ω; R3. La sucesión (vn)n∈N se define ahora en H1Ω;R3y satisface: ´ınf w∈Rig(Ω;R3)kvnkH1(Ω;R3)= 1 ∀n∈N CAPÍTULO 4. ECUACIÓN DE ELASTICIDAD 40 ke(vn)kL2(Ω;S3)→0cuando n→ ∞. Si Rig Ω;R3⊥denota el complemento ortogonal de Rig Ω; R3en el espacio de Hilbert H1Ω;R3, entonces vn=un+wncon un∈Rig Ω; R3⊥ywn∈ Rig Ω;R3. Las relaciones anteriores implican que: kunkH1(Ω;R3)= 1 for all n ke(un)kL2(Ω;S3)→0as n→ ∞. Entonces existe una función creciente σ:N→Ntal que la subsucesión uσ(n)es de Cauchy en H1Ω;R3. Este espacio por ser completo, existe u∈H1Ω;R3 tal que: uσ(n)→uen H1Ω;R3y el límite satisface que: e(u) = l´ımn→∞ euσ(n)= 0. Así que u∈Rig Ω; R3. Por otra parte, u∈Rig Ω;R3⊥dado que Rig Ω;R3⊥es cerrado y un∈Rig Ω;R3⊥. Entonces u= 0, lo que contradice uσ(n)H1(Ω;R3)= 1,∀n∈N. 4.2.2. Existencia, Estabilidad y Regularidad de la Solución Débil Se define una solución débil a la ecuación de elasticidad como la solución a las ecuaciones variacionales: ZΩ σ:∇vdx =ZΩ f·vdx +ZΓ1 g·vdx (4.23) ∀v, tal que satisfazca v: Ω →R3yv= 0 en Γ0, donde: σ=λ(tre(u))I+ 2µe(u)ye(u) = 1 2∇uT+∇u. Nótese que dado que σes simétrica, la integral de la izquierda puede ser escrita como σ:∇v=σ:e(v) donde e(v) := 1 2∇vT+∇v. La existencia y unicidad de solución a este problema variacional puede probarse con el Teorema de Lax-Milgram (Teorema 2.1). Distinguimos dos casos dependiendo de si la medida de Γ0es positiva o nula. Problema variacional elástico PEvar: PEvar Hallar u∈Vtal que RΩσ:∇vdx =RΩf·vdx +RΓ1g·vdx , ∀v∈V. CAPÍTULO 4. ECUACIÓN DE ELASTICIDAD 47 4.4.1. Caso elástico en 2D Como es de esperar, el código utilizado para las vigas de distintos materiales es el mismo con la pequeña variación de las constantes elásticas. En el código, comenzamos definiendo los parámetros necesarios para tratar un problema lineal elástico. Todos los declaramos reales, por lo que incluimos la orden -real- previo a la definición de los parámetros. En concreto, necesitamos: 1. rho, se refiere a la densidad del material. Acero: 7700 kg m2. Cobre: 8960 kg m2. Hormigón: 2300 kg m2. 2. µyλ, se obtienen a partir de las constante de Poisson (ν) y del módulo de Young (E), a partir de las expresiones: ν=λ 2(λ+µ)yE=µ(3λ+ 2µ) λ+µ. En concreto, los valores que hemos usado para cada material son: Acero: µ= 11,05e10 N m2yλ= 14,06e10 N m2. Cobre: µ= 15,110e10 N m2yλ= 32,110e10 N m2. Hormigón: µ= 1,125e10 N m2yλ= 0,75e10 N m2. 3. gravity, fuerza de la gravedad con valor −9,81 m s2que es la que actúa sobre el cuerpo. 4. twomul, definido como twomul=2µ+λnos servirá a la hora de escribir el problema variacional. 5. nymno describen más que el número de elementos en los que vamos a hacer el mallado regular en xeyrespectivamente. La idea es fijar un borde y estudiar el desplazamiento en el resto de la viga. Lo anterior se corresponde con condiciones de Dirichlet homogéneas en la parte fija y de tensión normal nula en el resto. En concreto, a que ∂ΩDsea la parte fija, y∂ΩNel resto (en el problema PN). Es por esto por lo que nos basta etiquetar tres de los bordes iguales (1) y uno diferente (2). Lo realizamos mediante el comando que aparece en la línea 18 del código de cualquiera de los tres materiales. A continuación, definimos la malla Th, que es una discretización en el número de elementos nym, indicados previamente. El dominio Ωocupado por la viga CAPÍTULO 4. ECUACIÓN DE ELASTICIDAD 48 Figura 4.2: Viga 2D de acero con un borde fijo y sometida a fuerza gravitatoria es un rectángulo definido en [0,10] ×[0,1] con distancias medidas en metros. Por otra parte, realizamos una discretización por EF tipo P1del espacio test. Definimos ahora el problema, escribiendo las integrales del problema variacional, y estableciendo que la solución para el desplazamiento uha de ser tal que la fuerza elástica determinada los los coeficientes antes definidos ha de compensar la fuerza correspondiente a la gravitacional que actúa sobre el cuerpo. Finalmente, realizamos una gráfica del espacio mallado inicial y final con el desplazamiento de la solución. La figura que se obtiene es el que aparece en la Figura 4.2. Los desplazamientos son similares en forma para los tres materiales considerados, pero no en magnitud. Es por esto por lo que sólamente incluimos el desplazamiento de uno de los materiales. Es interesante realizar una comparativa de los desplazamientos reales y máximos sufridos por las vigas de cada material. Así, obtenemos que el desplazamiento máximo se obtiene para el hormigón, constituyendo por tanto, el material más elástico. El valor que se obtiene del desplazamiento para la componente x en el caso del hormigón es de 0,0106751 m. El segundo más elástico sería el acero, con un desplazamiento en la componente x de 0,00324154 m. Por último, el cobre constituye el caso menos elástico con un desplazamiento longitudinal de 0,00247876 m. Estos resultados están totalmente de acorde con el material utilizado para las vigas de las construcciones, que son principalmente de hormigón, el encofrado y de acero, la estructura. 4.4.2. Caso elástico en 3D El interés de incluir el caso 3D del problema anterior (con las mismas condiciones al caso del acero) no es otro que ver cómo se modifica el código ya que es más complejo crear la estructura y construir el mallado que en el caso anterior. Para CAPÍTULO 4. ECUACIÓN DE ELASTICIDAD 49 Figura 4.3: Viga en 3D ello cargamos el código del cubo, obtenido de [18] , definimos sus dimensiones y coordenadas y construimos la malla regular sobre la viga que se forma. Además, hay que cambiar las funciones añadiendo una más correspondiente a la tercera dimensión. El resultado que se obtiene al ejecutarlo es el mostrado en la Figura 4.3. Capítulo 5 Ecuación de Difusión- Convección-Reacción La gran complejidad de ecuaciones de transporte, entre ellas las ecuaciones de Euler (fluidos perfectos) y de Navier-Stokes (viscosos incompresibles), ha provocado que desde principios de siglo XX el tratamiento de este tipo de ecuaciones haya sido un tema de permanente estudio en el Análisis Matemático y en el Análisis Numérico debido a la incapacidad de los métodos analíticos de obtener una solución cerrada. Desde el punto de vista arquitectónico, la ecuación de difusión-convección constituye una aproximación formal a la dinámica de los flujos (calor,aire...) en el interior/exterior de los espacios. Es por esto por lo que en este capítulo se modelará dicha ecuación, y se aplicarán los métodos de Petrov-Galerkin y Galerkin ya estudiados para prepararla para la búsqueda de una solución numérica. 5.1. Introducción al Problema En esta sección trataremos a la ecuación del calor desde el punto de vista físico. Realizaremos tanto la modelización del problema partiendo de leyes naturales como el tratamiento físico de la conducción del calor en sólidos cristalinos. 5.1.1. Modelización del Problema A continuación, realizaremos una modelización de la ecuación del calor huyendo en todo momento de la formalidad que a las Matemáticas caracteriza, ya que 50 CAPÍTULO 5. ECUACIÓN DE DIFUSIÓN-CONVECCIÓN-REACCIÓN 51 no anhelamos más que un acercamiento físico al comportamiento del calor en un medio. El modelado de la ecuación del calor se basa en una aproximación partiendo de una Ley de Conservación. Permítenos el lector que cometamos un abuso del lenguaje al escribir calor contenido en un sólido para referirnos al proceso de transferencia de calor en un sólido, siendo la energía térmica la magnitud que realmente se posee. En todo momento obviaremos la radiación como proceso de transferencia calorífica, considerando únicamente difusión y convección. Sea Ω⊂Rdun dominio para algún d∈N. Generalmente se tomará d= 4, que significará que estamos considerando un volumen y la evolución temporal del calor en él. Así construiremos nuestro modelo, suponiendo que Ω = V×S, V⊂Rd−1yS⊂R. Partimos de la hipótesis de que el calor contenido en un sólido o fluido por unidad de volumen puede modelarse por: H: Ω →R (x, t)7→ H(x, t) = σ(x, t)c(x, t)T(x, t) y donde: x∈Vyt∈Sson las variables espacial y temporal respectivamente. σ(x, t)es la densidad volumétrica del sistema donde se transporta. c(x, t)es el calor específico del material. T(x, t)es la temperatura. Observación 5.1. Es trivial ver que el calor total, N, contenido en el sistema viene dado por: N:S→R t7→ N(t) = ZV H(x, t)dV. (5.1) Definición 5.1. (Flujo) Se denomina como flujo −→ qa la cantidad de una determinada magnitud que atraviesa una superficie por unidad de tiempo, y tiene como dirección la de la normal a la superficie. Proposición 5.1. El flujo calorífico en un sólido viene dado por la expresión: −−−−→ q(x, T) = −κ(x, t)∇T= −−−→ ^ q(x, t), donde κ(x, T)>0se refiere a la conductividad térmica del material. Observación 5.2. La única forma de transferencia de calor en un sólido es la difusión, ya que basándonos en la Teoría Cinético Molecular, las partículas que conforman al sólido tienen como único grado de libertad la oscilación entorno a una posición de equilibrio. CAPÍTULO 5. ECUACIÓN DE DIFUSIÓN-CONVECCIÓN-REACCIÓN 52 Figura 5.1: Normales a la superficie del cono Observación 5.3. Dado que −−−−→ q(x, T) = −−−→ ^ q(x, t), al ser T=T(x, t), para aliviar notación voy a escribir a partir de ahora q(x, t). Por el mismo motivo, evitaremos escribir el argumento de las funciones en ciertos momentos que puedan sobreentenderse. Proposición 5.2. El flujo calorífico de un líquido viene dado por la expresión: −→ q(x, t) = −→ q(x, t)dif +−→ q(x, t)conv =−κ(x, t)∇T(x, t) + H(x, t)−→ β(x, t) donde −→ β(x, t)se refiere a la velocidad de convección del fluido. Observación 5.4. De acuerdo a la Teoría Cinético Molecular, las partículas que conforman al líquido están unidas entre ellas por fuerzas de cohesión más débiles que en el caso de los sólidos. Es por esto por lo que las moléculas de los líquidos, aparte de vibrar, pueden desplazarse relativamente. Tenemos, por tanto, dos formas de transferencia de calor: la difusión y la convección. Proposición 5.3. Sea −→ qel flujo de la magnitud H,−→ nel vector normal a∂V yf∈C0(Ω) las fuentes o sumideros de campo por unidad de volumen. Se tiene entonces que la variación temporal de la cantidad total de calor N(t) puede expresarse como: N0(t) = ZΓ −−→ q·−→ n dΓ + ZV f dV, (5.2) CAPÍTULO 5. ECUACIÓN DE DIFUSIÓN-CONVECCIÓN-REACCIÓN 53 donde Γse refiere a la frontera de Vque ha de ser de tipo C1. Además, se exige la regularidad de H∈C1(Ω) para que lo anterior tenga sentido. Observación 5.5. El producto −→ q·−→ nsale hacia afuera del volumen, por lo que ha de contribuir negativamente a la cantidad de calor contenido. Esto justifica la elección del signo negativo. Este integrando modela la cantidad de calor que atraviesa la superficie que encierra al volumen. Proposición 5.4. Supongamos que H, q yfson funciones con regularidad C1(Ω)∩C0(Ω) y que Γes de tipo C1. Entonces, la variación temporal del calor total contenido en el volumen Vviene dada por: N0(t) = −ZV (∇−→ q−f)dV. (5.3) Demostración. La prueba es consecuencia directa de aplicar a (5.2) la Fórmula de Green. N0(t) = ZΓ −−→ q−→ n dΓ+ZV f dV =−ZV ∇−→ q dV +ZV f dV =−ZV (∇−→ q−f)dV. Por otra parte, de derivar (5.1) se obtiene que: N0(t) = d dt ZV H(x, t)dV =ZV ∂H(x, t) ∂t dV, que igualando con (5.3) se llega a: ZV (∂H(x, t) ∂t +∇−→ q(x, t)−f(x, t)) dV = 0. Dado que ∂H(x,t) ∂t +∇−→ q(x, t)−f(x, t)∈L1(Ω), se tiene que ∂H(x,t) ∂t +∇−→ q(x, t)− f(x, t)=0e.c.t. Ω. Esto ecuación de continuidad constituye una Ley General de Conservación:∂H(x, t) ∂t +∇−→ q(x, t) = f(x, t). El proceso de difusión de calor en un sólido viene entonces dado por: ∂ ∂t(σ(x, t)c(x, t)T(x, t)) − ∇(κ(x, t)∇T(x, t)) = f(x, t).(5.4) CAPÍTULO 5. ECUACIÓN DE DIFUSIÓN-CONVECCIÓN-REACCIÓN 54 Si tratamos con materiales homogéneos en los que sus propiedades de κ(x, t), c(x, t) yσ(x, t)son constantes con respecto al punto material, llegamos a que la ecuación que modela el proceso de difusión se puede expresar como: ∂T(x, t) ∂t −eκ∇2T(x, t) = g(x, t),(5.5) donde eκ=κ σc yg(x, t) = 1 σc f(x, t). Para obtener la ecuación equivalente para el caso de los fluidos, basta tener en cuenta la forma del flujo dada en la Proposición 5.2. En este caso, la ecuación que modela la convección difusión en un fluido homogéneo viene dada por: ∂T(x, t) ∂t +∇ · (T(x, t)−−−−→ β(x, t)) −∇·(κ(x, t)∇T(x, t)) = f(x, t).(5.6) Observación 5.6. Tanto la ecuación de difusión en sólidos como la de convección-difusión en fluidos son ecuaciones elípticas, para las que se puede hacer un estudio temporal o estacionario anulando el término correspondiente a la derivada temporal de la temperatura. 5.1.2. Justificación física de transmisión del calor en sólidos Este apartado es para comentar brevemente cómo se estudia en Física un determinado tipo de los sólido en cuanto a su capacidad de transmitir el calor. Nos referimos a materia condensada a aquella caracterizada por el hecho de que la distancia media entre las entidades que la forma es de aproximadamente 1A . Lo que determina si un material es sólido o líquido es la energía cinética Ten relación a la energía potencial Uque mantiene cohesionados a los constituyentes de la materia. Líquidos,T≈U. Las partículas se mueven sin perder el contacto unas con otras pero sin la capacidad de mantenerse en una determinada posición. El hecho de que se observen ciertos grados de orden está relacionado con el morfismo de las partículas. Sólidos,T < U. Las partículas vibran respecto a una posición fija. La nula capacidad de desplazamiento relativo entre los entes que conforman al sólido provocan que sea más fácil encontrar estructuras con un elevado grado de orden. En cuanto a ésto, podemos clasificar los sólidos en: 1. Sólidos cristalinos. Son sólidos cuya estructura atómica se basa en la repetición de un patrón denominado: el retículo, formando las redes CAPÍTULO 5. ECUACIÓN DE DIFUSIÓN-CONVECCIÓN-REACCIÓN 55 de Bravais. Presentan, como es natural, alta simetría y orden de largo alcance. 2. Sólidos amorfos. Tienen un orden de corto alcance, con proceso de solidificación rápido y calor latente inexistente. Nos centramos particularmente en los sólidos cristalinos. Los átomos de los sólidos cristalinos a temperaturas superiores al cero absoluto (0K) están unidos a posiciones de equilibrio en torno a las cuales oscilan como consecuencia de tener cierta energía térmica. Debido a que los átomos están unidos a enlaces, las vibraciones se propagan por la red cristalina en forma de ondas elásticas. El tratamiento cuántico de las vibaraciones reticulares nos lleva a la necesidad de tener un ente corpuscular que satisfazca la dualidad ondacorpúslo para estos paquetes de ondas elásticas. Nace así el concepto de fonón, que tiene un papel esencial en la descripción de numerosas propiedades físicas, especialmente las térmicas y las eléctricas. En un sólido aislante, la temperatura se transmite únicamente a través de las vibraciones reticulares. Si el sólido es conductor, habría que considerar la contribución al desplazamiento electrónico por agitación térmica. En por esto, por lo que en los sólidos el concepto de fonón es esencial. En Física, la propagación de los fonones en su la red se modela con la teoría cinética de los gases, quedando la conductividad: κ=1 3cVvl, (5.7) donde cVse refiere al calor específico a volumen constante, va la velocidad de propagación de los fonones que es aproximadamente la velocidad del sonido en el medio y lal recorrido libre medio entre colisiones. Dado que vpuede considerarse constante, la dependencia de κcon Treside en los términos de cVyl. Para explicar la dependencia del calor específico cVcon la tempertatura Texisten numerosos modelos que funcionan mejor o peor en función del rango de temperaturas que estemos analizando. Quizás, una de las mejores aproximaciones teóricas a este comportamiento físico sea el Modelo de Debye ver [21]. Este modelo muestra una dependencia de cVcon altas y bajas temperaturas, muy parecido que en lo que en el laboratorio de forma experimental se obtiene. La premisa básica que asume es que la distribución de frecuencias en el sólido sigue la distribución de un sólido contínuo homogéneo y, con un tratamiento propio de física estadística, este modelo establece que el calor específico en un sólido cristalino cumple la siguiente ecuación: CV≈12π4 5NrkBT θD3 .(5.8) donde Nse refiere al número de celdas unidad, ral número de átomos por celda unidad, kBes la constante de Boltzman, Tes la temperatura y θDes la temperatura de Debye [21]. Por su parte, la dependencia de lcon Tpuede explicarse cuando consideramos CAPÍTULO 5. ECUACIÓN DE DIFUSIÓN-CONVECCIÓN-REACCIÓN 56 los térmicos anarmónicos del potencial de interacción entre los constituyentes del sólido. En concreto, hay que considerar procesos de interacción fonónica de tipo U, que (como idea general) son procesos en los que la energía de los fonones incidentes es alta y, como resultado cambian su sentido de propagación y se crea una resistencia que afecta a la conductividad térmica del sólido. Este modelo demuestra que a altas temperaturas, se tiene: l∝1 T,(5.9) y a bajas temperaturas l∝expθD T.(5.10) Cuando los sólidos cristalinos son conductores eléctricos, hay que tener en Figura 5.2: Teoría de Debye cuenta la contribución a la conductividad térmica del término relativo al calor transportado por los electrones de conducción. Para obtener la expresión de κ se utiliza el modelo de Sommerfeld y la distribución fermiónica que siguen los electrones por tener espín semientero. Se llega a: κ=π2nk2 BτFT 3m,(5.11) donde los términos que no aparecen en la κde vibraciones reticulares son τF que es la tempertaura de Fermi, nla densidad de electrones de conducción y m su masa. Además, en los metales, de acuerdo a la Ley de Wiedemann-Franz el cociente entre conductividad eléctrica y térmica a una temperatura absoluta fijada es aproximadamente constante para todos ellos, implicando en particular que ambas conductividades son proporcionales para cada metal. De acuerdo a la Teoría de Bandas, los metales son conductores. En los conductores, la banda de conducción y de valencia solapan, de forma que los electrones del último orbital (o los últimos dependiendo del solapamiento), quedan en la capa de conducción, significando esto que la energía potencial debida al núcleo atómico se apantalla con la energía cinética debido a la movilidad que tienen en esta capa. Una alta movilidad electrónica indica una alta capacidad de conducción de corriente. Para el lector interesado, consultar [21]. CAPÍTULO 5. ECUACIÓN DE DIFUSIÓN-CONVECCIÓN-REACCIÓN 63 Teorema 5.2. (Estimación H1)Sea Ωun poliedro en Rdy{Th}h>0una familia regular de mallas conformes de Ω. Sea Vhtal y como se define en (5.21). Se tiene entonces que ku−uhk1,Ω= 0. Además, si u∈Hs(Ω) con d 2< s ≤k+1, ∃c > 0tal que: ∀h, ku−uhk1,Ω≤chs−1|u|s,Ω.(5.24) Demostración. Para el lector interesado, consultar [16]. Teorema 5.3. (Estimación L2)Dadas las hipótesis del Teorema 5.2, asúmase que V=H1 0(Ω),Vh=Lk∩H1 0(Ω), y que el problema Paprox tiene propiedades de suavidad. Entonces, ∃c > 0tal que ∀h, ku−uhk0,Ω≤ch|u−uh|1,Ω.(5.25) Demostración. Para el lector interesado, consultar [16]. Para el caso del problema de Dirichlet no homogéneo, que recordemos tenía la forma: CDR.NH    Hallar u∈H1(Ω) tal que a(u, v) = RΩfv, ∀v∈H1 0(Ω) γ0(u) = g, en H1 2(∂Ω). con γ0el operador traza que aparece definido en el Teorema A.11. Supongamos que el problema está bien planteado en H1 0(Ω)×H1 0(Ω) gracias a las condiciones ofrecidas en el Teorema de BNB (Teorema 2.5). Tomamos a Vh=Lk c,h yN= dim Vhcon {ϕ1, ..., ϕN}como base y {a1, ..., aN}los nodos asociados. Por otra parte recordando que el interpolante de Lagrange de una función continua uen Ωestá definido por: Ihu= N X i=1 u(ai)ϕi,(5.26) tenemos que si ges una función continua en ∂Ωentonces, ha de ser: Ihg=X ai∈∂Ω g(ai)γ0(ϕi).(5.27) Veamos que la traza del interpolante de una función suficientemente regular coincide con el interpolante de su traza. Dado que {ϕ1, ..., ϕN}es una base nodal, (ai/∈∂Ω) =⇒(γ0(ϕi) = 0). CAPÍTULO 5. ECUACIÓN DE DIFUSIÓN-CONVECCIÓN-REACCIÓN 64 Por lo que si u∈C0(¯ Ω) ∩H1(Ω), γ0(Ihu) = γ0 N X i=1 u(ai)ϕi!= N X i=1 u(ai)γ0(ϕi) =X ai∈∂Ω u(ai)γ0(ϕi) = I∂ h(γ0(u)) (5.28) yγ0◦ Ih=I∂ h◦γ0. Considera el problema aproximado de CDR.NH: CDR.NHaprox    Hallar uh∈Vh(Ω) tal que a(uh, vh) = RΩfvh,∀vh∈Vh0 γ0(uh) = I∂ hg, en ∂Ω. Con Vh0={vh∈Vh, γ0(vh)=0} ⊂ H1 0(Ω) y asumiendo que la forma bilineal asatisface la condición BNB1hde la Proposición 3.1 en Vh0×Vh0. Proposición 5.8. Sea guna función lo suficientemente regular como para tener una función de levantamiento en C0(¯ Ω)∩H1(Ω). En tales condiciones, el teorema CDR.NHaprox está bien planteado. Demostración. Sea ugel levantamiento de gen C0(¯ Ω) ∩H1(Ω). Claramente, γ0(Ihug) = I∂ h(γ0(ug)) = I∂ h(g) = γ0(uh). Así que, θh=uh− Ihugnos lleva a θh∈Vh0ya(θh, vh) = RΩfvh−a(Ihug, vh) para todo vh∈Vh0. Dado que asatisface BNB1hen Vh0×Vh0, el problema CDR.NH (aprox) está bien planteado. Una vez probado el buen planteamiento del problema CDR.NHaprox, nuestro objetivo es estimar el error u−uh. Para ello, enunciamos brevemente y sin deternos demasiado los siguientes resultados extraídos de [16]. Lema 5.1. Sea guna función lo suficientemente regular como para tener una función de levantamiento en C0(¯ Ω) ∩H1(Ω), sea u una solución lo suficientemente regular para que su interpolador de Lagrange esté bien definido. Sea, además, kak:= kakH1(Ω),H1(Ω). Se tiene entonces que: ku−uhk1,Ω≤1 + kak ahku− Ihuk1,Ω.(5.29) Lema 5.2. Sean las condiciones del Lema 5.1, asúmase además que: El problema CDR.NHaprox tiene propiedades de regularidad. CAPÍTULO 5. ECUACIÓN DE DIFUSIÓN-CONVECCIÓN-REACCIÓN 65 La forma bilineal asatisface la siguiente propiedad de continuidad: ∃ctal que, ∀v∈H1(Ω) y∀w∈H2(Ω), |a(v, w)| ≤ c(kvk0,Ω+kγ0(v)k0,∂Ω)kwk2,Ω. Existe una constante de interpolación c > 0tal que: ∀h, ∀θ∈H2(Ω),kθ− Ihθk1,Ω≤chkθk2,Ω.(5.30) Existe entonces una constante ctal que: ∀h, ku−uhk0,Ω≤chkIhu−uk1,Ω+kIhu−uk0,Ω+kIhg−gk0,∂Ω. (5.31) Corolario 5.1. Sea Ωun poliedro, {Th}h>0una familia regular de mallados conforme de Ω, y Vhuna aproximación conforme de H1(Ω) mediante EF de Lagrange de grado k≥1. Satisfechas las hipótesis del Lema 5.2 yu∈HK+1(Ω). Entonces, se tiene que: ∀h, ku−uhk0,Ω+hku−uhk1,Ω≤chk+1kukk+1,Ω.(5.32) 5.3.2. La Aproximación no Conforme Crouzeix-Raviart El objetivo ahora es presentar una aproximación no conforme para la ecuación de convección-difusión. Vamos a realizar el análisis para el caso en el que únicamente hay difusión, lo que se traduce en trabajar con una forma bilineal que es un laplaciano. Es una simplificación importante en cuanto a lo que a facilitar al tratamiento del problema se refiere. Pero es real. Tenga el lector en cuenta que el mecanismo de transporte de calor en un sólido puede modelarse esencialmente por un laplaciano, ya que no existe transporte por convección en este caso. Bajo esta suposición, el operador Ldado en (5.12) queda como: Lu=σ∆u. En la expresión anterior también hemos supuesto que el coeficiente de conducción térmica σes constante. Vamos a hacer la aproximación usando EF de Crouzeix-Raviart. Para ello, consideremos Ωun poliedro en Rdyf∈L2(Ω). Bajo estas condiciones, en concreto por ser todo poliedro un espacio convexo, upuede verse como función de H2(Ω) [16], siendo usolución al problema variacional siguiente: Hallar u∈Vtal que a(u, v) = f(v),∀v∈V. CAPÍTULO 5. ECUACIÓN DE DIFUSIÓN-CONVECCIÓN-REACCIÓN 66 Consideremos {Th}h>0una familia de mallados de Ωregular afín y conforme. Sea P1 pt,h el espacio de EF de Crouzeix-Raviart definido por: P1 pt,h =vh∈L1(Ωh); ∀K∈ Th, vh|K∈P1 ∀F∈ Fi h,ZF [vh] = 0}.(5.33) Consideremos, además, el siguiente espacio: P1 pt,h,0=vh∈P1 pt,h;∀F∈ F∂ hZF , vh= 0, donde F∂ hdenota el conjunto de caras del mallado que están sobre el contorno de Ω. Es trivial que dim P1 pt,h,0=Ni f, el número de caras internas (superficies en dos dimensiones) de la triangulación. Dado que las funciones de P1 pt,h,0pueden ser discontinuas, RΩ∇u· ∇vhay que aplicarla a cada uno de los elementos por separado. Veámoslo: Hallar uh∈P1 pt,h,0tal que ah(uh, vh) = f(vh),∀vh∈P1 pt,h,0 , donde ahes la forma aproximada de a, ah(uh, vh) = X K∈ThZK ∇uh· ∇vhyf(vh) = ZΩ fvh.(5.34) Por otra parte, V(h) = P1 pt,h,0+H1 0(Ω) y para vh∈V(h)se define la seminorma en H1: |vh|h,1,Ω= X K∈Th k∇vhk2 0,K!1 2 . Equipamos al espacio V(h)con la norma k · kV(h)=k · k0,Ω+| · |h,1,Ω. Nuestro objetivo es investigar la convergencia de solución al problema aproximado en norma k · kV(h). Para ello, tenemos que mostrar la estabilidad, continuidad, consistencia y propiedades de proximidad vistas en el Capítulo 3. Para poder hablar de estabilidad, primero hay que probar la coercitividad de ahen P1 pt,h,0 que no es una propiedad trivial ya que P1 pt,h,06⊂ H1 0(Ω). Todas las pruebas de los resultados que aparecen a continuación no se exponen en este texto debido a complejidad y extensión. Para el lector interesado, consultar [16]. Tenemos que extender la desigualdad de Poincaré a esta norma. Veámoslo en el siguiente lema. Lema 5.3. Desigualdad de Poincaré extendida Existe c > 0dependien- do únicamente de Ωtal que, ∀h≤1, ∀u∈V(h), ckuk0,Ω≤ |u|h,1,Ω.(5.35) CAPÍTULO 5. ECUACIÓN DE DIFUSIÓN-CONVECCIÓN-REACCIÓN 67 Lema 5.4. Sea {Th}h>0una familia regular de mallas conformes y afines y m≥1un entero fijo. Para K∈ Th,ψ∈H1(K)my una cara F∈∂K, sea ψ=1 meas(F)RFψ. Entonces, existe c > 0tal que: ∀h, ∀K∈ Th,∀F∈∂K, ∀ψ∈H1(K)m,kψ−ψk0,F ≤ch1 2 K|ψ|1,K.(5.36) Demostración. Sea K∈ Thyψ∈H1(K)my considera una cara F∈∂K de la malla. Sea, además, ˆ Kser el simplex de referencia y TK:b K→Kla transformación afin con el Jacobiano JK. Haciendo b F=T−1 K(F), es trivial ver que debido al Teorema de la Traza(INCLUIR REF) se tiene que: kψ−ψk0,F ≤measF meas ˆ F1 2 kb ψ−b ψk0,b F≤cmeasF meas ˆ F1 2 kb ψ−b ψk1,b K. Del Teorema de Deny-Lions (ver [16]) se obtiene kb ψ−b ψk1,b K≤c|b ψ|1,b K. Volviendo al elemento Ky usando la regularidad de la malla finalizamos la prueba kψ−ψk0,F ≤cmeasF measF1 2J−1 Kd meas b K masK!1 2 |ψ|1,K ≤ch d−1 2 KhKh−d 2 K|ψ|1,K ≤ch1 2 K|ψ|1,K. Corolario 5.2. (Estabilidad) La forma bilineal ahes coercitiva en P1 pt,h,0 por simple aplicación de la Desigualdad de Poincaré extendida. Lema 5.5. (Consistencia asintótica) Sea usolución al problema de Dirichlet homogéneo con f∈L2(Ω) y asúmase además que u∈H2(Ω). Existe entonces c > 0tal que : ∀h, ∀wh∈P1 pt,h,0,|f(wh)−ah(u, wh)| kwhkV(h) ≤ch|u|2,Ω.(5.37) Demostración. Sea wh∈P1 pt,h,0. Dado que f=−∆u, ah(u, wh)−f(wh) = X K∈ThZK (∇u· ∇wh−fwh) = X K∈ThX F∈∂K ZF ∇u·nKwh. CAPÍTULO 5. ECUACIÓN DE DIFUSIÓN-CONVECCIÓN-REACCIÓN 68 Dado que cada cara Fdel elemento Klocalizada en el interior de Ωaparece doblemente en la suma, podemos restarle el valor medio de when la cara, denotémoslo por wh. Si Festá en δΩes claro ver que wh= 0. Así que: ah(u, wh)−f(wh) = X K∈ThX F∈∂K ZF ∇u·nK(wh−wh). Podemos además restar de ∇usu valor medio en la cara F,∇ullevándonos a: ah(u, wh)−f(wh) = X K∈ThX F∈∂K ZF (∇u− ∇u)·nK(wh−wh). La desigualdad de Cauchy-Schwarz implica: |ah(u, wh)−f(wh)| ≤ X K∈ThX F∈∂K k∇u− ∇uk0,F kwh−whk0,F y aplicando el Lema 5.4, llegamos a: |ah(u, wh)−f(wh)| ≤ PK∈Thch1 2 K|u|2,Kh1 2 K|wh|1,K ≤ch PK∈Th|u|2 2,K ·PK∈Th|wh|2 1,K1 2≤ch|u|2,ΩkwhkV(h). Lema 5.6. (Proximidad) Existe una constante c > 0tal que: ∀h, ∀u∈H2(Ω) ∩H1 0(Ω),´ınf vh∈P1 pt,h,0 ku−vhkV(h)≤ch|u|2,Ω.(5.38) Demostración. Ver [16]. Teorema 5.4. (Convergencia) Bajo las condiciones del Lema anterior, existe c > 0tal que: ∀h, ku−uhkV(h)≤ch|u|2,Ω.(5.39) Teorema 5.5. (Estimación de error L2)Bajo las condiciones del teorema anterior, asúmase que el laplaciano tiene propiedades de regularidad en Ω. Existe entonces c > 0tal que: ∀h, ku−uhk0,Ω≤ch |u−uh|h,1,Ω.(5.40) CAPÍTULO 5. ECUACIÓN DE DIFUSIÓN-CONVECCIÓN-REACCIÓN 69 Luego, con los resultados anteriores hemos conseguido ver mediante la Desigualdad de Poincaré extendida que la forma bilineal es coercitiva, resultado que implica estabilidad. Además, si f∈L2(Ω yu∈H2(Ω) nos aseguramos la consistencia asintótica. Mediante la proximidad vemos que el espacio Vhaproxima bien al espacio V. Finalmente, recordemos que tener consistencia y estabilidad nos asegura tener convergencia. 5.4. Implementaciones en FreeFem++ A continuación, vamos a realizar diferentes implementaciones de la ecuación de transporte. Por una parte, trataremos el caso de difusión pura, característico principalmente en el caso de los sólidos. Estudiaremos el caso difusivo para materiales con diferentes constantes térmicas, comparando los resultados que se obtienen. Por otra parte, estudiaremos un caso algo más real del comportamien- to del calor en las paredes de una habitación con un radiador. En este caso, las paredes están compuestas con distintos materiales. Como comentario, añadiremos que no se ha pretendido cubrir todos los casos antes expuestos en teoría, sólamente trazar los primeros pasos de la implementación de los casos más sencillos para afianzar el contenido estudiado y, al mismo tiempo, dar una visión más cotidiana de estas ecuaciones. 5.4.1. Primer Caso: Comparativa Materiales El primer problema que queremos resolver es el siguiente:                      ∂u ∂t −∇·(κ∇u)en (0, T )×Ω u(·, t = 0) = θen Ω, ∂u ∂n =−1sobre (0, T)×Γ1 u=θsobre [0, T)×Γ2, donde Γ1∪Γ2=∂Ωyθes la temperatura inicial a la que se encuentra el sistema Ω. Si el material donde se produce la difusión es aire, el anterior podría ser la simplificación del caso de una habitación Ωque se ha mantenido a una temperatura θy quitando la fuente que ha conseguido que la temperatura sea constante, se deja otra sobre Γ2a temperatura θy quiere estudiarse como evoluciona la temperatura considerando las pérdidas que hay a través de la pared Γ1. Suponiendo, además, que no hay transporte por convección o que es despreciable frente a la difusión. Esto podría darse en una habitación sin corrientes de aire, por ejemplo. Por los motivos arquitectónicos que conciernen a este trabajo, llamaremos CAPÍTULO 5. ECUACIÓN DE DIFUSIÓN-CONVECCIÓN-REACCIÓN 70 al sistema Ωhabitación y a Γ2radiador, independientemente de si el material donde se produce la difusión es aire u otro diferente. Mediante el Método de Euler podemos realizar la discretización temporal: un+1 −un ∆tn−∇·(κ∇un+1)=0.(5.41) Por su parte, para la discretización tenemos en cuenta que la frontera del recinto es, como se ha dicho previamente al composición de dos fronteras diferentes: Γ1frontera del recinto donde se pierde temperatura. Γ2frontera interior a una temperatura constante de θ. Por otra parte, vamos a realizar una discretización del espacio de soluciones H1(Ω) con traza igual a θsobre Γ2. Es decir, para cualquier v∈H1 2(Ω), definimos el espacio: Vg={v∈H1(Ω), v =gsobre Γ2}.(5.42) Obsérvese que en nuestro caso tenemos Vθ. Por otra parte, para el espacio test, consideremos V0. Transformamos nuestro problema anterior a su formulación variacional. Para ello, multiplicamos por v∈V0, aplicamos la fórmula de Green y las condiciones de contorno. Mediante el Teorema de Lax Milgram, se tiene que para cada un+1 existe y es única. El problema abstracto queda como:    Hallar u∈Vθtal que 1 ∆tnZΩ un+1vdx +ZΩ κ∇un+1 · ∇vdx +ZΓ1 κvdσ =1 ∆tnZΩ unvdx, ∀v∈V0. Discretizamos ahora los espacios mediante EF P1de Lagrange, quedándonos ahora: Vh={vh∈H1(Ωh), vh∈C0(Ωh), vh |K∈P1(K),∀K∈ Th}(5.43) Vh g={vh∈Vh, vh=ghsobre Γh 2}para gh∈H1 2(Γh 2).(5.44) Finalmente, el problema discreto es:    Hallar uh∈Vh θtal que 1 ∆tnZΩ uh,n+1vhdx +ZΩ κ∇uh,n+1 · ∇vhdx +ZΓ1 κvhdσ =1 ∂tnZΩ uh,nvdx, ∀vh∈Vh 0. Obsérvese que el problema discreto también tiene solución y es única por aplicación del Teorema de Lax Milgram. Hemos realizado el estudio para tres materiales bien diferentes en cuanto a su comportamiento difusivo: CAPÍTULO 5. ECUACIÓN DE DIFUSIÓN-CONVECCIÓN-REACCIÓN 71 Madera de Roble, con κ= 0,17 W m·oC. Plata, con κ= 419 W m·oC. Aire, con κ= 0,024 W m·oC. Figura 5.3: Mallado Hemos elegido un recinto con forma elíptica y un radiador interno a temperatura θ= 25oC, que es a su vez la temperatura inicial del sistema. El mallado que realizamos en FreeFem++ sobre la habitación es el que indica la Figura 5.2. Obsérvese que para el entorno del radiador hemos elegido una densidad de EF mayor, lo que es totalmente lógico ya que en su entorno es donde habrá un cambio de temperaturas mayor. En el código, nymse refieren justamente al número de EF y la malla se construye con el comando típico buildmesh. A continuación se muestran los resultados que se obtienen para cada material con un salto de tiempo dt = 1sdesde t= 1s. CAPÍTULO 5. ECUACIÓN DE DIFUSIÓN-CONVECCIÓN-REACCIÓN 72 (a) Madera 1 (b) Madera 2 (c) Madera 3 (d) Madera 4 Figura 5.4: Difusión en Madera Madera 1 2 3 4 T. mín oC24.5010 24.2080 23.9655 23.7493 Apéndice A Elementos del Análisis Funcional Con el objetivo de poder realizar un análisis profundo y pulcro desde el punto de vista matemático de la aproximación de las ecuaciones de convección-difusión y elasticidad mediante el método de los elementos finitos, en este capítulo se expondrán los aspectos necesarios de la teoría de espacios de funciones para los enunciados, demostraciones y resultados generales de los que constará la memoria que aqui se presenta. Los resultados correspondientes a esta capítulo se han tomado principalmente de [1], [2], [5], [3] y [11]. A.1. Espacios de Hilbert y de Banach Esta sección está enfocada a describir algunos conceptos y resultados sobre los espacios de Hilbert y de Banach. Los espacios de Hilbert constituyen la generalización más inmediata a espacios de dimensión infinita de los espacios euclídeos finitos. De hecho, la intuición geométrica desempeña un papel importante en muchos aspectos de su teoría, y además surgen de modo natural y frecuente en matemáticas, física e ingeniería; constituyendo herramientas esenciales en la teoría de ecuaciones en derivadas parciales, mecánica cuántica y procesamiento de señales. Por su parte, los Espacios Normados son espacios en los que aún sin producto escalar puede definirse una norma. En concreto, nos centraremos en los Espacios de Banach (espacios normados completos) por su importancia en el área del análisis numérico, donde lo que se construye son sucesiones convergentes. Definición A.1. Se dice que Hes un espacio prehilbertiano si Hes 79 APÉNDICE A. ELEMENTOS DEL ANÁLISIS FUNCIONAL 80 un espacio verctorial sobre un cuerpo Ken el que está definida una aplicación h·,·i :H×H→Kllamada producto escalar. Dicha aplicación ha de satisfacer: 1. hx, xi ≥ 0yhx, xi= 0 si y sólo si, x= 0 ∀x∈H. 2. hx, yi=hy, xi∗∀x, y ∈H, y ∗denotando el conjugado. 3. hλx, yi=λhx, yi ∀λ∈K,∀x, y ∈H. 4. hx+y, zi=hx, zi+hy, zi ∀x, y, z ∈H. Definición A.2. Si Hes un espacio prehilbertiano y x∈H, al número real no negativo phx, xise le llama norma cuadrática de xy se denota como ||x||. Diremos que la aplicación x∈H7→ ||x|| es una seminorma en caso de cumplir todas las propiedades de la norma salvo el punto 1. en la Definición 2.1.1. Por una parte, se deduce de forma inmediata que cualquier espacio prehilbertiano Hes un espacio topológico, en el que todo subconjunto abierto puede tomarse como unión de las bolas abiertas B(x, r) := {y∈H:d(x, y)< r}(x∈ H, r > 0) de la topología. Esto hace que, en particular, podamos hablar de aplicaciones H→T, con Tcualquier espacio topológico. Por otra parte, desde el punto de vista métrico, con el objetivo de poder definir el concepto de espacio de Hilbert, recordemos los siguientes conceptos relativos a la continuidad y convergencia en un espacio métrico general (X, d), donde d denota la distancia: Se dice que una sucesión xn⊂Xes de Cauchy cuando, dado  > 0, existe n0=n0()∈Ntal que d(xm, xn)<  ,∀m, n ≥n0. Toda sucesión convergente es de Cauchy (lo que no se cumple en sentido contrario a nivel general). El espacio métrico (X, d)se dice completo cuando toda sucesión de Cauchy es convergente. Sean A⊂Xyx0∈A. Se tiene que x0∈¯ Asi y solo si ∃xn⊂Atal que xn→x0.¯ Ase denomina como la clausura, cierre o adherencia de A. En particular, obtenemos que Aes cerrado si y solo si [xn⊂A, xn→x0=⇒ x0∈A]. Definición A.3. Se dice que Hun espacio prehilbertiano es, además, espacio de Hilbert si es completo. Definición A.4. Sea a:H×H→Runa forma lineal con Hespacio de Hilbert. Dicha forma lineal se dice: continua, si ∃c∈R+tal que |a(u, v)| ≤ c|u||v| ∀u, v ∈H coercitiva o elíptica, si ∃α∈R+tal que |a(v, v)| ≥ α|v|2∀v∈H APÉNDICE A. ELEMENTOS DEL ANÁLISIS FUNCIONAL 81 Definición A.5. Sea Xun espacio vectorial. Se dice que una aplicación || · || :X→Res una norma, si ∀x, y ∈Xy∀λ∈R: 1. ||x|| ≥ 0y||x|| = 0 si y solo si x= 0 2. ||λx|| =|λ|||x|| 3. ||x+y|| ≤ ||x|| +||y||. Definición A.6. Un espacio Xequipado con una norma, se dice espacio normado. Observación A.1. Como antes se vio, todo espacio prehilbertiano con la norma cuadrática es un espacio normado. Definición A.7. Un espacio normado es separable si contiene un subconjunto denso numerable. Definición A.8. Se llama espacio de Banach a un espacio normado que es completo para la distancia inducida por su norma. Definición A.9. Dados dos espacios normados V, W y un operador lineal A:V→W, se dice que el operado Aes: continuo, si dados x∈V, y{xn} ⊂ Vcon xn→x, se tiene que Axn→ Ax, con convergencia en W cerrado, si dado D⊂Vy para cada {xn} ⊂ Vcon xn→x, con x∈V tal que Axn→y∈W, se tiene que x∈DyAx =y. acotado si verifica que: sup x∈V ||Ax||W ||x||V <∞. Proposición A.1. Todo operador A:V→Wcontínuo entre espacios de Banach es cerrado. Definición A.10. El conjunto de los operadores lineales acotados de Ven Wse llama L(V, W ). Proposición A.2. El espacio L(V, W)es un espacio vectorial equipado con la norma: ||A||L(V,W )= sup x∈V ||Ax||W ||x||V . Definición A.11. Sea Vun espacio de Banach. Se denomina espacio dual de Val espacio L(V, R). Se denota por V0. Proposición A.3. Sea Vun espacio de Banach y Wespacio de Hilbert. Sus espacios duales V0yH0son respectivamente Banach y Hilbert. Proposición A.4. Dado Vespacio de Banach, existe un monomorfismo Γ : V→V00 que viene dado por: f(v) = hf, viV0,V =hΓv, fiV00,V ∀f∈V0. APÉNDICE A. ELEMENTOS DEL ANÁLISIS FUNCIONAL 82 Definición A.12. Se dice que Vespacio de Banach es reflexivo si el monomorfismo anterior Γes biyectivo. Teorema A.1. (De Representación de Riesz para espacios de Hilbert) Sea Vespacio de Hilbert y V0su espacio dual. Se tiene que los espacios V yV0son isomorfos como espacios de Hilbert mediante el isomorfismo canónico R:V→V0dado por: hu, viV=hR(u), viV0,V ∀v∈V. Además, se tiene que ||R||L(V,V 0)=||R−1||L(V0,V )= 1. Observación A.2. Todo espacio de Hilbert es, trivialmente, reflexivo. Definición A.13. Sean V, W espacios de Banach y A:V→Wun operador. Se define como operador traspuesto de Acomo AT:W0→V0, dado por la asignación: ATf, vV0,V =hf, AviW0,W . Definición A.14. Dado un espacio de Banach Vy un subconjunto C⊂V, se llama anulador de Cal subconjunto: Co={f∈V0/f(v)=0∀v∈C}. Proposición A.5. Dado un espacio de Banach Vy un subconjunto Cde V, si Co={0}entonces, Ces denso en V. Teorema A.2. (De la Inversa Acotada) Dado un espacio de Banach V , un espacio normado Wy un subespacio E⊂V, se tiene que si A:E→Wes operador cerrado y el rango de A,A(E)es de segunda categoría entonces: 1. Aes sobreyectivo, es decir, A(E) = W. 2. ∃m > 0tal que ∀w∈Wexiste algún v∈Econ Av =wy||v||V≤ m||w||W. 3. Si ∃A−1, entonces es continuo. Teorema A.3. (Hahn-Banach) Sea Vun espacio vectorial y Mun subespacio de V, si existen p:M→R≤0tal que p(v+w)≤p(v) + p(w),∀v, w ∈V p(tv) = tp(v),∀v∈V, t ∈R, yf:M→(R)lineal con f(x)≤p(x)∀x∈M, entonces ∃F:V→Rtal que: F(v) = f(v)∀v∈V y −p(−v)≤F(v)≤p(v)∀x∈M. APÉNDICE A. ELEMENTOS DEL ANÁLISIS FUNCIONAL 83 Observación A.3. Dado un espacio normado Vy un subespacio M⊂V, si f∈M0,∃F∈V0tal que: ||f||M0=||F||V0y f(v) = F(v)∀v∈M. Teorema A.4. (Banach-Steinhaus) Sean Xun espacio de Banach e Y un espacio normado. Dados Iun conjunto de índices de cardinalidad arbitraria y {Aα}α∈Iun subconjunto de L(X, Y ). Se tiene que si supα∈I||Aαx|| <∞,∀x∈ X, entonces supα∈I||Aα||L(X,Y )<∞. A.2. Espacios Lp Definición A.15. Se define Lp(µ)como el espacio de funciones f:X→C medibles tales que |f|pes integrable. Si f∈Lp(µ)se define: ||f||p=Z|f|pdµp || · ||pactúa como seminorma en Lp(µ). Teorema A.5. (Convergencia dominada) Sea fnuna sucesión en Lp(µ) que converge puntualmente a f. Si existe g∈Lp(µ)tal que ∀n,|fn| ≤ g; entonces f∈Lp(µ)yfnconverge a fen Lp(µ). Teorema A.6. (Riesz-Fischer) El espacio Lp(µ)es completo. Más precisamente si fnes una sucesión de Cauchy en Lp(µ), entonces existe una subsucesión fnkque converge en casi todo punto a una función f∈Lp(µ), además fnconverge a fen || · ||p. Definición A.16. (Espacio dual) Si Ees un espacio de Banach, el conjunto de todas las formas lineales continuas udefinidas sobre Eforman un espacio de Banach con la norma ||u|| =sup {|u(x)|:x∈E, con||x|| ≤ 1}. Este espacio normado se denota por E∗y se le llama espacio dual (topológico) de E. Definición A.17. (Exponente dual) Dado p∈N, se define como exponente dual o conjugado de p al único q∈N,q > 1, tal que 1 p+1 q= 1. Teorema A.7. Sea 1< p < ∞, y qel exponente conjugado. El dual de Lp(µ)se identifica con Lq(µ), en el sentido de que a cada forma lineal y continua u∈Lp(µ)∗se le puede asociar una función g∈Lq(µ)de forma que ∀f∈Lp(µ) u(f) = Zfgdµ. APÉNDICE A. ELEMENTOS DEL ANÁLISIS FUNCIONAL 84 Además ||u|| =||g||q. Recíprocamente, cada función g∈Lq(µ)determina mediante (2.3) una forma lineal y continua en Lp(µ). El caso p= 1 es más complicado, lo vemos a continuación. Definición A.18. Se define el espacio L∞(µ)como el conjunto de las funciones medibles y acotadas f:X→C. Es un espacio vectorial. Para f∈L∞(µ) se define la función || · ||∞como: ||f||∞=inf {C > 0 : |f(x)| ≤ C ect x ∈X} Dicha aplicación es una seminorma en L∞(µ). Teorema A.8. Si µes σ−finita el dual de L1(µ)se identifica con L∞(µ) de la misma forma que para cualquier otro caso Lp(µ),p > 1. Observación A.4. El dual de L∞(µ)no se identifica con L1(µ)más que en casos triviales, siendo fácil ver que cada función de L1(µ)sí induce una forma lineal y contúnua en L∞(µ). No obstante, a nivel general este tipo de formas lineales no son todas las que existen. A.3. Teoría de Distribuciones Laurent Schwartz (1915-2002) recibió la Medalla Fields en 1950 gracias a sus trabajos sobre la Teoría de Distribuciones, la cual nace de la necesidad de crear una forma de derivación más débil y general que la tradicional forma clásica para, entre otros, poder trabajar con EDPs y dar una base formal y coherente de los cálculos realizados por físicos como Paul Dirac en el ámbito de Mecánica Cuántica. El concepto de distribución es mucho más general que el de función, en el sentido que pueden ser infinitamente derivables incluso cuando se corresponden con funciones no derivables en el sentido clásico. Sus derivadas son, también, distribuciones. Esta es la razón por la que tiene sentido hablar de soluciones distribucionales a EDPs de cualquier orden. Además, introducir el concepto de distribuciones es esencial para poder hablar de los espacios en los que principalmente trabajaremos: los espacios de Sobolev. A.3.1. Sobre el concepto de solución de una EDP Comenzamos el texto expresando la necesidad de crear un nuevo concepto de solución a un problema de EDPs. Para ello considérese el problema de estudiar la flexión de una cuerda elástica de longitud L, sujeta a los extremos y sobre la que actúa una determinada fuerza f=f(x). La Figura 2.1 representa la situación. La ecuación diferencial que modela esta situación nace de la combinación de la Ley de Hooke de la elasticidad lineal, el principio de conservación de la cantidad APÉNDICE A. ELEMENTOS DEL ANÁLISIS FUNCIONAL 85 Figura A.1: Carga distribuida de movimiento y de la consideración de que el medio que tratamos es elástico. El problema se trata entonces de resolver:    Hallar u∈C2(0, L)∩C0[0, L] tal que −(κ(x)u(x)0)0+λ(x)u(x) = f(x),∀x∈(0, L) u(0)=u(L)=0 que es básicamente hallar la solución en el sentido clásico del problema. Considérese ahora que la fuerza en lugar de ser distribuida en la longitud de la cuerda, se aplica consituyendo una carga puntual como muestra la siguiente figura: Figura A.2: Carga puntual Ésta es una situación muy real en la que, no obstante, no puede pedírsele al desplazamiento ser C2(0, L), ya que la función no es diferenciable en x=L 2. Por claridad en la notación, no escribimos la dependencia con respecto a xde todas las funciones implicadas. Si la fuerza fproduce un desplazamiento v, entonces fv es el trabajo, que integrando en toda la cuerda tenemos que: ZL 0 fv APÉNDICE A. ELEMENTOS DEL ANÁLISIS FUNCIONAL 86 Procediendo de la misma manera, aplicando la integración por partes y considerando que los extremos sufren un desplazamiento nulo: ZL 0 −(κu0)0v+λuv =ZL 0 κu0v0+λuv Lo que reuniendo todo en una misma expresión nos lleva a que para user solución del problema, ha de satisfacer: ZL 0 κu0v0+λuv =ZL 0 fv No obstante, nos seguimos encontrando con los problemas de que u∈C1(0, L) y de que u= 0 es solución del problema ya que f(x)=0,∀x∈[0, L]yx6=L 2. Ambos resultados contradicen la experiencia física. En una primera corrección del modelo, consideremos que la función fuerza verifica: f(x) = (1 2 L 2−≤x≤L 2+ 0e.c.c. Aplicando el Teorema del Valor Medio y tomando →0, el trabajo producido por esta carga viene ahora dado por: ZL 0 fv =ZL 2+ L 2− v 2=v(ξ) = v(L 2). Luego, una nueva forma de tratar matemáticamente el desplazamiento producido por medio de una carga localizada sobre un punto x0, puede hacerse mediante la intruducción de una nueva aplicación, llámese δx0que actúa sobre una determinada clase de funciones y produce valores siguiendo la regla: δx0:v→ hδx0, vi:= v(x0). La anterior es la conocida función delta de Dirac, a la que en lo que sigue se le dará una caracterización más formal. Este ejemplo ilustra claramente el que fue el origen de la Teoría de las Distribuciones. A.3.2. Espacio de Distribuciones D0(Ω) En lo que sigue, sea Ωun subconjunto abierto, conexo y acotado de RN, cuya frontera será denotada por Γ. Definición A.19. Dada g: Ω →Runa función continua, llamamos soporte de gal conjunto: sop g ={x∈Ω : g(x)6= 0} APÉNDICE A. ELEMENTOS DEL ANÁLISIS FUNCIONAL 87 Definición A.20. Se denota por Ck(Ω) al conjunto de funciones definidas en Ωque son k∈Nveces derivables en el conjunto. Definición A.21. Se denota como Ck(¯ Ω) al espacio de restricciones a ¯ Ωde las funciones de Ck(Ω0)donde Ω0es un abierto que contiene a Ω. Definición A.22. Se denota como Ck c(Ω) para k∈N∪ {0,∞}, al conjunto de funciones de Ck(Ω) con soporte compacto en Ω. Definición A.23. Se define espacio de funciones test al espacio vectorial C∞ c(Ω) y se denota como D(Ω). Observación A.5. Las funciones test se anulan en un entorno de Γ(contenido en Ω). Figura A.3: Aspecto típico de Función test A continuación va a verse que pese a que el espacio D(Ω) no sea un espacio normado, éste puede dotarse de una topología, la la topología límite inductivo, que es compatible con la estructura de espacio vectorial del mismo. Definición A.24. Sean {φn}n≥1⊂D(Ω) una sucesión y φ∈D(Ω). Diremos que φnconverge a φen D(Ω) si existe un compacto K⊂Ωtal que: sop φn⊂K, ∀n≥1y Dαφn→Dαφuniformemente en Ω,∀α∈NN. Observación A.6. La noción de convergencia anterior es compatible con la estructura de espacio vectorial de D(Ω), es decir, si φnyψnson dos sucesiones de D(Ω) convergentes en D(Ω) aφyψrespectivamente y si α, β ∈R, entonces αφn+βψn→αφ +βψ en D(Ω) D(Ω) se inyecta de forma continua en Lp(Ω), lo que significa que D(Ω) ⊂ Lp(Ω) (i.e. inclusión de conjuntos) y que si φn→φen D(Ω), entonces φn→φen Lp(Ω) (i.e. inclusión de topologías). APÉNDICE A. ELEMENTOS DEL ANÁLISIS FUNCIONAL 88 Dado que la convergencia uniforme implica convergencia puntual, se tiene trivialmente que la función φtambién tiene soporte en K. Una vez definido el espacio vectorial D(Ω) dotado con la topología anterior, estamos en condiciones de definir su dual topológico, es decir, el espacio vectorial de las formas lineales y continuas sobre D(Ω). Dicho espacio es conocido como Espacio de las Distribuciones y se denota como D0(Ω). Acudiendo a la noción de convergencia en D(Ω), el espacio D0(Ω) puede definirse de la siguiente manera: Definición A.25. Se dice que Tes una distribución en Ωsi es una aplicación T:D(Ω) →Rlineal que cumple ser secuencialmente continua, i.e.: φn→φen D(Ω) =⇒T(φn)→T(φ)en R. Dicho conjunto, denotado como D0(Ω), constituye un espacio vectorial. Se utilizará la notación hT, φipara referirnos a T(φ). Figura A.4: Distribución de Dirac δ0como límite de funciones Gaussianas normalizadas. Ejemplo A.1. (La distribución de Dirac) Fijemos x∈Ω. Se define δx mediante: δx∈D0(Ω) : ψ∈D(Ω) → hδx, ψi=ψ(x)∈R. Es trivial que δx∈D0(Ω) pues es lineal, secuencialmente continua y no es la distribución nula pues ∃ψ∈D(Ω),(ψ= 1) tal que hδx, ψi=ψ(x)=1. Proposición A.6. Si consideramos el espacio L1 loc(Ω) de las funciones f: Ω→Rtales que RK|f|<∞para Kun compacto de Ω, es trivial ver que la aplicación: Tf:D(Ω) →R ψ→ hTf, ψi=ZΩ fψ APÉNDICE A. ELEMENTOS DEL ANÁLISIS FUNCIONAL 95 Dado que el operador traza es continuo y, por tanto, acotado, existe una constante cτ=cτ(Ω, p)tal que : kf|ΓkLp(Γ) ≤cΓkfkW1,p(Ω) para todo f∈W1,p(Ω). Además para dominios Lipschitz se tiene que: H1 0(Ω) = f∈H1(Ω) : f|Γ= 0. Por otra parte, si g∈L2(Γ) se tiene que: kgkL2(Γ) ≤cΓkukH1∀u:u|Γ=g. Por lo que, kgkL2(Γ) ≤cΓm´ın u|Γ=gkukH1=cΓkgkH1/2(Γ). Proposición A.11. (Desigualdad de Poincaré) Sea u∈W1,p(Ω), se tiene que ∃CΩ∈R+: kukLp(Ω) ≤CΩk∇ukLp(Ω).(A.1) Observación A.13. En virtud de la desigualdad de Poincaré, k · kH1 0(Ω) es una norma sobre H1 0(Ω), equivalente a k·kH1(Ω). A.5. Convergencias Sean X, Y ⊂Rnespacios de Banach y X∗el espacio dual de X. Definición A.34. (Convergencia débil) Sea {xn} ⊂ X. Se dice que {xn} ⊂ Xconverge débilmente a x∈Xsi l´ım n→∞ hf, xni=hf, xi,∀f∈X∗. Se denota con xn* x. Definición A.35. (Convergencia fuerte) Sea {xn} ⊂ X. Se dice que {xn} ⊂ Xconverge fuertemente a x∈Xsi l´ım n→∞ ||xn−x|| → 0. Proposición A.12. Sea {xn} ⊂ X. Se tiene que l´ım n→∞ ||xn−x|| → 0 =⇒xn* x. Teorema A.12. (Banach-Steinhaus) Sea (Ti∈I)una familia no necesariamente numerable de operadores lineales y continuos de X en Y tal que sup i∈I ||Tix|| <∞,∀x∈X. APÉNDICE A. ELEMENTOS DEL ANÁLISIS FUNCIONAL 96 Entonces, sup i∈I ||Ti||L(X,Y )<∞,∀x∈X i.e., ∃c∈Rtal que ||Tix|| < c||x||,∀x∈X, ∀i∈I. Definición A.36. Sea Hespacio de Hilbert, y xn⊂H, y x∈H. Entonces xn* x ⇔l´ım n→∞ (v, xn)→(v, x)∀v∈H. Además, si xn* u yvn⊂Hes tal que vn→v, entonces l´ım n→∞(vn, xn)→(v, x). Apéndice B Código FreeFem++ B.1. Simulaciones Elasticidad Código Viga 2D Acero 1// SUSANA PEREZ MARTINEZ 2//------------------------------------------------------- 3// Discretizacion por EF P1 ECUACION DE ELASTICIDAD LINEAL 4// Material acero 5//------------------------------------------------------- 6//------------------------------------------------------- 7 8// Parametros 9real rho = 7700; // rho =7700 kg/m ^3 $. 10 real mu = 11.05 e10; // mu = 11.05 e10 N/m ^2 11 real lambda = 14.06 e10 ; // lambda = 14.06 e10 N/m^2 , 12 real gravity = -9.81 , twomul =2* mu + lambda ; // fuerza de gravedad sobre el cuerpo (9 ,81 N) , Optimizacion 13 real n = 50; // numero de elementos para x 14 real m = 5; // numero de elementos para y 15 cout << " lambda ,mu , gravity =" << lambda << " " << mu << " " << gravity << endl;// indica la salida por defecto en pantalla o archivo 16 // 17 // Generacion del mallado 18 int [int] labs =[1 ,1 ,1 ,2]; //numera los bordes 19 mesh Th = square (n ,m ,[ x *10 ,y], label = labs ); // square (n , m , [ x0 +(x1 - x0 ) *x, y0 +( y1 - y0 )*y ]) 20 fespace Vh (Th ,[P1 , P1 ]); 21 // 22 // Problema 23 func A = [[ twomul , lambda , 0. ], 24 [lambda , twomul , 0. ], 25 [ 0. , 0. , mu ]] ; 26 27 macro epsV(u1 , u2 ) [dx(u1) ,dy(u2) ,dy(u1)+dx( u2 )] // EOM 28 macro div(u1 ,u2 ) ( dx(u1)+dy( u2) ) // EOM 29 Vh [ u1 , u2 ], [ v1 , v2 ]; 97 APÉNDICE B. CÓDIGO FREEFEM++ 98 30 problem Lame ([u1 , u2 ] ,[ v1 , v2 ])= int2d ( Th )( [ dx ( u1 ) ,dy (u2) ,dy( u1 )+dx(u2)] ’*A* [dx(v1) ,dy(v2) ,dy(v1)+dx( v2)] ) 31 - int2d ( Th) ( rho * gravity * v2 ) + on (2 , u1 =0 , u2 =0) ; 32 // Resuelve Problema 33 Lame; 34 // Dibujo del mallado 35 real dmax= u1 []. linfty , coef= 3/ dmax ; 36 cout << " max deplacement = " << dmax << " coef " << coef << endl; 37 mesh Thm = change ( movemesh (Th ,[ x+ u1 *coef ,y+ u2 *coef ]) ,fregion =1) ; 38 plot (Th , Thm ,wait =1 ,cmm =" coef amplification = "+ coef); Código Viga 2D Cobre 1// SUSANA PEREZ MARTINEZ 2//--------------------------------- 3// Discretizacion por EF P1 ECUACION DE ELASTICIDAD LINEAL 4// Material Cobre 5//----------------------------------- 6//----------------------------------- 7 8// Parametros 9real rho = 8960; // rho =8960 kg/m ^3 $. 10 real mu = 15.110 e10 ;// mu = 15.110 e10 N/m ^2 11 real lambda = 32.110 e10 ; // lambda = 32.110 e10 N/m ^2 , 12 real gravity = -9.81 , twomul =2* mu + lambda ; // fuerza sobre el cuerpo , Optimiszacion 13 real n = 50; // numero de elementos para x 14 real m = 5; // numero de elementos para y 15 cout << " lambda ,mu , gravity =" << lambda << " " << mu << " " << gravity << endl;// indica la salida por defecto en pantalla o archivo 16 // 17 // Generacion del mallado 18 int [int] labs =[1 ,1 ,1 ,2]; //numera los bordes 19 mesh Th = square (n ,m ,[ x *10 ,y], label = labs ); // square (n, m , [ x0 +(x1 - x0 ) *x, y0 +( y1 - y0 )*y ]) 20 fespace Vh (Th ,[P1 , P1 ]); 21 // 22 // Problema 23 func A = [[ twomul , lambda , 0. ], 24 [lambda , twomul , 0. ], 25 [ 0. , 0. , mu ]] ; 26 27 macro epsV(u1 , u2 ) [dx(u1) ,dy(u2) ,dy(u1)+dx( u2 )] // EOM 28 macro div(u1 ,u2 ) ( dx(u1)+dy( u2) ) // EOM 29 Vh [ u1 , u2 ], [ v1 , v2 ]; 30 problem Lame ([u1 , u2 ] ,[ v1 , v2 ])= int2d ( Th )( [ dx ( u1 ) ,dy (u2) ,dy( u1 )+dx(u2)] ’*A* [dx(v1) ,dy(v2) ,dy(v1)+dx( v2)] ) 31 - int2d ( Th) ( rho * gravity * v2 ) + on (2 , u1 =0 , u2 =0) ; 32 // Resolucion del Problema 33 Lame; 34 // Dibujo del mallado 35 real dmax= u1 []. linfty , coef= 3/ dmax ; APÉNDICE B. CÓDIGO FREEFEM++ 99 36 cout << " max deplacement = " << dmax << " coef " << coef << endl; 37 mesh Thm = change ( movemesh (Th ,[ x+ u1 *coef ,y+ u2 *coef ]) ,fregion =1) ; 38 plot (Th , Thm ,wait =1 ,cmm =" coef amplification = "+ coef); Código Viga 2D Hormigón 1// SUSANA PEREZ MARTINEZ 2//----------------------------------------- 3// Discretizacion por EF P1 ECUACION DE ELASTICIDAD LINEAL 4// Material Hormigon 5//------------------------------------------- 6//------------------------------------------ 7 8// Parametros 9real rho = 2300; // rho =2300 kg/m ^3 $. 10 real mu = 1.125 e10;// mu =1.125 e10 N/m ^2 11 real lambda = 0.75 e10 ; // lambda = 0.75 e10 N/m^2 , 12 real gravity = -9.81 , twomul =2* mu + lambda ; // fuerza sobre el cuerpo , Optimisation 13 real n = 50; // numero de elementos para x 14 real m = 5; // numero de elementos para y 15 cout << " lambda ,mu , gravity =" << lambda << " " << mu << " " << gravity << endl;// indica la salida por defecto en pantalla o archivo 16 // 17 // Generacion del mallado 18 int [int] labs =[1 ,1 ,1 ,2]; //numera los bordes 19 mesh Th = square (n ,m ,[ x *10 ,y], label = labs ); // square (n, m , [ x0 +(x1 - x0 ) *x, y0 +( y1 - y0 )*y ]) 20 fespace Vh (Th ,[P1 , P1 ]); 21 // 22 // Problema 23 func A = [[ twomul , lambda , 0. ], 24 [lambda , twomul , 0. ], 25 [ 0. , 0. , mu ]] ; 26 27 macro epsV(u1 , u2 ) [dx(u1) ,dy(u2) ,dy(u1)+dx( u2 )] // EOM 28 macro div(u1 ,u2 ) ( dx(u1)+dy( u2) ) // EOM 29 Vh [ u1 , u2 ], [ v1 , v2 ]; 30 problem Lame ([u1 , u2 ] ,[ v1 , v2 ])= int2d ( Th )( [ dx ( u1 ) ,dy (u2) ,dy( u1 )+dx(u2)] ’*A* [dx(v1) ,dy(v2) ,dy(v1)+dx( v2)] ) 31 - int2d ( Th) ( rho * gravity * v2 ) + on (2 , u1 =0 , u2 =0) ; 32 // Resolucion del Problema 33 Lame; 34 // Dibujo del mallado 35 real dmax= u1 []. linfty , coef= 3/ dmax ; 36 cout << " max deplacement = " << dmax << " coef " << coef << endl; 37 mesh Thm = change ( movemesh (Th ,[ x+ u1 *coef ,y+ u2 *coef ]) ,fregion =1) ; 38 plot (Th , Thm ,wait =1 ,cmm =" coef amplification = "+ coef); APÉNDICE B. CÓDIGO FREEFEM++ 100 Código Cubo 1load " msh3 " 2load " medit " 3// Funciones basicas para construir una malla regular de un cubo 4 5mesh3 Cube (NN ,BB ,L); 6-- construye la superficie mallada de una caja en 3D 7where : for exqmple: 8int [int] NN =[ nx ,ny , nz ]; // the number of seg in the 3 direction 9real [int,int] BB=[[xmin , xmax ],[ ymin , ymax ] ,[ zmin , zmax ]]; // bounding bax 10 int [int,int ] L=[[1 ,2] ,[3 ,4] ,[5 ,6]]; // the label of the 6 face left ,right , front , back , down , right 11 */ 12 func mesh3 Cube ( int [int] & NN , real [int,int] &BB ,int [int,int] & L) 13 { 14 // first build the 6 faces of the hex. 15 real x0 =BB (0 ,0) ,x1= BB (0 ,1); 16 real y0 =BB (1 ,0) ,y1= BB (1 ,1); 17 real z0 =BB (2 ,0) ,z1= BB (2 ,1); 18 19 int nx= NN [0] , ny= NN [1] , nz= NN [2]; 20 mesh Thx = square(nx ,ny ,[x0+(x1 -x0)*x,y0+(y1 -y0)*y ]) ; 21 22 int [int] rup =[0 , L (2 ,1) ], rdown =[0 ,L (2 ,0) ], 23 rmid =[1 ,L (1 ,0) , 2,L (0 ,1) , 3, L (1 ,1) , 4, L (0 ,0) ]; 24 mesh3 Th= buildlayers (Thx ,nz , zbound =[ z0 , z1 ], 25 labelmid =rmid , labelup = rup , 26 labeldown = rdown ); 27 28 return Th; 29 } Código Viga en 3D 1// Susana Perez Martinez 2// Septiembre de 2019 3//--------------------------------------------------------------- 4//--------------------------------------------------------------- 5// 6// Problema de Elasticidad Lineal en 3D ( Ecuaciones de Lame) 7 8load " medit " 9// Definimos la malla : 10 include "cube. idp";// Incluimos la definicion del cubo 11 int [int] Nxyz =[20 ,5 ,5]; 12 real [int,int] Bxyz=[[0.,5.],[0.,1.],[0.,1.]]; 13 int [int,int ] Lxyz =[[1 ,2] ,[2 ,2] ,[2 ,2]]; 14 mesh3 Th= Cube (Nxyz , Bxyz , Lxyz ); 15 // Definimos los parametros : 16 real E = 21.5 e4; APÉNDICE B. CÓDIGO FREEFEM++ 101 17 real sigma = 0.28; 18 real mu = E /(2*(1+ sigma )); 19 real lambda = E * sigma /((1+ sigma ) *(1 -2* sigma )) ; 20 real gravity = -9.81; 21 // Discretizamos el espacio con Elementos Finitos tipo P1 22 fespace Vh (Th ,[P1 ,P1 ,P1 ]); 23 Vh [u1 ,u2 ,u3], [v1 ,v2 ,v3]; 24 cout << " lambda ,mu , gravity =" << lambda << " " << mu << " " << gravity << endl; 25 26 real sqrt2 = sqrt (2.) ; 27 macro epsilon (u1 ,u2 , u3 ) [ dx (u1),dy(u2) ,dz(u3) ,( dz (u2)+dy(u3))/ sqrt2 ,( dz (u1) +dx( u3 ))/ sqrt2 ,( dy ( u1 )+dx( u2 ))/ sqrt2 ] // EOM 28 macro div(u1 ,u2 , u3 ) ( dx( u1 )+dy (u2)+ dz (u3) ) // EOM 29 // Resolvemos el problema 30 solve Lame ([u1 ,u2 ,u3] ,[v1 ,v2 ,v3])= 31 int3d ( Th )( 32 lambda *div (u1 ,u2 ,u3)*div(v1 ,v2 , v3 ) 33 +2.* mu *( epsilon (u1 ,u2 ,u3)’* epsilon (v1 ,v2 ,v3 ) ) // ’) 34 ) 35 - int3d ( Th) ( gravity *v3 ) 36 + on (1 ,u1 =0 ,u2 =0 ,u3 =0) 37 ; 38 real dmax= u1 []. max ; 39 cout << " max deplacement = " << dmax << endl; 40 // ut << " maximo desplazamiento ="< < dmax << endl ; 41 real coef= 0.1/dmax; 42 int [int] ref2 =[1 ,0 ,2 ,0]; 43 mesh3 Thm = movemesh3 (Th , transfo =[x +u1* coef ,y+ u2 * coef ,z+ u3 * coef ], label = ref2 ); 44 Thm = change ( Thm , label = ref2 ); 45 plot (Th ,Thm , wait =1 , cmm =" coef amplification = "+ coef ); B.2. Simulaciones Convección-Difusión Código Difusión Madera 1// Susana Perez Martinez 2// 3 4// Problema lineal de calor con condiciones de Dirichlet no homogeneas y 5condiciones de flujo tipo Neuman no homogeneas 6 7// 8// Parametros : 9real kappa = 0.17; // conductividad termica de Madera de Roble 10 real theta = 25; // temperatura inicial en el contorno 11 real n = 100; // numero de elementos finitos para Gamma1 12 real m = 60; // numero de elementos finitos para Gamma2 13 real dt = 1.0; // salto de tiempo 14 // 15 // Se definen los bordes 16 border Gamma1 (t =0 ,2* pi) {x =2.5* cos(t); y =1.75* sin( t);} APÉNDICE B. CÓDIGO FREEFEM++ 102 17 border Gamma2 (t=0 , 1) {x = -0.5+t; y =0.2;} 18 // Se define la malla 19 mesh Th = buildmesh ( Gamma1 (n)+ Gamma2 (m)); 20 // Se dibuja la malla 21 plot (Th , wait =1 , ps="Th. eps "); 22 // 23 // Se discretiza el espacio de funciones 24 fespace Vh(Th , P1 ); 25 Vh uh , uold , vh; 26 // 27 // Problema 28 uh = theta ; uold = uh; 29 problem heatstep (uh , vh , solver =LU ) = 30 int2d ( Th )(uh * vh / dt ) 31 - int2d ( Th )( uold * vh / dt ) 32 + int2d (Th)( kappa * dx (uh )* dx (vh) + kappa *dy ( uh )* dy (vh) ) 33 + int1d (Th , Gamma1 )( kappa *vh) 34 + on (Gamma2 , uh= theta); // Condiciones de Dirichlet en Gamma2 35 // Hallamos la solucion para 6 tiempos diferentes 36 for (int it =0;it <4; it ++) { 37 // Se renueva la temperatura con respecto al tiempo anterior 38 uold = uh; 39 // Se realiza un salto de tiempo y se resuelve 40 heatstep ; 41 // Se dibuja la solucion 42 plot (uh , nbiso =50 , fill =1 , value =1 , wait =1 , ps=" res"+it+". eps "); 43 44 } Código Difusión en Plata con salto de tiempo dt = 1 1// Susana Perez Martinez 2// 3 4// Problema lineal de calor con condiciones de Dirichlet no homogeneas y 5condiciones de flujo tipo Neuman no homogeneas 6 7// 8// Parametros : 9real kappa = 419; // conductividad termica de Plata 10 real theta = 25; // temperatura inicial en el contorno 11 real n = 100; // numero de elementos finitos para Gamma1 12 real m = 60; // numero de elementos finitos para Gamma2 13 real dt = 1; // salto de tiempo 14 // 15 // Se definen los bordes 16 border Gamma1 (t =0 ,2* pi) {x =2.5* cos(t); y =1.75* sin( t);} 17 border Gamma2 (t=0 , 1) {x = -0.5+t; y =0.2;} 18 // Se define la malla 19 mesh Th = buildmesh ( Gamma1 (n)+ Gamma2 (m)); 20 // Se dibuja la malla APÉNDICE B. CÓDIGO FREEFEM++ 103 21 plot (Th , wait =1 , ps="Th. eps "); 22 // 23 // Se discretiza el espacio de funciones 24 fespace Vh(Th , P1 ); 25 Vh uh , uold , vh; 26 // 27 // Problema 28 uh = theta ; uold = uh; 29 problem heatstep (uh , vh , solver =LU ) = 30 int2d ( Th )(uh * vh / dt ) 31 - int2d ( Th )( uold * vh / dt ) 32 + int2d (Th)( kappa * dx (uh )* dx (vh) + kappa *dy ( uh )* dy (vh) ) 33 + int1d (Th , Gamma1 )( kappa *vh) 34 + on (Gamma2 , uh= theta); // Condiciones de Dirichlet en Gamma2 35 // Hallamos la solucion para 4 tiempos diferentes 36 for (int it =0;it <4; it ++) { 37 // Se renueva la temperatura con respecto al tiempo anterior 38 uold = uh; 39 // Se realiza un salto de tiempo y se resuelve 40 heatstep ; 41 // Se dibuja la solucion 42 plot (uh , nbiso =50 , fill =1 , value =1 , wait =1 , ps=" res"+it+". eps "); 43 44 } Código Difusión en Aire 1// Susana Perez Martinez 2// 3// Problema lineal de calor con condiciones de Dirichlet no homogeneas y 4// condiciones de flujo tipo Neuman no homogeneas 5// 6real kappa = 0.024; // conductividad termica del Aire 7real theta = 25; // temperatura inicial en el contorno 8real n = 100; // numero de elementos finitos para Gamma1 9real m = 60; // numero de elementos finitos para Gamma2 10 real dt = 1; // salto de tiempo 11 // 12 // Se definen los bordes 13 border Gamma1 (t =0 ,2* pi) {x =2.5* cos(t); y =1.75* sin( t);} 14 border Gamma2 (t=0 , 1) {x = -0.5+t; y =0.2;} 15 // Se define la malla 16 mesh Th = buildmesh ( Gamma1 (n)+ Gamma2 (m)); 17 // Se dibuja la malla 18 plot (Th , wait =1 , ps="Th. eps "); 19 // 20 // Se discretiza el espacio de funciones 21 fespace Vh(Th , P1 ); 22 Vh uh , uold , vh; 23 // 24 // Problema APÉNDICE B. CÓDIGO FREEFEM++ 104 25 uh = theta ; uold = uh; 26 problem heatstep (uh , vh , solver =LU ) = 27 int2d ( Th )(uh * vh / dt ) 28 - int2d ( Th )( uold * vh / dt ) 29 + int2d (Th)( kappa * dx (uh )* dx (vh) + kappa *dy ( uh )* dy (vh) ) 30 + int1d (Th , Gamma1 )( kappa *vh) 31 + on (Gamma2 , uh= theta); // Condiciones de Dirichlet en Gamma2 32 // Hallamos la solucion para 4 tiempos diferentes 33 for (int it =0;it <4; it ++) { 34 // Se renueva la temperatura con respecto al tiempo anterior 35 uold = uh; 36 // Se realiza un salto de tiempo y se resuelve 37 heatstep ; 38 // Se dibuja la solucion 39 plot (uh , nbiso =50 , fill =1 , value =1 , wait =1 , ps=" res"+it+". eps "); 40 41 } Código Difusión en Habitación y Pared 1// Susana Perez Martinez 2// Temperatura en las paredes de una habitacion con radiador 3// La pared esta compuesta por dos materiales : uno caro y uno barato. 4// Parametros 5real Lx = 4, Ly =5 , lr =1 , lx =0.3; // longitud horizontal , y vertical de la planta , del radiador y anchura de la pared . 6real Tr =20 , Text = -5; // temperaturas del radiador y externa respectivamente 7real hy = 0.1; 8real hx = 0.1; 9real mu = 0.2; 10 real kappaa = 0.032; // conductividad del aire 11 real kappac = 0.017; // conductividad del material barato 12 real kappae = 0.008; // conductividad del material caro 13 real V = lx *Ly ; // volumen de la pared 14 real Ve = mu * V; // volumen del material caro 15 real xem = Lx + hx; 16 real deltaxe = Ve / (Ly -2* hy ); 17 // Definimos los bordes de la habitacion 18 border c1( t=0 , Lx ){x =t; y =0;} // borde inferior 19 border c1p (t=Lx , Lx +lx) {x=t; y =0;} // borde inferior pared 20 border cext(t=0 ,Ly){x=Lx+lx; y=t;} // borde exterior ( pared derecha ) 21 border c3p (t= Lx +lx , Lx ){x =t; y= Ly ;} // borde superior pared 22 border c3( t=Lx ,0) {x =t; y= Ly ;} // borde superior pared 23 border c4(t=Ly , Ly /2+ lr /2){x =0; y=t ;} // borde superior al radiador pared izquierda 24 border cr(t= Ly /2+ lr/2 ,Ly /2 -lr /2){x =0; y=t ;} // borde donde esta el radiador 25 border c5(t= Ly /2 -lr /2 ,0) {x=0; y=t ;} // borde inferior al radiador pared izquierda