pdetool-440

ESCUELA TÉCNICA SUPERIOR DE INGENIEROS INDUSTRIALES

Dpto. Lenguajes y Ciencias de la Computación

Ciencia de la Computación e Ingeniería Artificial

PROYECTO FIN DE CARRERA

MÉTODO DE VOLÚMENES FINITOS PARA

PROBLEMAS NO LINEALES

Autor: Pablo Udias de la Mora García Director: Dr. Francisco R. Villatoro Machuca Titulación: Ingeniero en Organización Industrial

MÁLAGA, mayo de 2006

Agradecimientos

A mis padres, Paul e Isabel.

A Clarisse.

Al Dr. Francisco Villatoro que, no solo ha guiado la realizacion de este proyecto sinoque me ha permitido dar mis primeros pasos en el mercado laboral y me ha servido deconstante estımulo para progresar.

Al Dr. Jose Aldana por brindarme la oportunidad de continuar mi formacion con unabeca TIC y explorar nuevos campos de investigacion.

A mis companeros del grupo KHAOS del Departamento de Lenguajes y Ciencias dela Computacion de la E.T.S.I. Informatica: Ma del Mar Roldan, Ma del Mar Rojano,Chema, Dani, Raul, Isma, Monsef, Ivan, Cesar y Sergio por los buenos 9 meses que mehan hecho pasar en el laboratorio.

MalagaMayo de 2006 Pablo Udias de la Mora

i

Indice general

1. Introduccion y objetivos 1

1.1. Preliminares . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

1.2. Objetivos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

1.3. Contenidos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

2. Desarrollo del metodo de los volumenes finitos con Matlab 5

2.1. Modelo de un problema generico en derivadas parciales parabolico . . . . 6

2.2. Metodo de los volumenes finitos . . . . . . . . . . . . . . . . . . . . . . . 6

2.2.1. Volumenes de control con circuncentro . . . . . . . . . . . . . . . 8

2.2.2. Volumenes de control con baricentro . . . . . . . . . . . . . . . . 11

2.3. Estructura de la informacion de la geometrıa . . . . . . . . . . . . . . . . 13

2.3.1. PDEToolbox de Matlab . . . . . . . . . . . . . . . . . . . . . . . 13

2.3.2. Nueva estructura de la informacion para FVM con circuncentro . 14

2.3.3. Nueva estructura de la informacion para FVM con baricentro . . 17

2.4. Implementacion en Matlab . . . . . . . . . . . . . . . . . . . . . . . . . . 17

2.4.1. Codigos Matlab para FVM con circuncentro . . . . . . . . . . . . 17

2.4.2. Codigos Matlab para FVM con baricentro . . . . . . . . . . . . . 22

2.5. Validacion del metodo . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

2.6. FVM baricentrico frente a FVM circuncentrico en problemas lineales ge-nerales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

3. Problema no lineal del calentamiento de preformas 34

ii

INDICE GENERAL iii

3.1. Propiedades termicas del aire y del PMMA . . . . . . . . . . . . . . . . . 37

3.1.1. Conductividad termica del PMMA . . . . . . . . . . . . . . . . . 38

3.1.2. Calor especıfico del PMMA . . . . . . . . . . . . . . . . . . . . . 41

3.1.3. Otras propiedades del PMMA . . . . . . . . . . . . . . . . . . . . 42

3.2. Formulacion del modelo del problema . . . . . . . . . . . . . . . . . . . . 43

4. Metodo de volumenes finitos para problemas termicos no lineales 46

4.1. Metodo de los volumenes finitos en problemas no lineales . . . . . . . . . 46

4.1.1. Ecuaciones para los VC pertenecientes a TΩ . . . . . . . . . . . . 48

4.1.2. Ecuaciones para los VC pertenecientes a T∂Ω . . . . . . . . . . . . 49

4.1.3. Resolucion del sistema . . . . . . . . . . . . . . . . . . . . . . . . 51

4.1.4. Metodo implıcito-explıcito . . . . . . . . . . . . . . . . . . . . . . 52

4.2. Codigos Matlab para la resolucion de PDEs no lineales mediante FVM . 53

5. Aplicacion y resultados de FVM en PDEs no lineales para produccionde PCF 56

5.1. Validacion del metodo desarrollado . . . . . . . . . . . . . . . . . . . . . 56

5.2. Comparativa de los resultados del modelo no lineal con los resultados ex-perimentales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

5.3. Analisis del calentamiento de la preforma . . . . . . . . . . . . . . . . . . 61

6. Conclusiones 66

6.1. Implementacion del metodo de volumenes finitos . . . . . . . . . . . . . . 66

6.2. Resultados obtenidos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67

6.2.1. Resultados en problemas lineales . . . . . . . . . . . . . . . . . . 67

6.2.2. Resultados en problemas no lineales . . . . . . . . . . . . . . . . . 68

6.3. Nuevas lıneas de investigacion . . . . . . . . . . . . . . . . . . . . . . . . 69

Bibliografıa 70

A. Contenido del CD 75

INDICE GENERAL iv

A.1. Contenidos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75

A.2. Codigos implementados . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76

B. Codigos en Matlab 77

B.1. Distancia mınima . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77

B.2. Circuncentro . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78

B.3. Sistema EDO baricentro . . . . . . . . . . . . . . . . . . . . . . . . . . . 79

C. Calidad del mallado 81

D. Artıculo enviado para publicacion 83

Capıtulo 1

Introduccion y objetivos

1.1. Preliminares

El mundo se mueve siguiendo leyes fısicas que han sido modelizadas durante siglospor los grandes pensadores. Multiples campos de investigacion intentan descifrar la com-prension de los sucesos por medio de las matematicas. En matematicas, una ecuacion enderivadas parciales (PDE) es una relacion que involucra a una funcion con varias variablesindependientes y sus derivadas parciales con respecto de estas variables. Las ecuacionesen derivadas parciales son utilizadas para formular y resolver problemas como pueden serla propagacion del sonido y del calor, electrostatica, electrodinamica, y elasticidad. Enterminos generales se utilizan para los procesos que se distribuyen en el espacio o en elespacio y tiempo. Problemas de muy distinto tipo pueden tener las mismas formulacionesmatematicas.

En economıa siguen sin existir principios fundamentales como en fısica, por lo queno se conocen las formas generales de las ecuaciones que describen la evolucion de lossistemas economicos [1]. Las ecuaciones son hipoteticas y requieren de la verificacion paracada caso particular. No obstante, grandes avances se estan haciendo en la modelizaciony es cada vez mas frecuente el encontrarnos con problemas que hacen uso de PDEs paradescribir estos sistemas. Es el caso del intercambio de capital y demanda del consumidoren la tecnologıa bajo la accion de cambios en los beneficios [2]. Tambien se modela coneste tipo de ecuaciones los incrementos de retorno, el tamano de la manufactura depoblaciones, y los costes de transporte en la formacion de una economıa espacial [3] eincluso para evaluar las style options europeas en multiples fondos [4].

El metodo de separacion de variables puede permitir obtener soluciones de PDEslineales sujetas a condiciones iniciales y/o de contorno. Si el dominio es finito o periodico,una suma infinita de soluciones como las series de Fourier pueden ser adecuadas. En elcaso de las PDEs no lineales no existen metodos generales que las resuelvan, no obstante

1

CAPITULO 1: INTRODUCCION Y OBJETIVOS 2

existen tecnicas especıficas que resuelven problemas concretos como el potente metodo delprincipio-h (homotopy-principle) para ecuaciones indeterminadas o la teorıa de Riquier-Janet [5] para obtener informacion acerca de sistemas sobredimensionados.

Como alternativa, se suele recurrir en la mayor parte de los problemas ingenieriles alanalisis numerico. Para la resolucion se hace uso de computadores e incluso de supercom-putadores para los problemas mas complejos que resuelven los metodos implementados yque pueden ir desde sencillos esquemas de diferencias finitas hasta los metodos recientesde elementos finitos multimalla. Entre estas tecnicas se encuentra el metodo de volumenesfinitos (FVM) que presenta ventajas que citamos a continuacion.

Los esquemas de volumenes finitos fueron desarrollados por ingenieros con el fin deestudiar complejos fenomenos fısicos acoplados, modelados por PDEs, en donde la con-servacion de las cantidades extensivas (como son la masa, la energıa, ...), han de serrigurosamente respetadas por la solucion aproximada. Una ventaja de este tipo de esque-mas es el uso de una amplia variedad de mallas que discretizan el dominio de estudio entrelos que podemos citar aquellos que hacen uso de los triangulos originales de la malla comovolumenes de control (VC) llamados metodos centrados en la celda (cell-centered) [6] ylos que recurren al uso de la malla dual como volumenes de control llamados metodoscentrados en los vertices (vertex-centered) [7]. Entre estos ultimos destacan la utilizacionde la malla Voronoi construida mediante teselacion del dominio [8, 9]. Indistintamente dela discretizacion espacial llevada a cabo, la idea basica del metodo es integrar la PDE encada volumen de control y aproximar los flujos a traves de los lımites de estos VC. Paraproblemas parabolicos el FVM es un metodo para discretizacion espacial de ecuaciones enderivadas parciales que combinado con una discretizacion temporal permite determinaruna aproximacion a la solucion resolviendo un sistema de ecuaciones algebraicas.

Varias variantes del FVM se pueden implementar usando un esquema centrado en lacelda. En efecto, dada una malla de triangulos no estructurada, como puede ser la mallade Delaunay que maximiza el angulo mınimo de los triangulos en el proceso de trian-gularizacion [10], se puede implementar el FVM mediante el uso del centroide tambienconocido por baricentro y centro de gravedad o el uso del circuncentro como nodo repre-sentativo de los VC. No obstante, la existencia de triangulos obtusos en la malla generaproblemas de computo [11]. En este marco se engloba el proyecto del que se proceden adescribir los objetivos.

1.2. Objetivos

El objetivo del presente proyecto es el desarrollo de un codigo en Matlab basado entecnicas matematicas que permitan resolver problemas en derivadas parciales no linealesen dominios complejos por su geometrıa y heterogeneidad. Centramos todo el estudio en el

CAPITULO 1: INTRODUCCION Y OBJETIVOS 3

metodo de volumenes finitos, con las distintas variantes que dispone, para la discretizacionespacial. De esta forma se resuelven problemas ingenieriles actuales.

En particular se procedera a calcular la respuesta termica de ciertos solidos no ho-mogeneos cuyas propiedades termicas son fuertemente dependientes de la temperatura yde la posicion y son por tanto no lineales. Para ello se hara una semidiscretizacion espacialmediante el metodo de los volumenes finitos y en tiempo mediante la regla del trapecio.Las no linealidades de la ecuacion obligaran a la resolucion mediante metodos iterativosque es de esperar converjan a la solucion. Los resultados se obtendran mediante un codi-go que partira de una malla de triangulos no estructurada, adaptada a las geometrıasa estudiar, que incluyen estructuras compuestas de objetos elementales como cırculos yparalelepıpedos.

En el desarrollo del metodo numerico se debera interpretar la codificacion de estageometrıa, en particular del mallado generado por el metodo de Delaunay (mediante lasherramientas de Matlab), y tratarla de forma adecuada ya que Matlab prepara dichageometrıa para una implementacion del metodo de los elementos finitos y en nuestrocaso estamos interesados en su aplicacion al metodo de los volumenes finitos. Se obten-dra como resultado la respuesta temporal en cada nodo de la malla. Mediante distintasrepresentaciones graficas se hara un estudio de la influencia de la variacion del valor de laspropiedades termicas de los materiales y de las condiciones de contorno. Se analizaran,igualmente, otros aspectos como la modificacion de la geometrıa del dominio de estudioy su influencia sobre la transferencia de calor en el medio no lineal.

Un objetivo subsidiario del proyecto es la elaboracion de varios artıculos de investiga-cion para publicar en revistas con ındice de impacto. Para ello se ha realizado un amplioestudio de la bibliografıa existente sobre temas especıficos como el metodo de volumenesfinitos con uso de centroide o de circuncentro en la resolucion de PDEs. Con los resultadosobtenidos en el proyecto, se ha enviado a publicar un articulo.

1.3. Contenidos

La estructura del proyecto es la siguiente: en el capıtulo segundo se procede a la reso-lucion numerica de PDEs con coeficientes constantes. En primer lugar se muestra el tipode problema a resolver mediante la descripcion del modelo. Posteriormente, en el aparta-do 2.2 se realiza el desarrollo matematico de las dos variantes de volumenes finitos parala resolucion del modelo. A continuacion se describe en profundidad el esquema de imple-mentacion del metodo para geometrıas de computo complejas, ası como los detalles de lamalla de Delaunay utilizada. Tambien se describe parte del codigo utilizado en lenguajeMatlab, afın de mostrar las complejidades de la implementacion. En el apartado 2.5, seprocede a validar el metodo desarrollado mediante la comparacion de los resultados conla solucion analıtica en un problema sencillo y la solucion dada por un software numerico

4

comercial. En el ultimo apartado del capıtulo 2, se realiza un analisis detallado de losresultados generados por las dos variantes FVM, es decir con uso de circuncentro y debaricentros como nodos de los volumenes de control, obteniendose datos no descritos enla literatura. A partir de los datos de la comparativa de las dos versiones de FVM seredacta un artıculo en colaboracion con el Dr. Francisco Villatoro.

En el capıtulo tercero, se describe una nueva tecnologıa que esta surgiendo en el mer-cado como son las fibras de cristales fotonicas (PCF) sobre las que vamos a implementarel metodo de volumenes finitos. En particular se analizan las propiedades termicas nolineales de estas fibras que influyen decisivamente en la etapa de calentamiento de lapreforma, etapa previa al estirado de fibras. Para ello se recopila la informacion existenteen la literatura en donde se describe de forma teorica y experimental estas caracterısticastermicas. Posteriormente, se plantea el modelo de un problema matematico bidimensio-nal no lineal a resolver, que describe el calentamiento de la preforma. El modelo se haadimensionalizado y cuenta con mejoras con respecto a los modelos lineales descritos enotros artıculos que no se ajustan con los resultados experimentales llevados a cabo endistintos laboratorios.

El metodo numerico iterativo basado en volumenes finitos es desarrollado en el capıtu-lo cuarto. En particular este metodo permite resolver problemas regidos por ecuacionesen derivadas parciales parabolicas con coeficiente y condiciones de contorno no lineales,como es el problema general del calor. Para ello, se linealizan a trozos las propiedadestermicas de los materiales y se procede de forma iterativa a su resolucion mediante uncontrol de la tolerancia del error. Otras tecnicas igualmente son descritas.

El capıtulo quinto trata sobre el analisis de resultados del calentamiento de preformasde PCF obtenidos mediante los metodos desarrollados en secciones anteriores. Se validael metodo de una forma clasica al no disponer de software que permita obtener resultadospara el problema no lineal con geometrıa compleja que se estudia. Los resultados obtenidosson comparados con los generados por otros autores que hicieron uso de modelos menosexactos. Tambien, se realiza un analisis del proceso de calentamiento de la preforma dePCF y se estudia la influencia de los distintos parametros, lo que permite identificar lasposibles mejoras en la produccion de esta nueva tecnologıa.

El capıtulo de conclusion muestra que el desarrollo de este PFC ha aportado un valoranadido a los estudios realizados sobre nuevas tecnologıas como son la produccion en ma-sa de PCF mediante la determinacion de un modelo mas semejante a la realidad fısica yel desarrollo de un metodo numerico de volumenes finitos adaptado a su resolucion. Asi-mismo se contribuye a la actividad investigadora mediante la elaboracion de dos artıculosque aportan nuevos conceptos no descritos anteriormente.

Los anexos hacen referencia a los codigos desarrollados, al contenido del CD-Rom quese adjunta y a la importancia de la calidad de la malla en el calculo numerico de FVM.Por ultimo, se adjunta un artıculo enviado para su posterior publicacion por el autor ydirector del presente proyecto.

Capıtulo 2

Desarrollo del metodo de losvolumenes finitos con Matlab

A diferencia de los complejos codigos utilizados en los paquetes informaticos comer-ciales, este proyecto proporciona una simple implementacion en Matlab para la resolucionde ecuaciones en derivadas parciales parabolicas por el metodo de los volumenes finitos.Se desarrollan dos esquemas de volumenes finitos de cuatro puntos en malla triangulartomando tanto el circuncentro como el baricentro como nodo central de los volumenes decontrol. Ambos, son metodos conservativos que presentan como ventajas la estabilidad yla convergencia. Para ello se hace una adaptacion de los datos obtenidos por la PDETool-box de Matlab en cuanto a codificacion de los datos de la discretizacion de la geometrıadel dominio de computo.

La PDEToolbox es una potente herramienta cuyos objetivos principales son (a) dis-cretizar un dominio 2D mediante triangularizacion y codificar la informacion en tresmatrices, y (b) la resolucion de PDEs mediante elementos finitos. La perspectiva de esteproyecto no es resolver todos los tipos de problemas mediante un unico codigo, sino pro-porcionar una herramienta facilmente entendible y reutilizable para otros problemas quesurjan, con el mismo espıritu de otros artıculos [12, 13].

El capıtulo se articula como sigue: el modelo de un problema regido por ecuacionesen derivadas parciales lineales y con condiciones de contorno tipo Dirichlet a resolverse describe en la seccion 2.1. El desarrollo del FVM y la discretizacion del modelo serealiza en la seccion 2.2 tanto con uso de circuncentro como de baricentro explicando lasprincipales ventajas obtenidas por cada variante. La seccion 2.3 se inicia con la descrip-cion de la estructura de la codificacion de la geometrıa por parte de la PDEToolbox deMatlab, y posteriormente se analizan las modificaciones a realizar para la adaptacion deesta informacion al FVM. El codigo implementado en Matlab para la resolucion de estemodelo mediante las dos variantes de FVM implementadas se detalla en la seccion 2.4.Posteriormente, en el apartado 2.5 se hace una comprobacion de la validez del metodo

5

CAPITULO 2: FVM con Matlab 6

desarrollado comparando las soluciones que proporciona con las soluciones analıticas ylas soluciones de otro metodo numerico implementado con un software comercial para unproblema sencillo que permita obtener soluciones analıticas. Por ultimo, la seccion 2.6muestra los resultados obtenidos mediante las dos versiones de FVM para unos parame-tros determinados del modelo en un dominio circular y se comparan con los generadospor la PDEToolbox mediante FEM. De esta forma se ajustan parametros basicos comoel numero de nodos de la malla a utilizar, comprobar que metodo es mejor, entre otros,para su posterior utilizacion en resolucion de problemas PDEs no lineales.

2.1. Modelo de un problema generico en derivadas

parciales parabolico

El programa de Matlab propuesto emplea el metodo de los volumenes finitos paracalcular una solucion numerica U , que aproxima la solucion u del problema bidimensionalde la ecuacion en derivadas parciales parabolica lineal (2.1), en el dominio computacionalΩ = Ω ∪ ∂Ω siendo Ω el interior y ∂Ω su frontera.

a∂u(~x, t)

∂t+∇ · (b∇u(~x, t)) + c u(~x, t) = d , t > 0, ~x ∈ Ω ⊂ IR2, (2.1)

Se trata de un problema de valores iniciales en donde a, b, c y d son parametros constantes.Como condicion inicial se tiene u(~x, 0) = u0. El problema se resuelve para condiciones decontorno constantes Dirichlet.

u(~x, t) = uext, t > 0, ~x ∈ ∂Ω, (2.2)

siendo uext la solucion constante en la superficie.

En la seccion 2.5 se resuelve haciendo a = 1, b = −1, c = 0 y d = 0, para intentarvalidar el FVM comparandolo con la solucion analıtica. En la seccion 2.6 se analizan endetalle los resultados de las dos variantes de FVM para este problema con valores nonulos para los parametros a, b, c y d.

Este modelo de problema PDE parabolico permite desarrollar y validar los metodosFVM implementados mediante otras opciones como son las soluciones analıticas de pro-blemas sencillos o las soluciones obtenidas mediante software comercial. Una vez validadose procedera en posteriores capıtulos a desarrollar FVM para problemas PDEs no lineales.

2.2. Metodo de los volumenes finitos

El metodo de volumenes finitos (FVM) fue introducido en el campo de la dinamicade fluidos numerica por McDonald (1971). Se da este nombre a la tecnica a traves de la


cual la formulacion integral de las leyes de conservacion son discretizadas directamenteen el espacio fısico. Segun algunos autores [14], este metodo puede considerarse como uncaso particular del esquema de diferencias finitas, para otros, lo es de elementos finitos.El dominio del problema se divide en volumenes de control triangulares con un malladode puntos asociado. Cada triangulo lleva asociado un nodo en el que se resolvera de formadiscreta el problema. Estos nodos pueden localizarse en distintas posiciones siendo lo mashabitual utilizar el circuncentro [15] o el baricentro [16]. La eleccion de estos puntos es deespecial importancia puesto que pueden simplificar de forma sensible la implementaciondel codigo.

En el dominio computacional Ω se define una malla formada por un conjunto detriangulos (volumenes de control) TΩ = T∂Ω+TΩ, donde T∂Ω es el subconjunto de triangu-los del contorno con al menos una arista perteneciente a ∂Ω y TΩ el resto de triangulos nosituados en el contorno. La nomenclatura para cada uno de los triangulos es la que sigue:los volumenes de control vienen definidos por tres vertices a, b y c; los tres trianguloscon arista comun bc, ca y ab al VC considerado tienen por centro los puntos A, B y C;las longitudes de las aristas son lab, lbc y lca;y por ultimo, las longitudes entre centros detriangulos son liA, liB y liC .

Se considera un volumen de control R (triangulo abc) con centro i, perteneciente a TΩ.Integrando la EDP en su forma conservativa, ecuacion (3.5), sobre el volumen de controlse obtiene

∫ ∫R

(a∂u

∂t+∇ · (b∇u) + cu− d)dR = 0, ∀ R ∈ TΩ. (2.3)

Aplicando el teorema de la divergencia se tiene la forma integral de la ley de conservacion∫ ∫R

(a∂u

∂t+ cu− d)dR +

∮S

b∇u · d~s = 0, (2.4)

donde d~s = ~nds, siendo ~n el vector unitario normal a la superficie (hacia fuera) y ds esel elemento diferencial de superficie.

El termino de la izquierda en la ecuacion (2.4) contiene la derivada en tiempo y sepuede evaluar asumiendo que la magnitud u en el punto i-esimo es el valor medio paratodo el volumen ∫ ∫

R

(a∂u

∂t+ cu− d)dR = (a

∂ui

∂t+ cui − d)Ai, (2.5)

siendo Ai el area del volumen de control. La segunda integral en la ecuacion (2.4) re-presenta el flujo saliente de la magnitud u para las tres fronteras del volumen de controldel punto i-esimo. Llegados aquı la eleccion del nodo de los VC que actuara como centrode la celda, en el FVM de cuatro puntos, cobra una gran importancia en cuanto a lasimplificacion del metodo. Optamos por seguir una estrategia cell-centered y por tantono hacemos uso de la malla dual. Los dos opciones elegidas son el uso sobre la malla


(a) (b)

Figura 2.1: Ejemplo de localizacion de los circuncentros de los volumenes de control (a)para triangulos agudos, (b) para triangulos obtusos.

original de Delaunay en primer lugar de nodos localizados en el circuncentro de los VCy posteriormente en el baricentro, con sus ventajas e inconvenientes respectivos.

2.2.1. Volumenes de control con circuncentro

El uso del circuncentro de los triangulos de la malla como nodo de calculo para cadavolumen de control permite una facil valoracion del miembro derecho de la ecuacion (2.4),es decir del flujo. En efecto como se sabe el circuncentro se localiza en la interseccion delas mediatrices de las aristas y por tanto la recta que une los nodos de dos VC vecinos esortogonal a una arista y por ello colineal al flujo, ver la figura 2.1. Se tiene entonces que elproducto escalar entre el gradiente de la magnitud, ∇u, y la normal a la superficie, ~n, essimplemente el valor ∇u. La eleccion del circuncentro simplifica logicamente el desarrollodel metodo que prosigue de la siguiente forma.

La integral de superficie cerrada de la ecuacion (2.4) se computa sobre cada una delas tres aristas del VC. Para la arista bc se puede expresar mediante∫

bc

b∇u · ~nds ' biA(uA − ui

liA)lbc, (2.6)

y se evalua de forma analoga en los segmentos ab y ca.

Ecuaciones para los VC pertenecientes a TΩ

Para todos los volumenes de control pertenecientes a TΩ, figura 2.2.a, la ecuacion (2.4)se transforma en


(a∂ui

∂t+ cui − d)Ai + b(

uA − ui

liA)lbc

+ b(uB − ui

liB)lca

+ b(uC − ui

liC)lab = 0, (2.7)

una ecuacion diferencial ordinaria para cada volumen de control. Siendo lab, lbc y lcalas longitudes de las aristas del volumen de control y, liA, liB y liC las distancias entre elcircuncentro del volumen de control del nodo i-esimo y los circuncentros de los volumenesde control vecinos.

Ecuaciones para los VC pertenecientes a T∂Ω

A los volumenes de control j-esimos pertenecientes a T∂Ω, se les aplica la condicionde contorno Dirichlet. Al ser la temperatura del contorno conocida, el flujo conductivo atraves de la arista frontera (bc) es

∫bc

b∇u · ~nds ' b(uext − uj

ljs1)lbc. (2.8)

siendo ljs1 la distancia mınima entre el circuncentro del volumen de control y el contorno(bc). A partir de la ecuacion (2.4), para los volumenes de control (j) con una unica aristacontorno, figura 2.2.b, se obtiene

(a∂uj

∂t+ cuj − d)Aj + b(

uC − uj

ljC)lab

+ b(uB − uj

ljB)lca − b

uj

ljs1lbc = −b

uext

ljs1lbc, (2.9)

puede ocurrir que un volumen de control tenga dos aristas contorno, figura 2.2.c, entoncesla ecuacion es

(a∂uj


uC − uj

ljC)lab

− (lbcljs1

+lcaljs2

)buj = −buext

ljs1lbc − b

uext

ljs2lca, (2.10)

siendo ljs1 y ljs2 las distancias mınimas entre el circuncentro y cada una de las aristascontorno.


(a) (b) (c)

Figura 2.2: Mallado para el metodo de volumenes finitos. Ejemplo de los tres tipos devolumenes de control, (a) interior (i), (b) con una frontera (j) y (c) con dos fronteras (j)

Resolucion del sistema

Por tanto se tiene un sistema de tantas ecuaciones diferenciales ordinarias comovolumenes de control tenga la malla. Un mallado estructurado para desarrollar el meto-do de volumenes finitos simplifica la notacion del metodo, no obstante, las geometrıascomplicadas requieren de un mallado poco regular.

El sistema EDO se representa de forma matricial

Adu

dt+ Bu = C, (2.11)

siendo u el campo escalar de la magnitud u de los volumenes de control codificado enforma de vector, A es una matriz diagonal; B es una matriz con dos o tres elementos porfila para los volumenes de control pertenecientes a T∂Ω en funcion de si tienen dos o unaarista frontera, respectivamente, y cuatro elementos para los pertenecientes a TΩ; C esun vector con valor 0 unicamente para las filas correspondientes a triangulos incluidos enTΩ.

El problema se resuelve mediante una semidiscretizacion: en espacio mediante FVMy en tiempo se hace uso de la regla del trapecio. Obteniendose para el instante de tiempo(n + 1) la solucion

un+1 = (A +4t

2B)−1(4t C + (A− 4t

2B)un), (2.12)

siendo u0 = u0 el valor inicial conocido.

El esquema desarrollado es incondicionalmente estable, esto no significa que se hallensoluciones fısicamente realistas independientemente del paso de tiempo tomado [17], porel contrario se pueden encontrar soluciones oscilatorias. La estabilidad asegura que estasoscilaciones se atenuan pero no garantiza soluciones fısicamente realistas. En efecto laregla del trapecio es un metodo A-estable pero no L-estable luego genera oscilaciones


espureas. Mediante el control del paso de tiempo, ∆t, se reduce la amplitud de estasoscilaciones. Para un pasos de tiempo inferior a un paso de tiempo crıtico, ∆t < ∆t∗, seanulan las oscilaciones. Por tanto es importante una correcta eleccion del paso de tiempo,que se toma logarıtmicamente espaciado.

Como inconvenientes a la implementacion de este metodo se tiene que para los triangu-los obtusos de la malla sus circuncentros se localizan fuera de ellos, ver la figura 2.1.b.Esta es una primera fuente de error que no afecta al resultado cuando la proporcion deestos triangulos es pequena. Tambien puede ocurrir que algun circuncentro este fuera deldominio de calculo con lo que los resultados si pueden diferir de forma mas importante.Dos sucesos fatales para el correcto funcionamiento de FVM son posibles (1) que doscircuncentros de dos triangulos vecinos se superpongan, y (2) que un circuncentro de lostriangulos pertenecientes al conjunto T∂Ω sea rectangulo y su circuncentro se localice sobrela arista frontera. Estos dos casos provocan indeterminacion en las ecuaciones (2.7), (2.9)y (2.10) siendo por tanto casos crıticos para la resolucion del problema por lo que hay quealterar de alguna forma la disposicion de los nodos de la malla. No obstante, el uso de unamalla de Delaunay genera triangulos de alta calidad, con alto ındice de equilateralidad ypor tanto solo una pequena proporcion de nodos de calculo se situan fuera de los VC. Sise presenta algun caso peligroso es suficiente con variar levemente la posicion de algunvertice de la malla original para subsanar este problema o calcular una nueva malla deDelaunay bajo otro criterio hasta que el triangulo conflictivo no vuelva a surgir.

2.2.2. Volumenes de control con baricentro

En este apartado se analiza la variante al desarrollo del metodo usando el baricentrocomo nodo de calculo de cada volumen de control en sustitucion del circuncentro. Estecambio repercute en la ecuacion (2.4) puesto que ahora el gradiente de la magnitud ∇uy el vector normal ~n ya no son colineales y por tanto su producto escalar se ve afectadopor un factor que es el coseno del angulo entre estos dos vectores∫

bc

b∇u · ~nds ' biA(uA − ui

liA) cos(αA)lbc. (2.13)


Por tanto para todos los volumenes de control pertenecientes a TΩ, la ecuacion (2.4)se transforma en

(a∂ui

∂t+ cui − d)Ai + b(

uA − ui

liA) cos(αA)lbc

+ b(uB − ui

liB) cos(αB)lca

+ b(uC − ui

liC) cos(αC)lab = 0, (2.14)


(a) (b)

Figura 2.3: Ejemplo de localizacion de los baricentros de los volumenes de control (a)para triangulos agudos, (b) para triangulos obtusos.


A los volumenes de control j-esimos pertenecientes a T∂Ω se les aplica la condicion decontorno Dirichlet. Conocida la temperatura del contorno, usamos el punto de la arista(bc) mas proximo al baricentro, ası ∇u y ~n son colineales verificandose la ecuacion (4.5) aligual que en el caso del circuncentro siendo en este caso ljs1 la distancia mınima entre elcentro de gravedad y la arista frontera. A partir de la ecuacion (2.4), para los volumenesde control con una arista contorno se obtiene

(a∂uj


uC − uj

ljC) cos(αC)lab

+ b(uB − uj

ljB) cos(αB)lca − b

uj

ljs1lbc = −b

uext

ljs1lbc. (2.15)

mientras que para los volumenes de control con dos aristas contorno, la ecuacion utilizadaes

(a∂uj


uC − uj

ljC) cos(αC)lab

− (lbcljs1

+lcaljs2

)buj = −buext

ljs1lbc − b

uext

ljs2lca. (2.16)

Se forma ası un sistema EDO que se resuelve mediante la ecuacion (2.12). El uso delbaricentro ciertamente de de mayor complejidad puesto que requiere para cada triangulode la malla del calculo de los angulos entre el gradiente de la magnitud u y la normal acada arista. Sin embargo, este metodo ofrece la garantıa de que todos los nodos de calculose localicen en el interior de su volumen de control correspondiente y por tanto dentrodel dominio de calculo. Ası mismo, el uso del metodo Delaunay para la discretizacion del


dominio asegura alta calidad de los triangulos y por tanto los cosenos de los angulos αX

seran cercanos a la unidad.

En el siguiente apartado se describe como se estructura y codifica la informacion dela geometrıa en el lenguaje de programacion utilizado para la resolucion de problemasPDEs, y en la seccion 2.6 se comparan los resultados obtenidos usando los distintos FVMimplementados.

2.3. Estructura de la informacion de la geometrıa

Como se ha mencionado anteriormente, la PDEToolbox de Matlab es una herramien-ta eficaz para la resolucion de PDEs en geometrıas complicadas por el metodo de loselementos finitos. Para ello dispone de funciones capaces de interpretar la geometrıa deldominio introducida por el usuario a traves de un interfaz o codificado mediante ma-trices [18]. A continuacion se presenta la estructura de la malla generada por Matlaby que puede ser exportada a la lınea de comandos. Posteriormente se vera como ha deser adaptada esta estructura para tener un facil acceso para el metodo de los volumenesfinitos implementado en la seccion 2.2.

2.3.1. PDEToolbox de Matlab

La geometrıa del dominio de computo puede ser generada mediante la GUI de formaintuitiva, sin embargo, es preferible para geometrıas complejas el uso de las funcionespropias de la toolbox. El proceso de codificacion de la geometrıa es el siguiente: en pri-mer lugar se construye la denominada matriz de descripcion de la geometrıa (GD). Cadacolumna de esta matriz representa un objeto basico del modelo de geometrıa. Existencuatro tipos de objetos. En la fila 1 se especifica el tipo de objeto:

Para un disco, la fila uno contiene un 1, la segunda y tercera fila contienen lascoordenadas del centro del disco y la fila cuarta determina el valor del radio deldisco.

Para un polıgono solido, la primera fila contiene un 2, la segunda fila contiene elnumero N de segmentos, las siguientes N filas dan la abscisa de los puntos de iniciode los segmentos y las siguientes N filas son para la ordenada de estos mismospuntos de inicio.

Para un rectangulo, la fila uno se pone a 3 y el formato es a partir de aqui identicoal del polıgono.

Para una elipse solida, la fila uno se pone a 4, la segunda y tercera fila contienen


las coordenadas del centro de la elipse. Las filas 4 y 5 contienen los ejes mayores ymenores de la elipse. La fila 6 almacena el angulo de rotacion de la elipse.

En segundo lugar se genera la llamada matriz de descomposicion de la geometrıa mediantela orden decsg. Se le pasan como parametros la matriz de descripcion GD, la formula deconjuntos que determina mediante los operadores ’+’, ’-’ y ’*’ el conjunto definitivoque conforma el dominio. El tercer parametro es el name space matriz que asocia a cadacolumna de GD el nombre del conjunto utilizado en la formula de conjuntos.

En tercer lugar mediante la orden initmesh, la PDEToolbox ejecuta el mallado detriangulos asociado por medio del metodo de Delaunay. Esta malla es codificada mediantetres matrices llamadas P, T y E cuya definicion es la siguiente:

Matriz P: o matriz de puntos con dimension (2×m). Posee dos filas para las coor-denadas (x e y) de los m puntos que forman los vertices de los triangulos (nodosde calculo para el FEM).

Matriz T: o matriz de triangulos con dimension (4×n). En las tres primeras filas sedisponen los ındices de los vertices de los n triangulos que componen la malla. Elorden de clasificacion de los vertices es contrario a las agujas del reloj. La cuarta filaes util para dominios no homogeneos ya que indica, mediante un numero natural,el subconjunto al que pertenecen los n triangulos.

Matriz E: o matriz de aristas. Esta matriz permite interpretar las aristas de lafrontera del dominio sin embargo no es de utilidad para la posterior aplicacion delmetodo de los volumenes finitos en este proyecto fin de carrera.

La matriz de puntos P y la matriz de triangulos T de dimension (2 × m) y (4 × n),respectivamente, se muestran en las tablas siguientes. En ellas se presenta un ejemplode codificacion numerica de la malla de la figura 2.4 generada por la PDEToolbox. Esteejemplo es utilizado en la seccion 2.3.2 para ilustrar el nuevo codigo implementado.

Matriz P Matriz TColumna 1 2 3 4 5 6 ... Columna 1 2 3 4 5 6 7 8 ...

Vertice 1 1 1 2 5 2 9 9 4 ...x -1 -0.2 0.2 -1 -0.6 -0.8 ... Vertice 2 4 9 5 9 7 4 8 6 ...y 1 1 0.77 0.4 1 0.06 ... Vertice 3 9 5 7 7 3 8 7 8 ...

Subconjunto 1 1 1 1 1 1 1 1 ...

2.3.2. Nueva estructura de la informacion para FVM con cir-cuncentro

El haber desarrollado el metodo de los volumenes finitos permite conocer exactamenteque datos son los necesarios para implementarlo. Solo nos hace falta definir una estructura


Figura 2.4: Ejemplo de numeracion de los nodos (numeros pequenos) y de los triangulos(numeros grandes) de una malla generada por la PDEToolbox de Matlab.

para poder clasificar esos datos de forma ordenada con lo que se facilita el posterioracceso a ellos. Se opta por reconstruir la informacion de la malla mediante dos nuevasmatrices P1 y T1. La matriz P1 sera una copia de la matriz P si bien en circunstanciasexcepcionales podra ver alterada las coordenadas de algunos de sus puntos. En efecto,el uso del circuncentro como nodo central de los volumenes de control presenta ciertosproblemas cuando la calidad de los triangulos de la malla no es alta, ver apendice C.

La matriz P1 altera la posicion de un vertice, para los triangulos cuyos circuncentrosson problematicos, en un entorno cercano por lo que sus circuncentros cambian de posiciony no provocan indeterminacion en las ecuaciones (2.7), (2.9) y (2.10). Para esta nuevamatriz P1 se recalcula la matriz T1 asociada y se comprueba que no ocurra nuevamenteel mismo problema. Esta matriz T1 anade informacion a la matriz T determinada por laPDEToolbox y pasa a ser de dimensiones (19 × n). En ella se almacenan la mayorıa delos datos relativos a la geometrıa como son las longitudes de las aristas, distancias entrepuntos, etc. A continuacion se detallan los datos disponibles en esta matriz.

El triangulo j-esimo de la malla lleva asociado la columna j-esima de la matriz T1

con las caracterısticas que se muestran en la tabla siguiente.


Matriz T1Fila no Nomenclatura Descripcion

1, 2, 3 y 4 - Identicas a las de la matriz T generada por la PDEToolbox

5 Aj Area del triangulo (m2)6 xcircj Abcisa (m) del circuncentro del triangulo (j)7 ycircj Ordenada (m) del circuncentro del triangulo (j)

8, 9 y 10 (lab)j, (lbc)j y (lca)j Longitud (m) de las aristas de vertices a-b, b-c y c-a.11, 12 y 13 ljC , ljA y ljB Distancia (m) del circuncentro del triangulo (j) a los

circuncentros de los triangulos anexos14 N Indica el numero de aristas contorno. Para los triangulos

pertenecientes a TΩ toma el valor 0.

15, 16 y 17 TjC , TjA y TjB Indice de los triangulos anexos al triangulo (j). Tomanel valor 0 en caso de que no existan

18 y 19 ljs1 y ljs2 Longitudes mınimas (m) entre el circuncentro de (j) ylas aristas contorno en caso de ser triangulo frontera

El area del triangulo j-esimo con vertices abc es igual a Aj = 12det( xa ya 1; xb yb 1; xc

yc 1).

Las coordenadas del circuncentro del triangulo j-esimo, xcircj e ycircj, se calcu-lan mediante la funcion circuncentro con los siguientes parametros de entrada(xa, ya, xb, yb, xc, yc), se detalla en el siguiente apartado.

Las longitudes de los segmentos arista del triangulo j-esimo, (lab)j, (lbc)j y (lca)j,se calculan facilmente al conocer las coordenadas de los vertices, p.e. (lab)j =√

(xa − xb)2 + (ya − yb)2. Del mismo modo se calculan las distancias entre circun-centros, ljC , ljA y ljB.

La determinacion de las filas 15, 16 y 17 que recogen el ındice de los triangulosvecinos (en sentido antihorario) a un determinado triangulo j es de especial im-portancia para el desarrollo de la implementacion. Para el ejemplo de malla de laseccion anterior se expone el proceso llevado a cabo para determinar los ındices delos triangulos vecinos al triangulo 2. En primer lugar se hallan los triangulos cuyosvertices son compartidos con el triangulo 2. Entre estos triangulos se determinancuales tienen en comun dos vertices con el triangulo 2. Para la arista (1-9) se tienenlos triangulos 1 y 2 por tanto la fila 15 de la matriz T1 toma el valor 1 por ser esteel ındice del triangulo vecino. Recorriendo las aristas en sentido antihorario se tieneque 4 es vecino de 2. Para la ultima arista (5-1) solo se tiene como triangulo candi-dato el mismo triangulo 2, por tanto la arista pertenece a la frontera del dominio.Se deja a 0 la fila 17 de la matriz T1.

Las distancias mınimas entre los circuncentros de los triangulos frontera y el con-torno, ljs1 y ljs2, se calculan mediante la funcion distancia_minima que se detallaen el apendice B.1.

En el apartado 2.4 se muestra la implementacion del codigo necesario para obtenerla matrix T1.


2.3.3. Nueva estructura de la informacion para FVM con bari-centro

Para la implementacion del metodo de volumenes finitos de cuatro puntos con usodel baricentro se recurre en este caso a clasificar toda la informacion necesaria en unamatriz nombrada T2, de dimensiones (22×n), que contiene como suplemento a la matrizT1 varios datos de los angulos.

Matriz T2Fila no Nomenclatura Descripcion1 a 5 - Identicas a las de la matriz T1

6 xbarij Abcisa (m) del baricentro del triangulo (j)7 ybarij Ordenada (m) del baricentro del triangulo (j)

8 a 19 - Identicas a las de la matriz T120 cos(αC) Coseno del angulo que forman el flujo conductivo con la normal

... a la arista (ab)21 cos(αA) ... a la arista (bc)22 cos(αB) ... a la arista (ca)

2.4. Implementacion en Matlab

2.4.1. Codigos Matlab para FVM con circuncentro

La resolucion del modelo mediante FVM haciendo uso de la PDEToolbox consta detres etapas. En primer lugar se ha de generar la matriz T1, que contiene las propiedadesde los VC vistos anteriormente, mediante la funcion matrizT1.m. Posteriormente, seconstruyen las matrices A, B y C del sistema de la ecuacion (2.11) mediante la funcionsistemaABC.m. Por medio de una tercera funcion llamada Resolucion.m se halla lasolucion del sistema de ecuaciones diferenciales ordinarias anterior. A continuacion, sedetalla el codigo implementado para cada una de ellas.

Codigo para la construccion de la matriz T1

La construccion de la matriz T1 requiere del uso de las matrices P y T generadas porla PDEToolbox. Hacemos uso de la notacion vista previamente.

En primer lugar, para cada uno de los n triangulos se obtienen las coordenadas de susvertices, su area mediante el calculo de un determinante, las coordenadas del circuncentroy las longitudes de las aristas.


matrizT1.m

function T1=matrizT1(p,t);

1 n=length(t);T1=sparse(zeros(19,n));T1(1:4,:)=t; imagi=i;

2 for i=1:n

3 xa=p(1,t(1,i));

4 xb=p(1,t(2,i));

5 xc=p(1,t(3,i));

6 ya=p(2,t(1,i));

7 yb=p(2,t(2,i));

8 yc=p(2,t(3,i));

9 T1(5,i)=1/2*det([xa ya 1;xb yb 1;xc yc 1]);

10 [T1(6,i),T1(7,i)]=circuncentro(xa,ya,xb,yb,xc,yc);

11 T1(8,i)=sqrt((xa-xb)^2+(ya-yb)^2);

12 T1(9,i)=sqrt((xb-xc)^2+(yb-yc)^2);

13 T1(10,i)=sqrt((xc-xa)^2+(yc-ya)^2);

14 end

Sigue . . .

La funcion circuncentro calcula las coordenadas del circuncentro de un triangulodadas las coordenadas de sus vertices. Para ello calcula las ecuaciones de las mediatricesde dos segmentos, y = max+pa e y = mbx+pb, y se halla su interseccion, ver apendice B.2

A continuacion se hace una busqueda de los tres triangulos anexos a cada triangulo.En caso de ser frontera, el triangulo tendra unicamente 1 o 2 triangulos vecinos, la variableN indica el numero de de aristas frontera para estos triangulos (fila 14 de T1). Con estosdatos se pueden calcular las distancias entre los circuncentros de los triangulos vecinos yla distancia mınima al contorno para los triangulos frontera.

matrizT1.m (cont.)

15 for i=1:n

16 set_a=find(t(1,:)==t(1,i)|t(2,:)==t(1,i)|t(3,:)==t(1,i));

17 set_b=find(t(1,:)==t(2,i)|t(2,:)==t(2,i)|t(3,:)==t(2,i));

18 set_c=find(t(1,:)==t(3,i)|t(2,:)==t(3,i)|t(3,:)==t(3,i));

---------- Arista (ab) ----------

19 C=setxor(intersect(set_a,set_b),i);

20 if ~isempty(C)

21 T1(11,i)=sqrt((T1(6,i)-T1(6,C))^2+((T1(7,i)-T1(7,C))^2);

22 T1(15,i)=C;

Sigue . . .


matrizT1.m (cont.)

23 else

24 T1(14,i)=1;

25 T1(18,i)=distanciaminima(T1(6,i),T1(7,i),p(1,T1(1,i)),...

p(2,T1(1,i)),p(1,T1(2,i)),p(2,T1(2,i)));

26 end

---------- Arista (bc) ----------

27 A=setxor(intersect(set_b,set_c),i);

28 if ~isempty(A)

29 T1(12,i)=sqrt((T1(6,i)-T1(6,A))^2+(T1(7,i)-T1(7,A))^2);

30 T1(16,i)=A;

31 else

32 T1(14,i)=T1(14,i)+1;

33 if T1(14,i)==1


p(2,T1(2,i)),p(1,T1(3,i)),p(2,T1(3,i)));

35 elseif T1(14,i)==2;


p(2,T1(2,i)),p(1,T1(3,i)),p(2,T1(3,i)));

37 end

38 end

---------- Arista (ca) ----------

39 B=setxor(intersect(set_c,set_a),i);

40 if ~isempty(B)

41 T1(13,i)=sqrt((T1(6,i)-T1(6,B))^2+(T1(7,i)-T1(7,B))^2);

42 T1(17,i)=B;

43 else

44 T1(14,i)=T1(14,i)+1;

45 if T1(14,i)==1


p(2,T1(1,i)),p(1,T1(3,i)),p(2,T1(3,i)));

47 elseif T1(14,i)==2;


p(2,T1(1,i)),p(1,T1(3,i)),p(2,T1(3,i)));

49 end

50 end

51 end

Fin matrizT1.m

Para la determinacion de los indices de los triangulos vecinos al triangulo i se pro-cede de la siguiente manera. En primer lugar se hallan los conjuntos de triangulos quetienen alguno de los tres vertices de i en comun (set_a, set_b, set_c). A continuacion


para la arista (ab) se halla la interseccion del conjunto de triangulos que tienen a y bpor vertice, excluyendo de ellos al propio triangulo i. Esta operacion se hace median-te C=setxor(intersect(set_a,set_b),i);, siendo C el ındice del triangulo vecino a imediante la arista (ab). Se realiza de la misma forma para las aristas (bc) y (ca).

Otra opcion es recorrer los triangulos uno por uno y comparar si se corresponde algunapareja de sus vertices con alguna pareja de vertices del triangulo i (equivale a tener unaarista en comun); y realizar este proceso para todos los triangulos i de la malla. Estemetodo tiene una complejidad de orden (n2) mientras que la implementacion llevada acabo es de orden (n) con el correspondiente ahorro de tiempo que conlleva al hacer usode mallas suficientemente finas y mejorar los resultados obtenidos.

Codigo para la construccion del sistema EDO (ec. (2.11))

La funcion sistemaABC.m devuelve las matrices masa y rigidez, A y B, respectiva-mente, y el vector independiente C del sistema Adu + Bu = C.

sistemaABC.m

function [A,B,C]=sistemaABC(T1,a,b,c,d,uext);

n=length(T1) A=spalloc(n,n,160000);

B=spalloc(n,n,160000); C=sparse(zeros(n,1));

1 for i=1:n

2 if T1(14,i)==0

Para los triangulos interiores se hace uso de la ecuacion (2.7)

3 B(i,T1(15,i))=b*T1(8,i)/T1(11,i);

4 B(i,T1(16,i))=b*T1(9,i)/T1(12,i);

5 B(i,T1(17,i))=b*T1(10,i)/T1(13,i);

6 B(i,i)=c*T1(5,i)-b*(T1(8,i)/T1(11,i)+T1(9,i)/T1(12,i)+...

T1(10,i)/T1(13,i));

7 C(i,1)=d*T1(5,i);

Para los triangulos con una arista contorno se hace uso de la ecuacion (2.9)

8 elseif T1(14,i)==1

9 ljs1=T1(18,i)

10 if T1(11,i)==0

11 C(i)=d*T1(5,1)-uext*b*T1(8,i)/ljs1;

12 B(i,T1(16,i))=b*T1(9,i)/T1(12,i)+c*T1(5,i);

13 B(i,T1(17,i))=b*T1(10,i)/T1(13,i)+c*T1(5,i);


Sigue . . .


sistemaABC.m (cont.)

T1(8,i)/ljs1);

15 elseif T1(12,i)==0

16 C(i)=d*T1(5,1)-uext*b*T1(9,i)/ljs1;

17 B(i,T1(15,i))=b*T1(8,i)/T1(11,i)+c*T1(5,i);

18 B(i,T1(17,i))=b*T1(10,i)/T1(13,i)+c*T1(5,i);


T1(9,i)/ljs1);

20 elseif T1(13,i)==0

21 C(i)=d*T1(5,1)-uext*b*T1(10,i)/ljs1;

22 B(i,T1(15,i))=b*T1(8,i)/T1(11,i)+c*T1(5,i);

23 B(i,T1(16,i))=b*T1(9,i)/T1(12,i)+c*T1(5,i);


T1(10,i)/ljs1);

25 end

Para los triangulos con dos aristas contorno se hace uso de la ecuacion (2.10)

26 elseif T1(14,i)==2

27 ljs1=T1(18,i)

28 ljs2=T1(19,i)

29 if T1(11,i)~=0

30 C(i)=d*T1(5,1)-uext*b*T1(9,i)/ljs1-uext*b*T1(10,i)/ljs2;

31 B(i,T1(15,i))=b*T1(8,i)/T1(11,i);

32 B(i,i)=c*T1(5,i)-b*(T1(9,i)/ljs1+T1(10,i)/ljs2+...

T1(8,i)/T1(11,i));

33 elseif T1(12,i)~=0


35 B(i,T1(16,i))=b*T1(9,i)/T1(12,i);


T1(9,i)/T1(12,i));

37 elseif T1(13,i)~=0

38 C(i)=d*T1(5,1)-uext*b*t1(8,i)/ljs1-uext*b*T1(9,i)/ljs2;

39 B(i,T1(17,i))=b*T1(10,i)/T1(13,i);


T1(10,i)/T1(13,i));

41 end

42 end

Matriz A de coeficientes, indistintamente de la localizacion del triangulo.

44 A(i,i)=a*T1(5,i);

45 end;

Fin sistemaABC.m


Codigo para la resolucion del sistema EDO

En segundo lugar llevamos a cabo la resolucion del sistema de ecuaciones diferencialesordinarias construido mediante la implementacion de la ecuacion (2.12).

Resolucion.m

function sol=Resolucion(A,B,C,tlist,uini);

1 u0=sparse(zeros(length(C),1)+uini);

2 n=length(tlist) sol=zeros(length(C),n); aux=u0;

3 sol(:,1)=u0;

4 for i=2:n

5 delt=tlist(i)-tlist(i-1);

6 sol(:,i)=(A/delt+B/2)\(C+(A/delt-B/2)*aux);

7 aux=sol(:,i);

9 end

Fin Resolucion.m

2.4.2. Codigos Matlab para FVM con baricentro

A continuacion se muestra el codigo modificado para implementar el FVM usando elbaricentro. En primer lugar construimos la matriz T2 mediante la funcion matrizT2.m.En negrita se muestran las lıneas donde se han realizado modificaciones respecto delprograma que hace uso del circuncentro.

matrizT2.m

function T2=matrizT2(p,t);

1 n=length(t);T1=sparse(zeros(22,n));T1(1:4,:)=t; imagi=i;

2 for i=1:n

3 xa=p(1,t(1,i));

4 xb=p(1,t(2,i));

5 xc=p(1,t(3,i));

6 ya=p(2,t(1,i));

7 yb=p(2,t(2,i));

8 yc=p(2,t(3,i));

9 T2(5,i)=1/2*det([xa ya 1;xb yb 1;xc yc 1]);

10 T2(6,i)=(xa+xb+xc)/3;

Sigue . . .


matrizT2.m (cont.)

11 T2(7,i)=(ya+yb+yc)/3;

12 T2(8,i)=sqrt((xa-xb)^2+(ya-yb)^2);

13 T2(9,i)=sqrt((xb-xc)^2+(yb-yc)^2);

14 T2(10,i)=sqrt((xc-xa)^2+(yc-ya)^2);

15 end

16 for i=1:n

17 set_a=find(t(1,:)==t(1,i)|t(2,:)==t(1,i)|t(3,:)==t(1,i));

18 set_b=find(t(1,:)==t(2,i)|t(2,:)==t(2,i)|t(3,:)==t(2,i));

19 set_c=find(t(1,:)==t(3,i)|t(2,:)==t(3,i)|t(3,:)==t(3,i));

Sigue . . .

A continuacion, para la arista (ab) se hallan todas las propiedades que la afectan talescomo longitud entre el baricentro i y el baricentro C del triangulo vecino, el coseno delangulo que forma el gradiente de la propiedad u con la normal, etc.

Para el calculo del coseno del angulo α se procede como sigue. En primer lugar,hallamos las coordenadas del vector normal, vec_n, externo a la arista. Posteriormentese calculan las coordenadas del vector flujo conductivo, vec_iC y se obtiene el coseno delangulo, cos_alfa_C como el coseno de la diferencia de los angulos de cada vector. Al sercos(α) = cos(−α), no hay problema con el signo del angulo.

matrizT2.m (cont.)

20 C=setxor(intersect(set_a,set_b),i);

21 if ~isempty(C)

22 T2(11,i)=sqrt((T2(6,i)-T2(6,C))^2+(T2(7,i)-T2(7,C))^2);

23 vec_n=p(2,t(2,i))-p(2,t(1,i))-(p(1,t(2,i))-p(1,t(1,i)))*imagi;

24 vec_iC=T2(6,C)-T2(6,i)+(T2(7,C)-T2(7,i))*imagi;

25 cos_alfa_C=cos(angle(vec_iC)-angle(vec_n));

26 T2(15,i)=C;

27 else

28 T2(14,i)=1;


p(2,T2(1,i)),p(1,T2(2,i)),p(2,T2(2,i)));

30 T2(11,i)=0; cos_alfa_C=1;

31 end

32 T2(20,i)=cos_alfa_C;

Fin matrizT2.m

De forma analoga se realiza para las dos aristas restantes (bc) y (ca), ver apendice ??,


siguiendo el desarrollo mostrado para la matriz T1 de circuncentros.

La funcion sistemaABC2.m devuelve las matrices masa y rigidez, A y B, respectiva-mente, y el vector independiente C del sistema Adu + Bu = C. Mientras que la funcionResolucion.m es valido para resolver el sistema creado para este FVM mediante bari-centro.

sistemaABC2.m

sistemaABC2.m function [A,B,C]=sistemaABC2(T2,a,b,c,d,uext);



1 for i=1:length(t1)

2 if t1(14,i)==0

Triangulos interiores, ecuacion (2.14)3 B(i,t1(15,i))=b*T2(8,i)*T2(20,i)/T2(11,i);

4 B(i,t1(16,i))=b*T2(9,i)*T2(21,i)/T2(12,i);

5 B(i,t1(17,i))=b*T2(10,i)*T2(22,i)/T2(13,i);

6 B(i,i)=c*T2(5,i)-b*(T2(8,i)*T2(20,i)/T2(11,i)+T2(9,i)*...

T2(21,i)/T2(12,i)+T2(10,i)*T2(22,i)/T2(13,i));

7 C(i,1)=d*T2(5,i);

Triangulos con una arista perteneciente al contorno, ecuacion (2.15)


9 ljs1=T2(18,i)

10 if T2(11,i)==0

11 C(i)=d*T2(5,1)-uext*b*T2(8,i)/ljs1;

12 B(i,T2(16,i))=b*T2(9,i)*T2(21,i)/T2(12,i)+c*T2(5,i);

13 B(i,T2(17,i))=b*T2(10,i)*T2(22,i)/T2(13,i)+c*T2(5,i);

14 B(i,i)=c*T2(5,i)-b*(T2(9,i)*T2(21,i)/T2(12,i)+...

T2(10,i)*T2(22,i)/T2(13,i)+T2(8,i)/ljs1);

15 elseif T2(12,i)==0

16 C(i)=d*T2(5,1)-uext*b*T2(9,i)/ljs1;

17 B(i,T2(15,i))=b*T2(8,i)*T2(20,i)/T2(11,i)+c*T2(5,i);

18 B(i,T2(17,i))=b*T2(10,i)*T2(22,i)/T2(13,i)+c*T2(5,i);

19 B(i,i)=c*T2(5,i)-b*(T2(8,i)*T2(20,i)/T2(11,i)+...

T2(10,i)*T2(22,i)/T2(13,i)+T2(9,i)/ljs1);

20 elseif t1(13,i)==0

21 C(i)=d*T2(5,1)-uext*b*T2(10,i)/ljs1;

22 B(i,T2(15,i))=b*T2(8,i)*T2(20,i)/T2(11,i)+c*T2(5,i);

Sigue . . .


sistemaABC2.m (cont.)

23 B(i,T2(16,i))=b*T2(9,i)*T2(21,i)/T2(12,i)+c*T2(5,i);

24 B(i,i)=c*T2(5,i)-b*(T2(8,i)*T2(20,i)/T2(11,i)+...

T2(9,i)*T2(21,i)/T2(12,i)+T2(10,i)/ljs1);

25 end

Fin sistemaABC2.m

El codigo se completa con la implementacion de la ecuacion (2.16) correspondientea los triangulos con dos aristas perteneciente al contorno, asociando cada termino a lasmatrices A, B y C. El codigo completo, muy similar al del caso FVM con circuncentro,se adjunta en el CD-Rom.

2.5. Validacion del metodo

Se procede a comprobar la solucion obtenida por el esquema FVM comparandolocon los resultados de un metodo numerico de elementos finitos generado por un softwarecomercial sobre la misma malla y por la solucion analıtica en un problema lineal basico.El objetivo es medir el error de dos metodos numericos (uno implementado por el autordel proyecto y otro implementado por un software comercial) con respecto a la solucionexacta obtenida para un problema que pueda resolverse analıticamente. Si estos erroresentran dentro de unos lımites aceptables se puede considerar que el metodo desarrolladopara resolver problemas regidos por PDEs lineales esta validado y a partir de aquı seprocedera a desarrollar un metodo tambien basado en FVM que resuelva problemas nolineales y por tanto mucho mas complejos.

El problema en cuestion es la ecuacion en derivadas parciales parabolica siguiente

∂u

∂t−∆u = 0, t > 0, (2.17)

que se resuelve en un dominio plano circular de radio R, con condicion de contorno linealde tipo Dirichlet u = uext, es decir valor constante de la magnitud u en la superficie, ycon condiciones iniciales nulas, u0 = 0.

Solucion analıtica

En este apartado se recoge la solucion del problema mencionado y obtenida de formaanalıtica mediante un metodo de transformaciones de Laplace [19]. Para ello reescribimosla ecuacion a coordenadas cilındricas quedando el problema siguiente

∂2u

∂r2+

1

r

∂u

∂r− ∂u

∂t= 0, 0 ≤ r < R, t > 0, (2.18)


con u(r = R) = uext.

Desarrollando el metodo [19], se llega a la solucion analıtica siguiente

u = uext −2uext

R

∞∑n=1

exp

(−α2

ntJ0(rαn)

αnJ1(Rαn)

)(2.19)

donde αn son las raıces de J0(Rα) = 0, J0 y J1 son las funciones de Bessel de primeraespecie de orden 0 y 1, respectivamente.

La solucion analıtica, hallada en este apartado, sirve para calcular la respuesta tempo-ral de la magnitud u en un solido macizo plano con condiciones de contorno constantes.Con la ayuda del resultado obtenido se pueden ajustar distintos parametros para losmetodos numericos de la PDEToolbox y de los volumenes finitos, expuestos en apartadosanteriores, y comprobar los errores cometidos por ellos.

Solucion mediante elementos finitos

El metodo de los elementos finitos (FEM) es un metodo de calculo utilizado en di-versos problemas de ingenierıa, que se basa en considerar al cuerpo o estructura divididoen elementos discretos, con determinadas condiciones de vınculo entre sı, generandose unsistema de ecuaciones que se resuelve numericamente. Se suele utilizar en matematicascomo metodo nodal aproximado para resolver ecuaciones diferenciales y PDEs en formanumerica. Existen dos tipos de caminos para su formulacion, basandose en el principio delos trabajos virtuales, es decir, formulaciones variacionales, o mediante el metodo de Gar-lekin. El metodo es muy poderoso debido a su generalidad y a la facilidad de introducirdominios de calculo complejos (en dos o tres dimensiones). Es muy utilizado en el calculode deformaciones (mecanica del continuo), campos de velocidades y presiones (fluido-dinamica) y transferencia de calor. Varios de estos problemas no tienen solucion analıticao es muy difıcil obtenerla, por lo que se convierte en la unica alternativa de resolucion.Para la resolucion de problema descrito anteriormente se hace uso de la PDEToolbox deMatlab, que integra un metodo de elementos finitos.

Comparacion de los resultados por los tres metodos

Con el fin de verificar la validez del metodo implementado se procede a la resoluciondel problema lineal planteado anteriormente. Se realizan las simulaciones para el dominiogeometrico establecido con una malla de 1058 triangulos y con 562 vertices. Por tantoel metodo numerico de los volumenes finitos tiene 1058 incognitas mientras que el de losvolumenes finitos tiene solo 562 y por tanto en este aspecto una menor precision. Losresultados se comparan con los de la solucion analıtica (que es posible determinar parael problema elegido) hallada mediante transformadas de Laplace. En la figura 2.5.a se


(a) (b)

Figura 2.5: Errores absolutos de la variable u obtenidos por los metodos FEM resueltomediante software comercial (puntos), FVM con circuncentros (+) y FVM con baricentro(x), en comparacion con la solucion analıtica, para nodos situados a una distancia radial(a) r=0, y (b) r=0.5.

observan los errores absolutos de la magnitud u cometidos por los distintos metodos enun nodo localizado en el centro del dominio circular y la figura 2.5.b recoge los erroresen otro nodo localizado a una distancia radial r = 0,5. En ambos casos los erroresmaximos absolutos estan limitados por debajo de 10−2 para una magnitud que varia entre0 y 1. El comportamiento del error en el FVM circuncentrico es analogo al del metodode elementos finitos, desarrollado mediante software comercial, obteniendose en amboserrores totalmente tolerables. No obstante, sorprende que el error del FVM baricentricodifiera tanto de los otros dos metodos a pesar de ser admisible. En todo caso, los resultadosmuestran inequıvocamente que el codigo desarrollado esta a la altura de uno comercial yque genera buenos resultados y por tanto se puede considerar validado sobre todo en elcaso del FVM circuncentrico. En el apartado siguiente se hace un analisis de los resultadosobtenidos por las dos variantes del FVM.

2.6. FVM baricentrico frente a FVM circuncentrico

en problemas lineales generales

Se esta interesado en comprobar los resultados generados por las dos variantes delmetodo de volumenes finitos desarrollados. Para ello se procede a resolver en un dominiocircular Ω, de un metro de radio, la ecuacion en derivadas parciales siguiente, mediante losdos metodos de volumenes finitos de cuatro puntos desarrollados y compararlo medianteFEM comercial,

100∂U

∂t−∇ · (∇U) + U = 10 , t > 0, ~x ∈ Ω ⊂ IR2, (2.20)


(a) (b) (c)

Figura 2.6: (a) Localizacion de los circuncentros (+) y de los baricentros (x) para unaseccion de malla obtenida por triangularizacion de Delaunay. (a) Respuesta transitoria dela variable u usando FVM circuncentrico para los nodos C1, C2, un nodo en r=0.5 y unnodo en r=0, representados mediante lıneas continua, a trazos, punteada y trazo-punto,respectivamente. (b) Respuesta transitoria de la variable u usando FVM baricentrico paralos nodos B1, B2, un nodo en r=0.5 y un nodo en r=0, representados mediante lıneascontinua, a trazos, punteada y trazo-punto, respectivamente.

con condicion de contorno Dirichlet U(~x) = 1 para ~x ∈ ∂Ω, y condicion inicial U(t =0) = 0.

La malla utilizada es generada por el codigo Delaunay implementado por Matlab conla restriccion de que los triangulos han de tener longitudes de aristas inferiores a 0.1 m.El dominio se discretiza en 1058 triangulos de alta calidad con un total de 562 verticesdiferentes. Se resuelve el problema mediante FEM con 562 nodos, y FVM con uso de cir-cuncentro y baricentro con 1058 nodos para ambos. Apreciamos en la figura 2.6.a comoa pesar de la buena calidad global de la malla existen triangulos comprometedores dondeel circuncentro no queda recluido en su interior (p.e. circuncentro C2) y se localiza en elVC vecino o incluso fuera del dominio de calculo (p.e. circuncentro C1) con los erroresque ello origina. Estos casos se contabilizan en 24 por tanto inferior a un 2,3% del totalde los triangulos presentan circuncentros problematicos de los cuales unicamente un cir-cuncentro esta fuera del dominio y perturba realmente el metodo. Mientras, logicamente,todos los baricentros quedan recluidos en el interior de sus volumenes de control.

Las figuras 2.6.b y 2.6.c muestran la respuesta transitoria de la variable u obtenida porel FVM con circuncentro y baricentro, respectivamente, para los nodos de varios VC dela malla. Los puntos cercanos a la superficie son los que evolucionan de forma mas rapiday alcanzan en regimen permanente valores cercanos a los de la condicion de contorno.Por otro lado, las zonas centrales del dominio tardan mas tiempo en evolucionar y es enellas en las que se alcanzan valores mayores para la variable u en regimen permanente.Se comprueba que los resultados de los dos metodos son similares teniendose un errormaximo entre los dos metodos para el nodo central inferior al 0.15%. No obstante, asimple vista se observa para FVM con circuncentro de la existencia de un resultado


Figura 2.7: Grafica 3D del valor de u en el instante de tiempo t = 10, resuelto medianteFVM circuncentrica y plano π de proyeccion.

erroneo con valor oscilante de gran amplitud. Este es el valor que toma la magnitud u enel nodo que se localiza fuera del dominio de calculo antes mencionado, es decir en C1.Por suerte este error es local, solo afecta a un pequeno numero de nodos vecinos paralos que se atenua rapidamente y puede controlarse mediante ajuste del paso de tiempo ydel parametro de malla. El FVM con baricentro es ajeno a estos problemas y muestra elcomportamiento esperado para los resultados.

Los resultados obtenidos tambien se pueden representar mediante graficas 3D, ver fi-gura 2.7. En ella se muestra la solucion obtenida por FVM con circuncentro, en el instantede tiempo t = 10, para todos los nodos de calculo. La solucion para el circuncentro C1, seve que no es del todo coherente con el resto de datos, esto es debido a los problemas yacomentados del circuncentro. En esta figura se muestra tambien un plano π sobre el quese van a proyectar los resultados de cada nodo obtenidos por las dos variantes de FVMy por FEM, para hacer un analisis mas claro y detallado.

Los resultados proyectados sobre el plano π, ver figura 2.8.a, representan el perfilradial de temperaturas en varios instantes de tiempo. Es de especial interes compararlos resultados de los metodos implementados con los de un software comercial como esla PDEToolbox de Matlab. En la figura 2.8.b se muestra la solucion obtenida por FEM(PDEToolbox), FVM con baricentro y FVM con circuncentro para cada nodo, proyecta-da en el plano π, representada mediante puntos negros, cruces azules y signos mas rojos,respectivamente. Se observa un aspecto inesperado puesto que a pesar de tener que existir


(a) (b)

Figura 2.8: (a) Perfil radial de la magnitud u proyectada en el plano π en distintosinstantes de tiempo. (b) Valor de la magnitud u proyectada en el plano π obtenida porlos metodos FEM (puntos negros), FVM con circuncentro (signo mas rojo) y FVM conbaricentro (cruces azules), en sucesivos pasos de tiempo.

simetrıa axial en los resultado, por ser el problema simetrico tanto en ecuaciones comoen dominio de calculo, el FVM con baricentro presenta cierta dispersion en los valoresque toma la magnitud u para distancias radiales muy proximas. En comparacion el FVMcon circuncentro y el FEM presentan de forma correcta la simetrıa al no exhibir disper-sion y sus resultados son muy proximos. Por tanto a pesar de los problemas que puedepresentar un metodo de volumenes finitos con circuncentro, los resultados obtenidos sonmas coherentes que los obtenidos mediante el uso del baricentro en mallados triangularesno estructurados.

Para ver la influencia del mallado en los resultados se representa el error relativo enregimen permanente de la solucion obtenida mediante las dos variantes de FVM en mallaDelaunay y en malla simetrica, respectivamente, ver figuras 2.9. El error se ha tomadorespecto a los resultados del FEM. El uso del circuncentro asegura una buena solucional ser el error admisible para las dos mallas utilizadas. No obstante, se ha de recordar laposibilidad de que surjan nodos complicados, como el anteriormente mencionado C1 quepueden generar soluciones locales oscilatorias con errores importantes. Los resultadosdel baricentro son mas complicados de analizar. En malla de Delaunay el error mediorelativo es admisible, no obstante un analisis de la desvicion tıpica de este error muestraque existe una dispersion en torno a la solucion exacta que si bien no es exagerada si quea priori no era de esperar. En malla simetrica, el resultado obtenido con baricentro es eneste caso totalmente inaceptable. El error medio es mucho mas elevado que el obtenidopor el circuncentro al igual que la dispersion. Por tanto, el efecto del mallado es muyimportante en el baricentro para el cual los resultados no son del todo fiables mientrasque el circuncentro se muestra mas insensible a la malla utilizada y sus resultados soncoherentes con los obtenidos mediante FEM.


(a) (b)

Figura 2.9: (a) Error relativo en regimen permanente de los FVM con baricentro (crucesazules) y circuncentro (signo mas rojo) en malla Delaunay. (b) Error relativo en regimenpermanente de los FVM con baricentro (cruces azules) y circuncentro (signo mas rojo)en malla Simetrica.

Por tanto el tipo de malla utilizada influye decisivamente en FVM con baricentro. En lafigura 2.10 se muestran la malla generada por el algoritmo de Delaunay y la malla llamadasimetrica generada mediante refinamientos sucesivos, sobre las que se van a realizar loscalculos. La malla de Delaunay presenta como ventajas el tener una calidad media detriangulos mayor mientras que la malla simetrica tiene una la calidad de triangulo mınimamuy superior. Para mas informacion sobre la calidad de los triangulos, q, ir al apendice C.El valor de la calidad de varias mallas es mostrado en la tabla siguiente, cuanto masproximo a la unidad sea q, mejor es la malla. Por tanto se tiene que la mallas de Delaunayson de mejor calidad pero existen ciertos triangulos de mala calidad. A continuacion vemoscomo afectan estas mallas a los dos FVM implementados.

FVM Error relativo respecto FEM Calidad q mallaTipo de malla No triangulos Circun. Bari. Max (%) Medio ( %) STD Media Mın

Delaunay X 2.58 1.73 0.19 0.78 0.95(h=0.3) 120 X 8.13 1.47 0.23

Delaunay X 0.99 0.44 0.10 0.73 0.96(h=0.15) 490 X 3.56 0.44 0.14Delaunay X 0.46 0.13 0.07 0.77 0.96(h=0.07) 2184 X 2.03 0.03 0.18Simetrica X 2.56 1.85 0.21 0.87 0.90

(3 refinados) 128 X 9.74 6.70 0.8Simetrica X 0.71 0.47 0.10 0.86 0.90

(4 refinados) 512 X 7.45 4.86 1.21Simetrica X 0.16 0.07 0.38 0.85 0.90

(5 refinados) 2048 X 7.64 4.40 2.22

Un analisis mas profundo del efecto del mallado se realiza con la ayuda de los datos


(a) (b)

Figura 2.10: (a) Ejemplo de malla de Delaunay. (b) Ejemplo de malla simetrica obtenidapor refinamiento sucesivo.

de la tabla anterior. Para mallas de tipo Delaunay tanto el error relativo medio como elerror relativo maximo del FVM con circuncentro, en regimen permanente disminuye deforma sensible con el aumento del numero de nodos o lo que es lo mismo al disminuirel parametro de malla h. En mallas simetricas el comportamiento es analogo siendo losresultados totalmente admisibles llegando a ser levemente mejores en mallas simetricasque en las de Delaunay. Esto sin duda es debido a que a pesar de ser la calidad media delos triangulos mejor en mallas Delaunay la calidad de algunos triangulos puede llegar aser mala lo que distorsiona muy localmente la solucion en ciertos nodos.

Cuantos mas puntos tengan las mallas de Delaunay, mejor comportamiento tieneFVM con baricentro. Este comportamiento no es tan claro en mallas simetricas, en dondelos errores cometidos no son tolerables. El uso de baricentro y malla Delaunay permiteobtener errores relativos medio del 0.03%, siendo este valor ficticio puesto que existeuna amplia dispersion de estos errores como muestran los elevados valores que toma ladesviacion tıpica. Este comportamiento no era esperado y sorprende, la alta dependenciadel tipo de malla utilizada puede dar alguna pista de donde surgen los problemas.

Se ha resuelto una PDE parabolica que rige un problema general de calor median-te FEM y dos variantes FVM sobre una malla de Delaunay. No se hace uso de malladual (teselacion de Voronoi) y por tanto reducimos el numero de operaciones a realizar.Globalmente los resultados obtenidos por los tres metodos son coherentes y muy simila-res, sin embargo se tienen pequenas problematicas locales para los tres metodos. FEMpresenta pequenas oscilaciones en los pasos iniciales de calculo que provoca que ciertosnodos superficiales tengan soluciones fuera de la realidad fısica. El metodo FVM por serconservativo evita estos problemas no obstante sus dos variantes sobre la malla originaltambien generan resultados locales erroneos. FVM con uso del circuncentro como nodo delos VC provoca que estos se situen para los triangulos obtusos fuera de sus VC pudiendoseproducir nodos de calculo fuera del dominio de computo con resultados incorrectos. Estos


casos son excepcionales y se pueden corregir reajustando ciertos nodos de la malla origi-nal. El problema que presenta FVM con uso de nodos en los baricentros no es tan obvio.Los nodos siempre quedan recluidos en el interior de los VC, no obstante, la resoluciondel problema del calor en un dominio cilındrico no genera resultados simetricos respectoal eje existiendo una pequena dispersion. Este resultado inesperado no sucede en los otrosdos metodos desarrollados.

Obviamente este proyecto no trata de resolver la generalidad de todas las PDEs si-no que muestra distintas vıas numericas para la resolver problemas que se presentan amenudo en la ensenanza cientıfica. Con esta vision se ha resuelto un problema basico deingenierıa termica mediante tres metodos numericos y mostrado algunas ventajas e incon-venientes para cada metodo. A partir de estos resultados se puede proceder a desarrollarun metodo que resuelva PDEs no lineales.

Capıtulo 3

Problema no lineal delcalentamiento de preformas

Las fibras de cristales fotonicos (PCF) tambien llamadas fibras microestructuradas ofibras de huecos, son fibras en las cuales huecos de aire las recorren longitudinalmente. Lageometrıa de la distribucion de los huecos permite el control de las principales propiedadesopticas como son la dispersion [20], la birrefringencia [21] o las no linealidades [22] graciasal mecanismo de reflexion total interna. Otro tipo de fibra de huecos optica es la fibra debanda fotonica prohibida (photonic bandgap) [23], la cual conduce la luz haciendo usode las bandas prohibidas fotonicas que se generan en las estructuras periodicas.

Desde 2001, se esta fabricando un nuevo tipo de PCF, llamada fibra optica micro-estructurada de polımero (MPOF) [24, 25, 26]. Las temperaturas relativamente bajasde estirado de los polımeros, tales como el polimetilmetacrilato (PMMA), permiten alas MPOF ser una alternativa a las PCF de sılice para aplicaciones especıficas [27]. Lacalidad de los resultados del proceso de estirado es gobernado por la dependencia dela viscosidad y la tension superficial de la fibra con la temperatura. A temperaturas deestirado, la viscosidad del PMMA y del sılice son ambas del mismo orden de magnitud,mientras que la tension superficial del PMMA es un orden de magnitud inferior al delsılice [28, 29, 30]. Por tanto, la distorsion de la distribucion espacial de huecos provocadapor las tensiones superficiales pueden ser minimizadas reduciendo la temperatura de esti-rado, aun provocando un aumento de la viscosidad, permitiendo la fabricacion de MPOFde diametro micrometrico.

La fabricacion de MPOF suele realizarse en varias etapas de estirado de la preformahasta que se obtienen dimensiones de fibra, como se muestra en la figura 3.1. La estruc-tura de huecos de la preforma primaria puede ser obtenida por distintos procedimientoscomo apilado de capilares, extrusion, polimerizacion en molde, perforacion o fundiendoe inyectando en un molde. El proceso mas usado es perforar el patron de huecos en uncilindro de PMMA usando un taladrador controlado numericamente. La preforma prima-

34

CAPITULO 3: FORMULACION DEL PROBLEMA 35

Figura 3.1: Principales etapas de estirado en el proceso de fabricacion de fibras de cristalesphotonicos.

ria se estira obteniendose una preforma secundaria con la misma estructura de huecos.La preforma secundaria se introduce en un tubo de refuerzo de PMMA y es introducidaverticalmente en el horno. Tras calentarse, la preforma forma un cuello al estar sometidaa tensiones axiales de estirado obteniendose fibra de dimensiones micrometrica que saledel horno a una velocidad de estirado especificada. Fuera del horno, por debajo de este,un sistema de medicion de diametro por laser es usado como entrada para un sistema decontrol automatico para la posicion horizontal y la velocidad de estirado de la preforma.La fibra obtenida se enrolla en una bobina.

Los polımeros se caracterizan por tener una temperatura de transicion vıtrea, Tg,que determina sus propiedades mecanicas. El estirado de la preforma no es posible paratemperaturas por debajo de Tg por comportarse los polımeros como materiales fragiles ycomo resultado la fibra puede romperse. Para temperaturas por encima de Tg, el polımeroes un material plastico y por tanto las fuerzas a las que se somete durante el estiradoprovocan deformaciones permanentes. La temperatura de transicion vıtrea del PMMAdepende del tipo de polimerizacion llevada a cabo durante su preparacion, p.e., la sa-turacion y la estabilidad de las cadenas de polımero resultantes. Las varas de PMMAobtenidas por extrusion, pueden caracterizarse por una Tg de 95oC o 115oC segun que lapolimerizacion haya sido desasistida o asistida [27], respectivamente.

El proceso de calentamiento en el horno debe asegurar que la totalidad de la preformaalcance una temperatura superior a Tg, antes de comenzar el proceso de estirado, evitandoası o mas precisamente reduciendo el numero de colapsos de huecos y su distorsion,preservando de esta forma el patron de huecos disenado en el origen. El analisis dela evolucion de la temperatura en la preforma durante el calentamiento en el hornoes necesario para entender el comportamiento mecanico del material y las condicionesoperativas durante el estirado.

La estructura de huecos de aire del polımero influye en la distribucion de la temperatu-ra durante el calentamiento, introduciendo tensiones superficiales que pueden distorsionarla estructura resultante, especialmente para las fibras de banda prohibida fotonica [31].El impacto de las propiedades de los materiales y de las condiciones de estirado sobre el


tamano y la forma de los huecos de aire ha sido modelado para el transcurso del procesode estirado [32].

Las propiedades opticas de las PCF dependen enormemente del tamano y localizacionde los huecos de aire, siendo de vital importancia determinar los efectos de las condicionesde estirado en la seccion transversal de la fibra. Para ello, se han realizado numerososmodelos y estudios teoricos para fibras con un unico hueco asumiendo flujos isotermos.Muchos de estos modelos emplean una viscosidad constante y consideran fluidos Newto-nianos tambien con una viscosidad dinamica constante [30, 33, 34, 35] o una viscosidadque sigue una ley de tipo Arrhenius [36]. Reeve et al. [37], estudiaron el calentamientode preformas de fibras por medio de una formulacion bidimensional para los gases en elhorno teniendo en cuenta los intercambios de calor radiantes entre la preforma, el hornoy los irises. Ramos [38], estudio la fluido dinamica y los procesos de transferencia de calorque ocurren en la fabricacion de una fibra de material compuesto con un unico hueco so-metida a intercambios de calor de tipo convectivo y radiante entre las paredes del hornoy el conjunto preforma y fibra, por medio de un modelo acoplado.

Este estudio se centra en la transferencia de calor para preformas de fibras microes-tructuradas de material PMMA. La dependencia de las propiedades termicas del aire querellena los huecos y del polımero con la temperatura son tenidas en cuenta. Se trata demejorar el modelo de Lyytikaınen et al [39], y proponer un metodo numerico basado envolumenes finitos, que permita mejores aproximaciones a los resultados experimentalesobtenidos para PCF con disposiciones de huecos en anillos hexagonales. El modelo vienedeterminado por la ecuacion de conduccion del calor en un medio no homogeneo, con co-eficientes variables dependientes de la temperatura y con condicion de contorno no linealdebido a la radiacion.

Este capıtulo se estructura de la siguiente forma. En el apartado 3.1 se recopila lainformacion sobre las propiedades termicas del aire y se hace un estudio pormenorizadode los datos publicados en artıculos en los que se han intentado evaluar las caracterısticastermicas del PMMA. En el apartado 3.2 se describe el modelo matematico no lineal delcalentamiento de fibras de cristales fotonicos de PMMA que es estudiado en posteriorescapıtulos.

A continuacion se muestra una tabla en la que se recoge la nomenclatura utilizada eneste capıtulo y posteriores para describir los aspectos matematicos y termicos del modelodel sistema a estudiar.


Nomenclatura Sımbolos griegos

cp calor especıfico a presion constante α difusividadf, F funcion de la conductividad dimensional ε emisividad

y adimensional, respectivamente ρ densidadg,G funcion producto densidad por calor especıfico σ constante de Stefan-Bolzmann

dimensional y adimensional, respectivamente Θ temperatura adimensionalh coeficiente de transferencia de calor por Ω dominio computacional

por conveccion efectivo Ω region interna del dominio computacionalk conductividad termica ∂Ω frontera del dominio computacionall longitud ∆t escalon de tiemposM funcion adimensional que engloba los ∆Θ incremento de temperatura adimensional

terminos convectivo radiantes entre dos iteracionesn normal hacia fuera η normal hacia fuera adimensionalizadaq flujo de calorR radio de la preforma Subındicest tiempoT temperatura A, B, C VC adyacentes al VC consideradoTg temperatura de transicion vıtrea cr convectivo radianteVC volumen de control j VC perteneciente a ∂Ωx, X abscisa dimensional y adimensional (k) iteracion para un tiempo dadoy, Y ordenada dimensional y adimensional (n) instante de tiempo

i VC perteneciente a Ω0 valor referencia

3.1. Propiedades termicas del aire y del PMMA

El problema con el que nos enfrentamos viene determinado por la complejidad queaporta la no linealidad de la ecuacion debido tanto a la discontinuidad en espacio comoa la dependencia con la temperatura de las propiedades termicas del polımero y delaire contenido en los huecos de la preforma. Se hace necesario recopilar informacionsobre las caracterısticas de estas propiedades para poder desarrollar el metodo numerico.Existen modelos teoricos y datos experimentales para determinar la dependencia de laspropiedades con la temperatura, no obstante, lo que a priori se presumıa iban a serresultados similares dista de ser una realidad.

La conductividad termica, es una propiedad de los materiales que, excepto en el casode los gases a bajas temperaturas, no es posible predecir analıticamente; la informaciondisponible esta basada en medidas experimentales y en ciertas leyes propuestas. Laspropiedades termicas del aire a saber conductividad termica, calor especıfico y densidadestan bien documentadas y se pueden hallar en numerosos libros de texto, ver tablasiguiente. No ocurre lo mismo con las propiedades de los polımeros.


Propiedades termicas del aire [40]Temper. (oC) Dens. (kg/m3) Calor espec. (Jkg−1K−1) Conduct. (Wm−1K−1)

20 1.2042 1006.1 0.0256440 1.1273 1006.8 0.0271060 1.0596 1008.0 0.0285280 0.9996 1009.5 0.02991100 0.9460 1011.3 0.03127120 0.8979 1013.4 0.03261140 0.8544 1015.9 0.03394160 0.8150 1018.6 0.03525

3.1.1. Conductividad termica del PMMA

La conductividad termica es una propiedad importante en muchas operaciones de pro-cesado de polımeros, puesto que en la mayorıa de los procesos se hace necesario anadiro restar calor a estos polımeros. El control y optimizacion de las tecnicas de producciondependen en gran medida en un amplio conocimiento sobre la transferencia de calor. Nu-merosos investigadores han centrado sus estudios en determinar la conductividad termicade los polımeros teoricamente y experimentalmente. Caruthers et al. [41] revisaron losdatos, las teorıas y correlaciones existentes hasta 1998. No, obstante, a pesar de los es-fuerzos realizados, los metodos para predecir la conductividad de los polımeros no sonsatisfactorios.

Los posibles problemas para la evaluacion experimental de la conductividad termicade los polımeros es probablemente debido a los estrechos margenes de temperatura en losque se producen cambios en las propiedades fısicas del material. En efecto los rangos detemperatura para las distintas fases de los polımeros llamadas vidriosa, plastica, viscoso-plastica y lıquida, provocan que el mismo metodo experimental pueda no ser adecuadosobre todo el rango de interes. El desarrollo teorico tambien se enfrenta a varios problemasal depender la conductividad de los polımeros amorfos de numerosos factores, como sonlos constituyentes quımicos, las condiciones de procesado, tipo de estructura, distribucionde la densidad molecular, la temperatura, etc [42]. Este tipo de informacion no suele serdisponible con lo que no se pueden hacer estudios teoricos rigurosos.

Modelos teoricos

Un primer modelo para el calculo teorico de la conductividad en la region por debajode la temperatura de transicion vıtrea viene dado para los polımeros amorfos por la leyde Matthiessen [43], expresada por

1

k= A +

D

T(3.1)


donde A y D son parametros dependientes del tipo de polımero. El lımite de validez parala expresion anterior es relativamente alejada de la temperatura Tg. En esta region detemperaturas la dependencia de k viene determinada por la llamada variacion de la rutamedia libre del fonon. A temperaturas cercanas a la de transicion vıtrea y superiores,la conductividad como funcion de la temperatura para los polımeros amorfos no linealesse puede clasificar en tres grupos [42]: (i) los que tienen un claro valor maximo de laconductividad para la temperatura de transicion vıtrea Tg; (ii) los que presentan unameseta entorno a Tg; y (iii) los que muestran un comportamiento practicamente lineal enregiones vıtreas y elasticas con un cambio de pendiente en Tg.

Otro modelo teorico usado para la determinacion de la conductividad en los polımerosamorfos es la ecuacion de Debye [44]. Debido a la practica ausencia de efectos electronicosen la mayorıa de los polımeros, la conduccion de calor ocurre como resultado de lasvibraciones de la estructura reticular, por tanto, la conductividad termica teorica sepuede expresar mediante

k = (1/3)Cvl, (3.2)

donde C es la capacidad de calor especıfico por unidad de volumen, v es el la velocidadpromedio del fonon, y l es el camino libre del fonon. El valor de l para un polımero amorfoes una constante extremadamente pequena, debido a que un estado amorfo tiene multitudde defectos que producen dispersion de los fonones. Es conocido que la conductividad delos polımeros amorfos aumenta con la temperatura hasta Tg, mientras que esta decrecepor encima de Tg [44]. Esto se produce debido a que la conductividad termica de lospolımeros amorfos obedece a una dependencia de C con la temperatura por debajo deTg, mientras que las distancias de las cadenas de polımero por encima de Tg son demasiadolargas para que la conductividad termica respete la dependencia de v con la temperatura.Por tanto, el efecto de v en la conductividad termica es menor para cadenas de polımeroen un rango de distancias cortas, mientras que para rangos de distancia largas su efectosupera al de C [45].

Un tercer modelo para la estimacion de la conductividad de polımeros amorfos es eldado por el modelo Bicerano [46]. Este modelo reproduce mediante dos ecuaciones lacurva generalizada de van Krevelen [47], que relaciona el ratio de conductividades en unatemperatura dada y a Tg con el ratio de temperaturas a esa misma temperatura y Tg.Las ecuaciones del modelo para polımeros amorfos son las siguientes,

k(T ) = k(Tg)

(T

Tg

)0,22

, T ≤ Tg (3.3)

k(T ) = k(Tg)

[(1,2− 0,2

T

Tg

)]0,22

, T > Tg. (3.4)

No obstante, se requiere de la determinacion de la conductividad a la temperatura detransicion vıtrea que se puede aproximar mediante correlaciones en funcion de la estruc-tura del polımero [48] o experimentalmente.


Datos experimentales

El material utilizado en la fabricacion de las MPOF es el PMMA caracterizado por serun material termoplastico amorfo, incoloro de excelente transparencia optica, de trans-mision luminosa de aproximadamente el 92%, resistente a la abrasion y estable dimen-sionalmente aunque es fragil. Las medidas de su conductividad termica son complejasy se pueden realizar mediante distintas tecnicas. Usualmente, un metodo de medida es-pecıfico y sus instrumentos de medida tienen aplicacion para un determinado rango detemperaturas y de conductividades, en los que la incertidumbre ha de ser detallada. Otroaspecto importante a contemplar es la rigurosidad de las preparaciones pues es esencialya que pueden provocar errores de medida importantes.

Un total de 18 centros de investigacion europeos estudiaron, en 2003, la dependenciade la conductividad con la temperatura para una determinada muestra de PMMA me-diante diversos metodos de medicion, siendo los resultados recogidos por Rudtsch [49]. Lastecnicas usadas para medicion en regimen permanente fueron guarded hot-plate (GHP)y guarded heat-flow (GHF). Mientras que las medidas transitorias fueron llevadas a cabomediante las tecnicas transient hot-strip (THS) y transient plane-source (TPS). La den-sidad del polımero utilizado fue ρmean = 1185,0 kg ·m−3 con una incertidumbre ampliadapara las medidas (k=2) de U(ρmean) = 2,0 kg ·m−3. Los resultados obtenidos por los dis-tintos laboratorios para la conductividad presentan una alta variabilidad tanto en valorescomo en tendencia de las curvas en funcion de la temperatura. La tecnica utilizada parala determinacion de esta influye de forma decisiva en los resultados. Se comprueba portanto el grado de complejidad para obtener resultados normalizados y su importancia.

Es conocido que los polımeros pueden exhibir un comportamiento anisotropo en suconductividad termica, incluso en el caso de polımeros de estructura amorfa, caso delPMMA, cuando sus cadenas estan parcialmente alineadas entre sı [50]. Esto es debido aque el transporte de energıa termica se realiza mas eficientemente a traves de las cadenasdel polımero, las cuales consisten en fuertes enlaces covalentes carbono-carbono, queperpendicularmente a las cadenas de polımero donde la energıa termica es dirigida porlas debiles interacciones de van der Waals entre moleculas. Para polımeros sometidos aprocesos de estirado la anisotropıa de la conductividad termica aumenta.

Los resultados obtenidos para la conductividad termica del PMMA, por Rudtsch yHammerschmidt [49], en el rango de temperaturas entre 0-70 oC dados con un gradode incertidumbre de 1.7% son coherentes con los obtenidos por el NPL [51], mediantevacuum guarded hot-plate (VGHP) para el material de referencia Perspex (nombre decomercializacion del PMMA en Inglaterra). El VGHP permite medir la conductividadtermica para muestras de menos de 65 mm de grosor, conductividad menor a 2 W/(mK)y resistencia termica superior a 0.025 m2K/W para un rango de temperaturas compren-didads entre -50 oC y 70 oC. La tendencia de la curva de conductividad es monotonaascendente y practicamente lineal. No obstante, el rango de temperaturas no engloba alas temperaturas cercanas a la zona crıtica de transicion vıtrea donde el comportamiento


(a) (b)

Figura 3.2: Calor especifico del PMMA en funcion de la temperatura (a) Experimental [52]y Teorico [53]; (b) Recomendado por ATHAS databank [54].

muy probablemente es no lineal.

Por tanto, se comprueba que la medida de la conductividad para los polımeros amorfosno es una tarea sencilla por depender de multitud de factores. Una tecnica de medicionestandarizada debe de imponerse. Queda claro, no obstante, la fuerte dependencia dela conductividad con la temperatura pudiendo aproximarse con unos amplios margenesde error a un comportamiento lineal a bajas temperaturas que extrapolaremos hasta lastemperaturas levemente superiores a Tg.

3.1.2. Calor especıfico del PMMA

La determinacion del calor especifico del PMMA puede hacerse tanto experimental-mente como teoricamente. El concepto de las dislocaciones basado en el modelo meanderaplicado a polımeros amorfos, permite el calculo teorico del calor especifico del PM-MA [53]. Mientras que la determinacion experimental para polımeros a bajas tempera-turas suele realizarse mediante el uso del differential scanning calorimeter (DSC).

La figura 3.2.a muestra los valores experimentales obtenidos por Schawe et al. [52],para PMMA con peso molecular de 150kg/mol, ası como los obtenidos mediante la teorıade las dislocaciones ajustando los parametros a los resultados experimentales. Ambos re-sultados son muy similares en un amplio rango de temperaturas y muestran un compor-tamiento claramente no lineal conforme nos aproximamos a la zona de transicion vıtrea.Hay que destacar que el comportamiento del calor especıfico tiene histeresis puesto quesu valor difiere a una temperatura determinada en funcion de si se ha llegado por mediode calentamiento o enfriamiento, para temperaturas proximas a la transicion vıtrea.

Agari et al. [45] hallaron el calor especifico del PMMA mediante DSC obteniendosela misma morfologıa que los presentados en la figura 3.2.a, sin embargo, cuantitativa-


mente los valores son mayores, debido seguramente al uso de muestras de PMMA decaracterısticas distintas.

Resultados suplementarios obtenidos para la determinacion del calor especıfico delPMMA en la literatura se hallan en la base de datos ATHAS databank [54], y son cohe-rentes con los obtenidos por Theobald et al. [53] para el rango de temperaturas en el quese desarrolla el calentamiento de la preforma. En las zonas por debajo y por encima dela temperatura de transicion vıtrea, se puede aproximar el calor especıfico linealmente,mientras que a Tg ATHAS presenta una discontinuidad, ver figura 3.2.b. Se opta por unaaproximacion exponencial en el rango de temperaturas proximos a Tg, en consonanciacon los resultados presentados en los otros artıculos.

3.1.3. Otras propiedades del PMMA

El National Physical Laboratory ha realizado medidas de la emitancia proxima a lanormal para una muestra de PMMA mediante una tecnica propia de espectrofotometrıa yconvirtiendola a emitancia total hemisferica. Para un estrecho intervalo de temperaturasentre -25 oC y 33 oC se tienen valores en la emitancia para el PMMA que van desde 0.896hasta 0.899 al aumentar la temperatura [51], y por tanto muy debil dependencia respectode la temperatura.

Para el tipo de PMMA utilizado en la fabricacion de PCF, los experimentos indicanque las caracterısticas de calentamiento de las preformas de polimetilmatracrilato puedenser adecuadamente simuladas asumiendo que la preforma manifiesta un comportamientoproximo al de un cuerpo gris de emisividad ε = 0,96 al ser expuesto a las temperaturasde horno utilizado en el proceso de fabricacion [37].

Otra propiedad del PMMA como es la densidad tiene una dependencia con la tem-peratura. Esta decrece debilmente con la temperatura para temperaturas entre 130 oC y180 oC [45].

Nuevas tecnicas de medicion de las propiedades termicas van apareciendo conse-guiendose mayor exactitud en la valoracion de estas. No obstante, sigue siendo una tareaardua obtener resultados con un grado de exactitud elevado. No se puede obtener un re-sultado exacto para la variacion de la conductividad, del calor especıfico y de la densidaddel PMMA en funcion de la temperatura por las multiples variables que afectan a sudeterminacion. Obviamente, las muestras de PMMA no tienen todas las mismas propie-dades fısicas, ya sea por la polimerizacion llevada a cabo, las condiciones de procesado,los defectos de la muestra, etc. Sin embargo, con los datos recopilados tenemos una claraidea de los margenes en los que se mueven estos valores para el rango de temperaturasen el que estamos interesados.


(a) (b)

Figura 3.3: (a) Ejemplo de seccion transversal de la preforma de PCF. (b) Aproximaciondel Calor Especifico del PMMA en funcion de la temperatura (punteada [54], a trazos detipo exponencial) utilizados para el metodo numerico.

3.2. Formulacion del modelo del problema

El calentamiento de la preforma cilındrica de fibras de cristales fotonicos se realizaen un horno a su vez cilındrico cuyas paredes se encuentran a una temperatura Text. Lapreforma de PCF esta compuesta por una matriz de PMMA en la que se insertan huecosde aire dispuestos paralelos al eje longitudinal de la preforma, de forma que forman anilloshexagonales, ver figura 3.3.a. El aire localizado entre la preforma y las paredes del hornose considera que esta a la misma temperatura que las paredes del horno.

La difusividad del aire es dos ordenes de magnitud inferior a los del PMMA por tantoel PMMA es el medio controlante de la transferencia de calor. Por ello se desprecia laconveccion y la radiacion en los huecos de aire en los que el mecanismo considerado detransferencia de calor es unicamente la conduccion.

La solucion del problema bidimensional se computa en el dominio computacionalΩ = Ω ∪ ∂Ω siendo Ω el interior y ∂Ω su frontera. La seccion transversal de la preformade cristales fotonicos no tiene simetrıa axial por tanto se hace uso de ecuaciones en coor-denadas cartesianas. La ecuacion de conduccion del calor que modela el comportamientotermico en el interior de la preforma es

ρ(~x, T ) cp(~x, T )∂T (~x, t)

∂t−∇ · (k(~x, T )∇T (~x, t)) = 0 , t > 0, ~x ∈ Ω ⊂ IR2, (3.5)

donde k, ρ y cp son la conductividad, la densidad y el calor especıfico del material,respectivamente, T (~x, t) es el campo de temperatura, ~x son las coordenadas en el dominioΩ bidimensional y t es el tiempo. La condicion inicial es T (~x, 0) = T0.

Mientras la conduccion y la conveccion tienen lugar solo a traves de un medio material,la radiacion termica puede transportar la energıa en forma de calor a traves de un fluido


o del vacio, en forma de ondas electromagneticas que se propagan a la velocidad de laluz. Se tiene por tanto la condicion de contorno no lineal expresada por la ley de Stefan-Boltzmann propia de la radiacion entre las paredes del horno y las paredes de la preforma,y la conveccion dada por un coeficiente global

−k(~x, T )∂T (~x, t)

∂n= σε (T (~x, t)4 − T 4

ext) + h(T (~x, t)− Text), ~x ∈ ∂Ω, (3.6)

siendo ε la emisividad del material supuesta constante, σ la constante de Stefan-Bolzmann,h el coeficiente de transferencia de calor efectivo por conveccion y Text tanto la tempera-tura del aire que rodea a la preforma como la temperatura de las paredes del horno.

Tenemos por tanto una PDE con coeficientes variables y con condicion de contorno nolineal. Las propiedades termicas del material se modelan mediante funciones no linealesdependientes de la temperatura y de la posicion,

k = f(~x, T ), ρcp = g(~x, T ). (3.7)

Las propiedades termicas relativas al aire se aproximan mediante regresion lineal parala conductividad y regresion cuadratica para el producto densidad por calor especıfico apartir de la serie de datos dados en la tabla del apartado 3.1. En cuanto a las intrınsecasal PMMA se modelan mediante diversas funciones dependientes de la temperatura, co-herentes con lo visto en el apartado anterior. Las funciones utilizadas para aproximar laspropiedades termicas del PMMA son lineales para la conductividad, mientras que para elproducto densidad por calor especıfico se definen a trozos. En particular, esta funcion atrozos se modela para temperaturas en las cuales se produce la transicion vıtrea median-te una funcion exponencial, y por debajo y por encima de estas temperaturas mediantefunciones lineales, figura 3.3.b.

A continuacion se especifican las funciones que aproximan a las propiedades termicastanto del aire como del PMMA, usadas para el modelo.

f(T)=

6,905 ×10−5 T + 5, 479× 10−3 aire, 297 ≤ T ≤ 403K (3.8a)

4,528 ×10−4 T + 2, 550× 10−2 PMMA, 297 ≤ T ≤ 403K (3.8b)

g(T)=

8,540 ×10−3 T 2 − 8, 872 T + 3076,8 aire, 297 ≤ T < 403K (3.9a)

4,668 ×103 T + 2, 219× 105 PMMA, 297 ≤ T < 350K(3.9b)

1,639 ×105exp6,930×10−3 T PMMA, 350 ≤ T < 387K(3.9c)

2,831 ×103 T + 1, 302× 106 PMMA, 387 ≤ T ≤ 403K(3.9d)


Se adimensionalizan las coordenadas con respecto a R0, siendo R0 el radio carac-terıstico, X = x/R0, Y = y/R0 y η = n/R0 ; la temperatura Θ = (T − T0)/∆T siendoT0 la temperatura inicial y ∆T = Text − T0; las propiedades termicas son funcion de latemperatura y de la posicion (aire o PMMA) por tanto F ( ~X, Θ) = f(~x, ∆TΘ + T0) y

G( ~X, Θ) = g(~x, ∆TΘ + T0).

El modelo del problema adimensional esG( ~X, Θ) R2

0

∂Θ

∂t−∇ · (F ( ~X, Θ)∇Θ) = 0, para t > 0 y | ~X| < 1, (3.10a)

−F ( ~X, Θ)∂Θ

∂η= M(Θ), para t > 0 y | ~X| = 1, (3.10b)

siendo M(Θ) = σ εR0((∆TΘ + T0)4 − T 4

ext)/∆T + hR0(Θ − 1), y la condicion inicialΘ0 = 0.

Capıtulo 4

Metodo de volumenes finitos paraproblemas termicos no lineales

En el campo de la ingenierıa gran cantidad de procesos fısicos vienen gobernados porecuaciones en derivadas parciales. El uso de FVM para la resolucion de estas resultaaconsejable por las propiedades de estabilidad y convergencia que hacen del esquemade volumenes finitos un metodo robusto si la malla cumple ciertos criterios [15]. Enparticular se desarrolla el FVM de cuatro puntos en malla triangular no estructurada paraadaptarse a la complicada geometrıa de las PCF y resolver el problema del calentamientode preformas sometidas a un entorno convectivo radiante dado por el modelo adimensionaldescrito en el capıtulo anterior, ec (3.10).

Se ha comprobado el buen comportamiento de FVM, especialmente con uso de cir-cuncentro, para la resolucion de PDEs con coeficientes constantes. En el apartado 4.1 sedesarrolla el metodo de volumenes finitos para el modelo no lineal propuesto y se planteandos vıas de resolucion, la primera linealizando de la ecuacion a trozos vıa series de Taylory la segunda mediante un metodo implicito-explicito, o de retardo. En el apartado 4.2 semuestra el codigo implementado en Matlab y con el que se procede en sucesivos capıtulosa analizar el comportamiento termico de las PCF durante las primeras fases del estirado.

4.1. Metodo de los volumenes finitos en problemas

no lineales

En esta seccion se resuelve mediante FVM el modelo adimensional, ec. (3.10), dadopor una ecuacion en derivadas parciales parabolica cuyos coeficientes son variables y concondiciones de contorno no lineales. Se requiere para ello de la linealizacion en tiempomediante un proceso iterativo. Para obtener la solucion Θ(t + ∆t) a partir de Θ(t),

46

CAPITULO 4: FVM EN PROBLEMAS NO LINEALES 47

hacemos Θk=0(t) = Θ(t), y resolvemos en la iteracion (k+1)-enesima la siguiente ecuacion

G( ~X, Θk+1) R20

∂Θk+1

∂t−∇ · (F ( ~X, Θk+1)∇Θk+1) = 0, Θk+1(t) = Θk(t), (4.1)

con la condicion de contorno linealizada dada por

−F ( ~X, Θk+1)∂Θk+1

∂η= M(Θk) +

∂M(Θk)

∂Θk+1

(Θk+1 −Θk), (4.2)

obtenida por series de expansion de Taylor a partir de la ecuacion (3.10b) hasta terminosde segundo orden, p.e., O(Θk+1 − Θk)

2. Se trata de una condicion de contorno mixta ode Robin facilmente tratable. Tras la convergencia se tiene Θ(t + ∆t) = Θk+1(t + ∆t).

Mediante la siguiente tabla se muestra la nomenclatura que se va a utilizar posterior-mente para hacer mas legibles las expresiones matematicas que intervienen en el desarrollomatematico.

NomenclaturaFunciones de las propiedades termicas

Fi(k+1) = F ( ~Xi, Θi(k+1)) Derivadas parciales de las funciones

FiA(k+1) = (F ( ~Xi, Θi(k+1)) + F ( ~XA, ΘA(k+1)))/2 ∂Fi(k)/∂Θi = ∂Fi(k+1)(Θi(k))/∂Θi(k+1)

Gi(k+1) = G( ~Xi, Θi(k+1)) ∂FiA(k)/∂ΘA = ∂FiA(k+1)(Θi(k), ΘA(k))/∂ΘA(k+1)

Ms,j(k+1) = M( ~Xs,j, Θs,j(k+1)) ∂Gi(k)/∂Θi = ∂Gi(k+1)(Θi(k))/∂Θi(k+1)

Incremento de temperatura adimensional ∂Ms,j(k)/∂Θs,j = ∂Ms,j(k+1)(Θs,j(k))/∂Θs,j(k+1)

∆ΘA = ΘA(k+1) −ΘA(k)

En el dominio computacional Ω se define una malla formada por un conjunto detriangulos (volumenes de control) TΩ = T∂Ω+TΩ, donde T∂Ω es el subconjunto de triangu-los del contorno con al menos una arista perteneciente a ∂Ω y TΩ el resto de triangulosno situados en el contorno.

Integrando la forma conservativa de la ecuacion del calor (4.1) sobre un VC pertene-ciente a TΩ y aplicandole el teorema de la divergencia se tiene∫ ∫

R

G( ~X, Θk+1) R20

∂Θk+1

∂tdR−

∮S

F ( ~X, Θk+1)∇Θk+1 · d~s = 0, ∀ R ∈ TΩ, (4.3)

donde d~s = ~η ds, ~η es el vector normal unitario exterior al volumen R, y ds es el elementodiferencial de superficie.

El termino de la izquierda se evalua asumiendo que la temperatura en el nodo i-esimoes el valor medio para todo el volumen∫ ∫

R

G( ~X, Θk+1) R20

∂Θk+1

∂tdR = Gi(k+1) R2

0 Ai

∂Θi(k+1)

∂t, (4.4)

siendo Gi(k+1) = Gi( ~Xi, Θi(k+1)) y Ai el area del volumen de control. La segunda integralen la ecuacion (4.3) representa el flujo de calor saliente para las tres fronteras del volumen


de control de circuncentro (i). Este vector flujo es colineal al vector normal a la superficie,y para la frontera (bc) se representa por

−∫

bc

F ( ~X, Θk+1)∇Θk+1 · ~ηds ' −FiA(k+1)(ΘA(k+1) −Θi(k+1)

liA)lbc, (4.5)

siendo FiA(k+1) el valor medio de la conductividad adimensionalizada entre el volumen decontrol i y el volumen de control A para la iteracion k + 1.

4.1.1. Ecuaciones para los VC pertenecientes a TΩ

Para todos los volumenes de control pertenecientes a TΩ, es decir aquellos sin aris-tas frontera denotados por el circuncentro (i) en la figura 4.1.(b), la ecuacion (4.3) setransforma en

Gi(k+1) R20 Ai

dΘi(k+1)

dt− FiA(k+1)(

ΘA(k+1) −Θi(k+1)

liA)lbc

− FiB(k+1)(ΘB(k+1) −Θi(k+1)

liB)lca

− FiC(k+1)(ΘC(k+1) −Θi(k+1)

liC)lab = 0, (4.6)

una ecuacion diferencial ordinaria para cada volumen de control.

Tomando diferencias finitas en el tiempo y linealizando la ecuacion (4.6) (truncandoen terminos de segundo orden), se obtiene la siguiente expresion para el primer termino:

Gi(k+1) R20 Ai

dΘi(k+1)

dt' Gi(k+1) R2

0 Ai (Θi(k+1) −Θi(n)

∆t)

= Gi(k) R20 Ai (

Θi(k) −Θi(n)

∆t)

+

[∂Gi(k)

∂Θi

(Θi(k) −Θi(n)) + Gi(k)

]R2

0 Ai

∆t∆Θi, (4.7)

siendo ∆Θi = Θi(k+1) − Θi(k) el incremento de la temperatura adimensional entre dositeraciones sucesivas en un instante de tiempo dado y Gi(k) evaluada en (Θi(n) + Θi(k))/2.

La notacion Θi(k+1) es equivalente a Θ(n+1)i(k+1) y representa la temperatura adimensional en

la iteracion k + 1 para el paso de tiempo n + 1.


La ecuacion (4.6) para todos los VC de TΩ queda al linealizar los terminos conductivos[∂Gi(k)

∂Θi

(Θi(k) −Θi(n)) + Gi(k)

]R2

0 Ai

∆t∆Θi

−[∂FiA(k)

∂Θi

(ΘA(k) −Θi(k))− FiA(k)

]lbcliA

∆Θi

−[∂FiB(k)

∂Θi

(ΘB(k) −Θi(k))− FiB(k)

]lcaliB

∆Θi

−[∂FiC(k)

∂Θi

(ΘC(k) −Θi(k))− FiC(k)

]lab

liC∆Θi

−[∂FiA(k)

∂ΘA

(ΘA(k) −Θi(k)) + FiA(k)

]lbcliA

∆ΘA

−[∂FiB(k)

∂ΘB

(ΘB(k) −Θi(k)) + FiB(k)

]lcaliB

∆ΘB

−[∂FiC(k)

∂ΘC

(ΘC(k) −Θi(k)) + FiC(k)

]lab

liC∆ΘC

= −Gi(k) R20 Ai (

Θi(k) −Θi(n)

∆t) + FiA(k)

ΘA(k) −Θi(k)

liAlbc

+ FiB(k)

ΘB(k) −Θi(k)

liBlca

+ FiC(k)

ΘC(k) −Θi(k)

liClab. (4.8)

4.1.2. Ecuaciones para los VC pertenecientes a T∂Ω

A los volumenes de control pertenecientes al contorno T∂Ω, es decir aquellos con unaarista frontera denotados por el circuncentros j en la figura 4.1.b, se les aplica la con-dicion de contorno no lineal. Para evaluar el flujo de calor en la arista contorno (bc) esnecesario conocer la temperatura superficial. Por ello se introducen nuevas incognitasllamadas Θs,j(k+1) y localizadas en los puntos de los segmentos que conforman el perıme-tro de la preforma situadas a distancia mınima del circuncentro del volumen de controlconsiderado, ası se consigue que el vector flujo conductivo sea colineal con la normal alsegmento (bc).

La condicion de contorno linealizada (4.2), que representa el flujo de calor convectivoradiante, qcr, integrada en (bc), es

qcr = −∫

bc

Fj(k+1)∇Θk+1 · ~nds

= (Ms,j(k) +∂Ms,j(k)

∂Θs,j

∆Θs,j)lbc. (4.9)


(a) (b)

Figura 4.1: Ejemplo de volumen de control perteneciente al contorno T∂Ω y sujeto acondiciones de contorno conductivas hasta la superficie y convectivo radiantes en ella. (b)Notacion del mallado para VC de circuncentro (i) en TΩ y de circuncentro (j) en T∂Ω

Sea m1 el numero de volumenes de control y m2 el numero de temperaturas superficia-les incognita. Se tienen (m1 +m2) incognitas para solo m1 ecuaciones. Es necesario hacerbalance de flujos de calor en la frontera con lo que se obtienen m2 ecuaciones adicionales.El flujo conductivo entre el circuncentro de un volumen de control y la superficie es igualal flujo convectivo radiante en la superficie, figura 4.1.a. Por tanto se tiene al linealizarlola siguiente expresion

qcr ' −Fj(k+1)

Θs,j(k+1) −Θj(k+1)

ljslbc

= −Fj(k)

Θs,j(k) −Θj(k)

ljslbc −

[∂Fj(k)

∂Θj

(Θs,j(k) −Θj(k))− Fj(k)

]lbcljs

∆Θj

−Fj(k)lbcljs

∆Θs,j. (4.10)

Igualando las ecuaciones (4.9), (4.10) y haciendo

Φ(k) =

−∂Fj(k)

∂Θj

(Θs,j(k) −Θj(k)) + Fj(k)

∂Ms,j(k)

∂Θs,j

ljs + Fj(k)

, (4.11)

Ψ(k) =− Fj(k)(Θs,j(k) −Θj(k))−Ms,j(k)ljs

∂Ms,j(k)

∂Θs,j

ljs + Fj(k)

, (4.12)

se obtiene el incremento de temperatura adimensionalizada, para la iteracion (k + 1)-esima, del nodo superficial de un volumen de control (j) en funcion de la temperatura


del circuncentro del propio volumen de control (j)

∆Θs,j = Φ(k)∆Θj + Ψ(k). (4.13)

A partir de la ecuacion (4.3), para cada volumen de control de T∂Ω se obtiene lasiguiente ecuacion evaluando la incognita Θs,j(k+1) mediante la ecuacion (4.13).[

∂Gj(k)

∂Θj

(Θj(k) −Θj(n)) + Gj(k)

]R2

0 Aj

∆t∆Θj +

∂Ms,j(k)

∂Θs,j

Φklbc∆Θj

−[∂FjB(k)

∂Θj

(ΘB(k) −Θj(k))− FjB(k)

]lcaljB

∆Θj

−[∂FjC(k)

∂Θj

(ΘC(k) −Θj(k))− FjC(k)

]lab

ljC∆Θj

−[∂FjB(k)

∂ΘB

(ΘB(k) −Θj(k)) + FjB(k)

]lcaljB

∆ΘB

−[∂FjC(k)

∂ΘC

(ΘC(k) −Θj(k)) + FjC(k)

]lab

ljC∆ΘC

= −Gj(k) R20 Aj (

Θj(k) −Θj(n)

∆t) + FjB(k)

ΘB(k) −Θj(k)

ljBlca

+ FjC(k)

ΘC(k) −Θj(k)

ljClab − (Ms,j(k) +

∂Ms,j(k)

∂Θs,j

Ψk)lbc.

(4.14)

4.1.3. Resolucion del sistema

Mediante el desarrollo matematico anterior se ha llegado a un sistema de m1 ecuacio-nes lineales, formado por las ecuaciones (4.8) y (4.14) para cada volumen de control nofrontera (i) y frontera (j), respectivamente. Este sistema se puede representar de formamatricial mediante la siguiente expresion:

A(k)∆Θ = B(k), (4.15)

siendo ∆Θ el campo escalar del incremento de temperaturas adimensionalizadas de losvolumenes de control codificado en forma de vector, para la iteracion (k +1), A(k) es unamatriz cuadrada (m1 × m1), tipo dispersa (sparse), con tres y cuatro elementos por filapara los volumenes de control pertenecientes a T∂Ω y TΩ, respectivamente, y Bk un vectorcolumna (m1 × 1).

Se hace uso del mallado de Delaunay para conseguir VC con la mayor calidad posible,entendiendo por calidad su proximidad a la equilateralidad. Con ello se consigue quela practica totalidad de los circuncentros queden recluidos en el interior del VC que le


corresponde. Casos excepcionales de triangulos obtusos cuyos circuncentros se localizaranfuera del propio triangulo no alteran los resultados.

Este sistema se resuelve, mediante factorizacion LU resultando un metodo incondicio-nalmente estable con el valor inicial Θ0,0 = 0. La resolucion de la ecuacion del calor concondiciones de contorno no lineales expuesto en la seccion anterior se realiza mediante unproceso iterativo, a traves del cual se calcula la solucion para cada paso de tiempo delvector de tiempos.

El proceso iterativo es el siguiente, primeramente se inicializan las temperaturas detodos los nodos, tanto de los volumenes de control como de la superficie a su valor inicialΘ0. Para el instante (n) y la iteracion (k + 1) se evaluan Φ(k) y Ψ(k), funciones de Θs,j(k)

y por tanto de valor conocido, mediante las ecuaciones (4.11) y (4.12). Posteriormentese calculan los elementos de la matriz Ak y del vector Bk, dependientes de Φ(k) y Ψ(k),mediante las ecuaciones (4.8) y (4.14). A continuacion se construye el sistema matricial yresolviendo la ecuacion (4.15) se obtiene el vector solucion ∆Θk+1 para la iteracion (k+1).Por ultimo se calcula ∆Θs,j(k), que es funcion del vector Θk+1, con la ecuacion (4.13). Sila norma ||∆Θ||2 es inferior a una tolerancia, impuesta arbitrariamente de valor 10−6, seprocede a calcular la solucion para un instante de tiempo (n + 1). En caso contrario, seincrementa k y se repite el mismo proceso.

4.1.4. Metodo implıcito-explıcito

Se implementa otro metodo numerico para la resolucion del problema, simplificandolousando un retardo. La ecuacion que rige la dinamica del modelo es lineal, sin embargotiene coeficientes variables, y la condicion de contorno es no lineal por tanto es en generalimposible obtener soluciones analıticas para la ecuacion. Para su resolucion, se puedenevaluar las propiedades termicas, cp(T ) y k(T ), congelandolas en el paso de tiempo an-terior (n) al igual que la condicion de contorno y por tanto siendo las temperaturasconocidas. Se tiene ası un metodo implıcito-explıcito, tambien llamado congelado en eltiempo, que evita el tener que linealizar la ecuacion del contorno (3.10b). El modelo delproblema se simplifica por tanto, y la resolucion se realiza al igual que en el apartado an-terior aplicando el teorema de la divergencia y posteriormente FVM al sistema retardadosiguiente

G( ~X, Θ(n)) R20

∂Θk+1( ~X, t)

∂t−∇ · (F ( ~X, Θ(n))∇Θk+1( ~X, t)) = 0, | ~X| < 1(4.16a)

−F ( ~X, Θ(n))∂Θk+1

∂n= M(Θ(n)), | ~X| = 1. (4.16b)

La condicion de contorno es ahora lineal, de tipo Robin, y por tanto no hace falta iterar.Se resuelve para cada paso de tiempo y los resultados son utilizados para el paso detiempo posterior. En el proximo capıtulo se analiza la influencia de la simplificacion de


este metodo en los resultados obtenidos para problemas concretos. Obviamente, el usarel retardo va a generar errores que se cuantificaran.

4.2. Codigos Matlab para la resolucion de PDEs no

lineales mediante FVM

La implementacion del codigo para los metodos numericos desarrollados se describea grandes rasgos en este apartado mediante los diagramas 4.2 y 4.3.

Como etapa comun tanto al metodo iterativo apartados 4.1.1-4.1.3, como para elmetodo implıcito-explıcito del apartado 4.1.4, se recoge toda la informacion necesaria dela malla como son las longitudes de las aristas, las longitudes entre circuncentros, lasareas de los VC, los VC frontera, etc. Esto se realiza siguiendo el esquema desarrolladoen la seccion 2.3 con el programa llamado matrizT1 que almacena y organiza toda estainformacion de forma estructurada en una matriz de dimensiones (19× n).

Para el metodo iterativo empleado en la resolucion de PDEs no lineales, la segundaetapa de la implementacion consiste en construir el sistema matricial A(k)∆Θ = B(k) paracada iteracion y resolverlo bajo unas restricciones de tolerancia a especificar. Mientrasque el metodo implıcito-explıcito no requiere de iteraciones y calcula directamente lasolucion en cada instante de tiempo.

Los codigos de ambos metodos hacen llamadas a subprogramas que permiten evaluarel valor de las propiedades termicas como son el calor especıfico, la conductividad, y ladensidad, expresadas mediante las funciones F y G. Para ello se les pasa como parametroel tipo de compuesto (aire o PMMA) y el valor de la temperatura para cada iteracion y/opaso de tiempo ası como para cada volumen de control. En el caso del metodo iterativotambien se han de calcular las derivadas parciales de las funciones F y G respecto de lastemperaturas. En el CD-Rom que se adjunta con el proyecto se tiene acceso a todo loscodigos de todas las funciones implementadas.


Figura 4.2: Implementacion codigo iterativo.


Figura 4.3: Implementacion codigo implıcito-explıcito.

Capıtulo 5

Aplicacion y resultados de FVM enPDEs no lineales para produccion dePCF

El modelo matematico no lineal presentado en el capıtulo 4 y el metodo de resoluciondesarrollado permiten determinar el campo escalar de temperaturas de una preforma dePCF en la etapa inicial de estirado de fibras. Centramos el estudio en preformas cilındri-cas, con tres anillos hexagonales de huecos (figura 3.3.a) con distintas disposiciones, so-metidas a un calentamiento de tipo convectivo radiante en un horno. La dependencia dela conductividad y del calor especıfico con la temperatura, tanto para el aire como parael PMMA, son tenidas en cuenta por el modelo desarrollado, y presenta una compleji-dad respecto de otros modelos lineales utilizados en varios artıculos [39, 55]. Se estudiancomo afectan estos cambios introducidos comparando los resultados obtenidos con losexistentes en la literatura, tanto numericos como experimentales, para el problema con-creto de calentamiento de preformas de PCF. Posteriormente, se detalla este proceso decalentamiento y la influencia de las propiedades termicas del PMMA.

5.1. Validacion del metodo desarrollado

Antes de proceder a obtener resultados, es necesario comprobar que el metodo escorrecto. Partimos de una version del metodo basado en volumenes finitos que ha sidovalidado para su uso en PDEs lineales, capıtulo 3. En el caso del metodo basado envolumenes finitos adaptado para la resolucion de PDEs no lineales no se puede calcularla solucion analıtica, ni tampoco se dispone de un programa informatico que lo resuelva.Por tanto se recurre a una forma clasica de validacion del metodo.

El metodo desarrollado es consistente y convergente como se dijo previamente, por

56

CAPITULO 5: RESULTADOS FVM PARA PRODUCCION DE PCF 57

(a) (b)

Figura 5.1: Error absoluto maximo cometido por el metodo en funcion del parametro demalla, h, tomando como solucion numerica valida la obtenida para una malla caracteri-zada por (a) h=0.07 (b) h=0.04

tanto existe una ley que ha de seguir el error en funcion del parametro de malla [56], h,tomado. En efecto, el error es proporcional al parametro de malla elevado a una potencia,como muestra la expresion siguiente

Er ∝ A hq, (5.1)

y por tanto en escala logarıtmica se tiene una expresion lineal

ln Er ∝ ln A + q ln h. (5.2)

Para el calculo del error se requiere de la solucion exacta y de la solucion numerica.El resultado analıtico no esta disponible para este problema no lineal. Por tanto se haceuso de una solucion numerica calculada en un mallado extremadamente fino para hallarlos errores cometidos por las soluciones en otras mallas con menos nodos. Se toma comoprimera aproximacion por solucion valida la obtenida para una malla caracterizada porel parametro h=0.07 con 2184 volumenes de control y como segunda aproximacion laobtenida en una malla h=0.04 con 6730 llegando a los lımites tolerables de computo parael ordenador utilizado.

Los errores absolutos maximos obtenidos por mallas mas sencillas, tomando como re-ferencia las dos mallas mencionadas anteriormente, se muestran en la figura 5.1. Medianteregresion lineal se obtiene el valor de los parametros A y q que toman los valores que semuestran en la figura. En particular q(h=0,07)=2.334 y q(h=0,04)=2.418 con coeficientes decorrelacion para los dos casos de regresion superiores a 0.97.

El valor del parametro q da una idea del comportamiento del metodo. En amboscasos q es mayor que 2, por tanto se trata de una caso supercuadratico y la convergencia


(a) (b)

Figura 5.2: (a) Evolucion temporal de la temperatura del centro de la preforma en elmodelo lineal (lınea a trazos), y en el modelo no lineal (lınea continua). (b) Temperaturadel centro de la preforma medida experimentalmente (lınea continua) y calculada con elmodelo lineal (lınea a trazos) por Lyytikaınen [39]

del metodo parece demostrada. El metodo de los volumenes finitos es consistente enel sentido utilizado por el metodo de diferencias finitas solo respecto a los flujos perono a la solucion. No obstante, el metodo desarrollado proporciona resultados que sepueden considerar validos y a partir de los cuales comentaremos el comportamiento delas preformas de PCF en el estirado de fibras.

5.2. Comparativa de los resultados del modelo no li-

neal con los resultados experimentales

Se describen en primer lugar los resultados obtenidos con el modelo desarrollado,el cual incluye intercambios de calor radiantes en el horno, conduccion en un dominioheterogeneo y propiedades termicas variables, dependientes de la temperatura. El malladotriangular usado es de Delaunay y esta compuesto por aproximadamente 9000 triangulos.El paso de tiempos es variable, estando adaptado para obtener mayor resolucion en laszonas con mayores cambios en la solucion.

Los calculos realizados en esta seccion corresponden a preformas de fibras de cristalesfotonicos, fabricados en material PMMA, con las siguientes caracterısticas, R0 =0.025 m,con radio de huecos de 1 mm, radio de los tres anillos 0.009, 0.014 y 0.019 m, T0= 297 K yText=403. La conductividad del aire y del polımero vienen determinadas en el capıtulo 3.La emisividad del polımero tiene un valor ε =0.96, el coeficiente de transferencia de calorexterna de valor 8 Wm−2K−1.

La fase inicial del estirado de fibras concluye cuando la preforma ha sido calentada


en su totalidad por encima de Tg. El conocimiento de la evolucion de la temperaturadel nodo central de la preforma es de especial importancia puesto que es el de evolu-cion mas lenta y marca el tiempo crıtico a partir del cual se alcanza la temperatura detransicion vıtrea en toda la preforma. Se hace un estudio de la respuesta transitoria dela temperatura del centro de la preforma de PCF durante los primeros 160 minutos decalentamiento calculada mediante dos modelos distintos, uno lineal y otro no lineal, verfigura 5.2.a. Los resultados del modelo lineal, representado mediante lınea a trazos, consi-dera las propiedades termicas del aire y del PMMA constantes y por tanto independientesde la temperatura. Los resultados generados por el modelo no lineal desarrollado en es-te proyecto, representados mediante lınea continua, considera las propiedades termicasdependientes de la temperatura. Ambos resultados son practicamente identicos durantelos primeros 45 minutos a partir de los cuales se tiene una evolucion mas rapida de latemperatura para el modelo lineal con valores constantes de las propiedades del material.

Lyytikainen et al [39], obtuvieron resultados experimentales para este mismo problemade estirado de fibras y propusieron un modelo matematico lineal para simularlo. Losresultados experimentales y numericos obtenidos por estos investigadores se muestran enla figura 5.2.b. Para condiciones analogas a las presentadas anteriormente se muestranen la figura los resultados numericos del modelo lineal mediante lınea a trazos y losresultados obtenidos experimentalmente mediante lınea continua. Se observa, al igualque en la figura 5.2.b, una correspondencia practicamente total entre los resultadosexperimentales y los simulados durante los primeros 45 minutos. No obstante, a partir deese instante con temperaturas en el centro de la preforma de aproximadamente 90 oC, seobserva una clara separacion para los resultados. El modelo construido por Lyytikainenet al [39] fue ajustado mediante la seleccion de la conductividad del polımero, dando lasmejores aproximaciones para un valor de k=0.15 W/(mK).

El haber considerado las propiedades termicas del material de la preforma constan-tes no permite generar mejores ajustes a los resultados experimentales. La comparativaentre las figuras 5.2.a y 5.2.b muestra que los resultados del modelo de coeficientes varia-bles predice mejor los datos conseguidos en laboratorio, sin embargo, siguen existiendodiferencias remarcables.

Es posible modificar varios parametros de nuestro modelo no lineal para conseguirresultados mas acordes con los datos experimentales [39]. Entre las diversas opciones seprocede a variar la conductividad utilizada para el PMMA, dentro de los lımites recogi-dos en la literatura, que como se vio en el capıtulo 3 puede tener varias tendencias enfuncion de diversos factores. Se simula el calentamiento para tres modelos de conductivi-dad distintos que son linealmente ascendentes hasta la temperatura de transicion vıtrea,estimada en torno a los 90 oC. A partir de esta temperatura siguen tres distintos compor-tamientos lineales: ascendente, descendente o constante, ver figura 5.3.a. Los resultadosde las simulaciones para estas tres tendencias de la conductividad del polımero se mues-tran en la figura 5.3.b, donde se observa que la temperatura del centro de la preformapara el modelo con conductividades decrecientes a partir de Tg es menor a la del modelo


(a) (b)

Figura 5.3: (a) Posibles funciones de la conductividad utilizadas, influenciadas por latemperatura de transicion vıtrea (b) Evolucion temporal de la temperatura del centro dela preforma estructurada para distintas funciones de la conductividad.

con conductividad constante que es a su vez menor que el modelo con conductividadcreciente. Esto es logico puesto que la transferencia de calor disminuye al disminuir laconductividad termica. Por tanto el modelo con conductividad termica del polımero des-cendiente a partir de Tg es el que mejor se ajusta a los datos experimentales. Existen,no obstante, pequenos errores posiblemente provocados al no considerarse en el modeloperfiles de temperatura parabolicos en el horno lo que posiblemente disminuye levementela temperatura real de las paredes del horno.

Otras modificaciones del modelo, como pueden ser variar la curva del producto densi-dad - calor especifico en funcion de la temperatura del PMMA, tambien nos aproximan alos resultados experimentales pero a costa de alejarnos de las propiedades termicas de losmateriales. Se vio en otro apartado que los datos existentes acerca de estas propiedadesno son del todo concordantes entre los distintos estudios y por tanto se pueden variar ennuestro beneficio para mejorar el modelo.

Los responsables de la ralentizacion del proceso de calentamiento de la preforma,respecto del modelo determinado por la PDE lineal de coeficientes constantes, son elaumento del calor especıfico de forma exponencial en las cercanıas de Tg ası como eldescenso de la conductividad del PMMA. Esto es debido a una mayor necesidad deaporte de calor para satisfacer los cambios en las propiedades termicas del polımero.Sin embargo, la variacion de las propiedades termicas del aire, tenidas en cuenta en elmodelo no lineal, apenas altera los resultados puesto que es el polımero el compuestocontrolante del tiempo de calentamiento. Es decir, la dinamica de la transferencia decalor es controlada por el polımero de la matriz de la preforma por ser su difusividadmucho menor a la de los huecos de aire.

La figura 5.4.a muestra la temperatura adimensional para tres radios de preforma,R0=0.02, 0.025 y 0.03 m. Cuanto menor es el radio, mayores son las temperaturas tanto


(a) (b)

Figura 5.4: (a) Evolucion temporal del centro de la preforma y de la superficie. Lineasolida: R0=0.025; a trazos: 0.02; punteada: 0.03.(b) Evolucion temporal del centro de lapreforma y de la superficie. Linea solida: Text=400; a trazos: 450; punteada: 500.

en la superficie como en el centro de las distintas preformas. Esto responde a una mayornecesidad de aporte de calor cuanto mayor sea el volumen de material a calentar. Logi-camente, esto tambien significa que las preformas de menores dimensiones alcanzan Tg

antes y el proceso de produccion de fibras se acelera.

La temperatura adimensional para calentamientos de la preforma para tres tempera-turas de horno Text=403, 450 y 500 K, se visualizan en la figura 5.4.b. La desviacion entrelos tres resultados es debida a la no linealidad de las propiedades termicas del polımero,ya que de ser lineales las tres curvas tendrıan que ser identicas. Se obtienen en la su-perficie temperaturas adimensionales mayores cuanto mayor sea Text, mientras que en elinterior de la preforma los resultados son muy similares aunque siempre con evolucionesde temperatura mas rapidas al aumentar la temperatura del horno.

5.3. Analisis del calentamiento de la preforma

El calentamiento de la preforma depende en gran medida del volumen ocupado por loshuecos de aire puesto que la temperatura en ellos evoluciona de forma mucho mas rapidaque en el polımero. Este hecho se puede observar para los nodos que recorren un eje de lapreforma que atraviesa tres huecos de aire, figura 5.5.a, para los cuales se muestra el perfilradial de temperaturas en distintos instantes de tiempo, figura 5.5.b. Las temperaturasevolucionan de forma rapida cerca de la superficie mientras que hasta que la conduccionde calor hacia el centro de la preforma de forma mas lenta. El gradiente de temperaturasen los huecos es muy importante en comparacion con las zonas anexas del polımero. Estoes debido a que la conductividad del aire es aproximadamente cinco veces menor que ladel PMMA. No obstante, el producto densidad por calor especıfico es tres ordenes demagnitud inferior para el aire por lo que la difusividad (α = k/ρcp) en los huecos de aire


(a) (b)

(c) (d)

Figura 5.5: (a) La figura muestra los nodos en los que se representa la solucion. Los nodosrecorren un eje y pasan por tres huecos de aire. (b) Distribucion radial de temperaturascomo funcion del tiempo. Las zonas sombreadas representan las localizaciones de loshuecos. (c) Flujo de calor aportado de tipo radiante (continua) y de tipo convectivo (atrazos) en funcion del tiempo. (d) Error cometido entre el metodo retardado y el metodolinealizado a trozos

es mucho mayor que en el polımero y por tanto en los huecos la transferencia de calor seproduce de manera mas efectiva.

No solo la proporcion de cada compuesto determina el comportamiento frente al ca-lentamiento de la preforma, tambien es importante el estudio de la localizacion de loshuecos teniendose evoluciones mas rapidas para las temperaturas cuanto mas proximosal centro se localicen los huecos [55]. Se ha comprobado que para determinadas estructu-ras de preformas de PCF de anillos hexagonales el calentamiento es mas lento que parauna preforma de dimensiones equivalente sin huecos de aire [39, 55], resultado a prioriinesperado teniendo en cuenta las difusividades relativas de ambos compuestos. En efec-to, resulta sorprendente la evolucion de la temperatura para una preforma compacta dePMMA sea mas rapida que para una preforma con huecos de aire cerca de la superficie.No obstante, al evolucionar muy rapidamente la temperatura cerca de la superficie debido


a la accion de los huecos, el flujo de calor radiante incidente en la superficie disminuyepor ser el gradiente de temperaturas con el horno menor lo que redunda en un menoraporte de calor. La mejor difusividad no compensa esta disminucion de flujo y por tantose retarda el calentamiento.

Dos mecanismos de transferencia de calor aportan la energıa a la preforma: la radia-cion y la conveccion. La figura 5.5.c muestra el flujo de calor aportado por ambos enfuncion del tiempo. Se observa que la radiacion aporta mas calor que la conveccion queha sido modelada mediante un coeficiente global h. En efecto un 61.3%, del calor totalcuantificado en 37.4 kJ/m2, es aportado por la radiacion. Obviamente, cuanto mayor seala temperatura de las paredes del horno, mayor sera el ratio de calor aportado por laradiacion puesto que este aumenta con la temperatura a la cuarta potencia mientras quela conveccion aumenta linealmente con la temperatura.

En el apartado 4.1.4 se desarrollo una variante para el metodo numerico que llamadoimplıcito-explıcito en tiempos. Es interesante evaluar el error que produce la resolucionde un problema real, como es el caso del calentamiento de PCF, mediante este metodo encomparacion con el metodo iterativo de linealizacion a trozos implementados en la seccionanterior. La figura 5.5.d muestra el error en la temperatura adimensional para el nodocentral de la preforma. Se observan errores maximos de entorno al 2% alcanzandose tras51 minutos de calentamiento que es cuando la mayor parte de la preforma se encuentra atemperaturas proximas a Tg. En regimen permanente, logicamente, los errores entre losdos metodos son nulos para este problema de calor sin generacion.

Las figuras 5.6.a y 5.6.b muestran la distribucion de la temperatura, de la conductivi-dad, del producto calor especifico por densidad y de la difusividad para los tiempos t=20min y t=40 min respectivamente. Se tiene para ambos tiempos de calentamiento unadistribucion de temperaturas en donde los mayores gradientes se localizan en los huecosde aire. Por otro lado, la conductividad y el producto densidad por calor especıfico en loshuecos es mucho menor a la del resto del polımero mientras que la difusividad en ellos esmucho mayor.

La localizacion de las zonas en las cuales se produce la transicion vıtrea se apreciaclaramente en las graficas de la conductividad puesto que se observa un maximo de estapropiedad en | ~X| ' 0,85 y | ~X| ' 0.4, a los 20 minutos y 40 minutos, respectivamente.Las conductividades para las regiones externas a la zona de transicion vıtrea disminuyenlo cual retarda el proceso de calentamiento. De igual modo el producto de la densidad porel calor especifico al aumentar de forma exponencial su valor en el PMMA, en regionesde transicion vıtrea, provoca una mayor demanda de energıa por lo que se retarda elcalentamiento de la preforma.

Por ultimo, se observa a los 40 minutos valores de difusividad menores que a los 20minutos tanto en zonas internas de la preforma como en zonas externas. Esto indica que ladependencia de las propiedades termicas con la temperatura provocan una ralentizacionde todo el proceso de calentamiento. Sobre todo debido a los importantes cambios fısicos


(a) Tiempo t=20 min

(b) Tiempo t=40 min

Figura 5.6: Distribucion espacial de temperatura oC (arriba izquierda), conductividadW/(m K)(arriba derecha), producto densidad por calor especıfico J/(m3 K)(abajo iz-quierda) y difusividad m2/s (abajo derecha) en la preforma a los (a) 20 minutos, (b) 40minutos.


que se producen por la transicion vıtrea.

En resumen, se ha desarrollado un nuevo modelo para calentamientos de PCF basadoen PDEs no lineales. El modelo tiene en cuenta la complicada geometrıa de la prefor-ma estructurada por numerosos huecos de aire ası como la variacion de las propiedadestermicas con la temperatura. Tanto la conductividad como el calor especifico del materialpolımero constituyente de la preforma tiene comportamiento claramente no lineal en lasproximidades de la transicion vıtrea. Los resultados obtenidos por este modelo se ajustande manera mucho mejor a los resultados experimentales realizados por otros investiga-dores y permiten controlar de forma mas eficaz todo el proceso de calentamiento de lapreforma de PCF.

Capıtulo 6

Conclusiones

En este capıtulo se presentan las principales conclusiones sobre la resolucion de proble-mas regidos por ecuaciones en derivadas parciales tanto lineales como no lineales medianteel metodo de volumenes finitos. En primer lugar se exponen los aportes del proyecto encuanto a desarrollo de una tecnica numerica de discretizacion espacial como es el FVMy la forma de implementacion llevada a cabo. Posteriormente se muestran los resultadosobtenidos mediante el codigo implementado para distintos problemas.

6.1. Implementacion del metodo de volumenes fini-

tos

Usualmente se viene utilizando el metodo de los elementos finitos en los softwarescomerciales para resolver los distintos tipos de problemas gobernados por ecuacionesen derivadas parciales. En efecto, FEMLAB o MATLAB son dos ejemplos de paquetesinformaticos que integran el FEM para obtener resultados de este tipo de problemas entoda clase de regiones de calculo. En el caso de este proyecto se ha desarrollado un metodoalternativo basado en volumenes finitos. Se han mencionado las ventajas que ofrece estatecnica en cuanto a convergencia, estabilidad y conservacion.

El metodo requiere en primer lugar de la discretizacion del dominio de computo quehemos considerado bidimensional, aunque el metodo es perfectamente aplicable para mor-fologıas tridimensionales. La discretizacion idonea para FVM en recintos complejos es latriangularizacion de Delaunay que genera volumenes de control de alta calidad. A partirde estos datos, se ha disenado una estructura de almacenamiento de los datos geometricosindispensables para implementar el metodo. Entre ellos se encuentran la forma de codi-ficar cada uno de los triangulos y nodos de estos, ası como las longitudes de las aristas,area de cada VC, distancia entre nodos de los VC, etc. Para ello partimos de las mallasgeneradas por la PDEToolbox de Matlab, pero preparadas para un metodo de elementos

66

CAPITULO 6: CONCLUSIONES 67

finitos y por una tecnica vertex-centered mientras que la variante elegida de FVM es detipo cell-centered. De las distintas formas de reutilizar la informacion se elige una querequiera un menor coste computacional.

Con toda informacion se puede implementar de forma eficiente el FVM. Las formulasque se obtienen en el desarrollo de este metodo permiten construir un sistema matricial detipo disperso. Debido a la notacion inherente de la malla no estructurada estas matricesno tienen sus elementos no nulos en varias diagonales, incrementandose la complejidad dela implementacion. Ademas existen dos tipos de VC, los que poseen una arista fronteray los que carecen de ella. El tratamiento de cada uno de ellos es distinto. Para el tipo deproblemas de PDEs parabolicos que hemos resuelto, una vez realizada la discretizacionespacial, se realiza una discretizacion temporal que nos permite obtener la solucion paracada VC en cada paso de tiempo elegido.

Dentro de las variantes existentes en FVM se han implementado con uso de baricentroy de circuncentro como nodo representativo de cada VC. Las expresiones matematicas seven afectadas por el cambio al igual que los resultados. En efecto, los errores numericosobtenidos por el FVM baricentrico se encuentran por encima de lo esperado. Mas auncuando el comportamiento de los resultados generados por el FVM circuncentrico tienenla misma calidad que los de un programa comercial, tal como se ha demostrado en variascomparativas.

Por tanto se ha creado un codigo alternativo al existente de FEM implementadopor programas comerciales que permite resolver problemas PDEs tanto lineales como nolineales.

6.2. Resultados obtenidos

6.2.1. Resultados en problemas lineales

Una primera lınea del trabajo ha sido desarrollar el FVM para problemas parabolicosde tipo lineal. Este tipo de problemas que presentan la ventaja de obtener la solucionanalıtica, bajo ciertas condiciones, permiten determinar la exactitud del metodo numeri-co. Por ello ha sido posible comprobar el funcionamiento de varias versiones de FVM, yajustar los distintos parametros que influyen de forma decisiva en la obtencion de buenosresultados. Entre estos parametros se encuentran el parametro de malla que es el tamanomaximo de arista de los triangulos de la malla, y el paso de tiempos en el cual se integrala PDE. Una seleccion incorrecta del valor de estos parametros provoca oscilaciones nodeseadas en los resultados.

Se ha realizado un analisis detallado de los resultados del problema del calor concondiciones de contorno tipo Dirichlet. Mas concretamente de los resultados obtenidos porlas variantes circuncentricas y baricentricas del FVM en comparacion con los resultados


obtenidos mediante la PDEToolbox de Matlab. El FVM es un metodo conservativo y portanto presenta ventajas respecto al FEM. Los resultados obtenidos para la ecuacion delcalor mediante FVM con uso del circuncentro generan resultados con la misma precisionque los obtenidos mediante software comercial. No ocurre lo mismo si se hace uso delbaricentro. En este caso, los resultados dependen fuertemente del tipo de malla usada yse genera una dispersion en los resultados que es exagerada. A partir de estos resultadosse elabora el artıculo adjunto en el apendice.

6.2.2. Resultados en problemas no lineales

Si bien existen problemas gobernados por PDEs lineales, suele ser una realidad queuna gran parte de los procesos fısicos tienen componentes no lineales. En el proyectose describe el proceso de estirado de preformas de fibras de cristales fotonicos, en elcual una preforma cilındrica de polımero es introducida en el horno y calentada hastatemperaturas por encima de la transicion vıtrea. Este proceso viene dado por la PDEparabolica de la conduccion del calor. Pero es importante resenar que las condiciones decontorno son de tipo no lineal debido a la accion de la radiacion expresada mediante laley de Stefan - Boltzman. Tambien las propiedades de los materiales son dependientes dela temperatura de forma no lineal. Por tanto la resolucion del problema ha de recurrir aun proceso iterativo que permita linealizar el modelo construido.

Una vez disenado e implementado el metodo numerico se procede a realizar compara-ciones de los resultados obtenidos para el problema del calentamiento de preformas conlos resultados existentes en la literatura. En particular se comparan con los resultadosexperimentales obtenidos por un grupo de investigacion en la Universidad de Sydney. Elmodelo no lineal construido permite un ajuste excelente a los resultados experimentales,mientras que el modelo lineal que ha desarrollado el grupo de Sydney y que no contemplala no linealidad en las propiedades termicas del PMMA no genera resultados tan precisos.

La simulacion del calentamiento bajo distintas condiciones permite un mejor entendi-miento del proceso de calentamiento de la preforma. En particular se han variado tantolas dimensiones de la preforma como la temperatura del horno y estos resultados puedenincidir en la mejora del proceso de produccion. Asegurando un menor tiempo de produc-cion y una mejor calidad de las fibras en las que a menudo se dan colapsos de los huecosde aire por no cumplirse los requisitos del calentamiento. La etapa en la que se producela transicion vıtrea en el polımero es crucial por las enormes no linealidades que produceen el modelo matematico. Nuestro modelo las contempla y por tanto asegura mejoresresultados que aquellos descritos en otros trabajos.


6.3. Nuevas lıneas de investigacion

Las fibras de cristal fotonico son una nueva tecnologıa que se viene desarrollando estosultimos anos y con expectativas realmente optimistas para dar un salto en el control dela luz. Por ello numerosas investigaciones vienen desarrollandose sobre la influencia de loshuecos en la transmision de la luz. Nuevos materiales y sobre todo nuevas geometrıas ydisposiciones de los huecos son las principales innovaciones. Resulta esencial seguir inves-tigando sobre los procesos termicos que se producen en el horno con el fin de conseguirfibras de alta calidad que soporten el proceso de estirado en la cual los grandes tamanosde huecos y la poca cantidad de material (PMMA o sılice) confieren poca resistencia alas fibras.

Los metodos implementados en este proyecto aportan un valor anadido respecto deotros ya existentes por tener en cuenta las no linealidades de las propiedades termicas delmaterial. No obstante, la complejidad de la realidad fısica del proceso ha sido limitadaen nuestro modelo. En efecto futuras lıneas investigacion tendran que contemplar las tresdimensiones de la preforma. Para preformas compactas y propiedades termicas constantesya se han obtenido resultados tridimensionales [37], pero sin tener en cuenta la estructurainterna de las PCF. Igualmente se ha obviado la radiacion interna que se produce en loshuecos de aire. Aun queda un largo camino de investigacion por recorrer, tanto dentrodel area de la optica como de la termica que permitan hacer de esta nueva tecnologıaun estandar para la comunicacion. Este proyecto ha sentado unas bases que para laresolucion de problemas no lineales gobernados por PDEs mediante FVM mediante elque se han obtenido resultados acordes con la experimentacion. El metodo desarrolladopuede extrapolarse a diversos campos tanto ingenieriles, como biologicos o economicos.

Bibliografıa

[1] T. Puu, Nonlinear Economic Systems, Mir, Moscow (1999).

[2] N.A. Magnitskii and S.V. Sidorov, “Distributed model of a self-developing marketeconomy,”Computational Mathematics and Modeling, Vol. 16, pp. 83–97 (2005).

[3] M. Fujita, P. Krugman, and T. Venables, The spatial economy: cities, regions, andinternational trade, MIT Press (1999).

[4] I. Sapariuc, M.D. Marcozzi, and J.E. Flaherty, “A numerical analysis of variationalvaluation techniques for derivative securities,”Applied Mathematics and Computa-tion, Vol. 159, pp. 171–198 (2004).

[5] G.J. Reid, A.D. Wittkopf, and A. Boulton, “Reduction of systems of nonlinear partialdifferential equations to simplified involutive forms,”European Journal of AppliedMathematics, Vol. 7, pp. 604–635 (1996).

[6] E. Bertolazzi and G. Manzini, “A cell-centered second-order accurate finite volu-me method for convection-diffusion problems on unstructured meshes,”MathematicalModels and Methods in Applied Sciences, Vol. 14, pp. 1235–1260 (2004).

[7] M. Ohlberger, “A posteriori error estimates for vertex centered finite volume appro-ximations of convection-diffusion-reaction equations,”Mathematical Modelling andNumerical Analysis, Vol. 35, pp. 355–387 (2001).

[8] Q. Du and L. Ju, “Numerical simulations of the quantized vortices on a thin super-conducting hollow sphere,”Journal of Computational Physics, Vol. 201, pp. 511–530(2004).

[9] Q. Du and D. Wang, “The optimal centroidal Voronoi tessellations and the gers-ho’s conjecture in the three-dimensional space,”Computers and Mathematics withApplications, Vol. 49, pp. 1355–1373 (2005).

[10] F. Aurenhammer and R. Klein, Handbook of Computational Geometry, North-Holland, Amsterdam (2000).

70

71

[11] Z.C. Li and S. Wang, “The finite volume method and application in combina-tions,”Journal of Computational and Applied Mathematics , Vol. 106, pp. 21–53(1999).

[12] J. Alberty, C. Carstensen, and S. Funken, “Remarks around 50 lines of Matlab: shortfinite element implementation,”Numerical Algorithms, Vol. 20, pp. 117–137 (1999).

[13] J. Alberty, C. Carstensen, S. Funken, and R. Klose, “Matlab implementation of theFinite Element Method in Elasticity,”Computing, Vol. 69, pp. 239–263 (2002).

[14] C. Hirsch, Numerical Computation of Internal and External Flows, Vol. 1, Wiley(1988).

[15] R. Herbin and O. Labergerie, “Finite volume schemes for elliptic and elliptic-hyperbolic problems on triangular meshes,”Computational Methods in Applied Me-chanics and Engineering, Vol. 147, pp. 85–103 (1997).

[16] S. Chou and S. Tang, “Comparing Two Approaches of Analyzing Mixed Finite Vo-lume Methods,”Journal of the KSIAM, Vol. 5, pp. 55–78 (2001).

[17] S. V. Patankar, Numerical heat transfer and fluid flow, Taylor & Francis (1980).

[18] L. Langemyr et al., Partial Differential Equation Toolbox User”s Guide , The MathWorks (1996).

[19] H C Carslaw and J C Jaeger, Conduction of Heat in Solids, Clarendon Press, Oxford(1986).

[20] A. Ferrando, E. Silvester, and P. Andres, “Designing the properties of dispersion-flattened photonic crysal fibers,”Optics Express, Vol. 9, pp. 687–697 (2001).

[21] A. Ortigosa-Blanch, J.C. Knight, W.J. Wadsworth, J. Arriaga, B.J. Mangan, T.A.Birks, and P.S. Russel, “Highly birefringent photonic crystal fibers,”Optics Letters,Vol. 25, pp. 1325–1327 (2000).

[22] P. Petropoulos, T.M. Monro, W. Belardi, K. Furusawa, J.H. Lee, and D.J. Ri-chardson, “2R-regenerative all-optical switch based on a highly nonlinear holey fi-ber,”Optics Letters, Vol. 26, pp. 1233–1235 (2001).

[23] J.C. Knight, “Photonic crystal fibres,”Nature, Vol. 424, pp. 847–851 (2003).

[24] M.C.J. Large, M.A. van Eijkelenborg, A. Argyros, J. Zagari, S. Manos, N.A. Issa,I. Bassett, S. Fleming, R.C. McPhedran, C.M. de Sterke, and N.A.P. Nicorovici,Microstructured polymer optical fibres: a new approach to POFs, POF Conference,Amsterdam (2001).

72

[25] M.A. van Eijkelenborg, M.C.J Large, A. Argyros, J. Zagari, S. Manos, N.A. Issa, S.C.Fleming, R.C. McPhedran, C. M. Sterke, and N. A. P. Nicorovici, “Microstructuredpolymer optical fibre,”Optics Express, Vol. 9, pp. 319–327 (2001).

[26] A. Argyros, I.M. Bassett, M.A. van Eijkelenborg, M.C.J. Large, J. Zagari, N.A.P.Nicorovici, R.C. McPhedran, and C.M. Sterke, “Ring structures in polymer micros-trutured polymer optical fibres,”Optics Express, Vol. 9, pp. 813–820 (2001).

[27] G. Barton, M.A. van Eijkelenborg, G. Henry, M.C.J. Large, and J. Zagari, “Fabri-cation of microstructured polymer optical fibres,”Optical Fiber Technology, Vol. 10,pp. 325–335 (2004).

[28] N.P. Bansal and R.H. Doremus, Handbook of Glass Properties, New York Academic(1986).

[29] S. Wu, “Surface and interfacial tensions of polymer melts. IIPoly(methylmethacrylate), poly(n-butyl-methacrylate), and polystyrene,”Journal ifPhysical Chemistry, Vol. 74, pp. 632–638 (1970).

[30] A.D. Fitt, K. Kurusawa, T.M. Monro, and C.P. Please, “Modeling the fabricationof hollow fibers: capillary drawing,”Journal Lightwave Tecnology, Vol. 19, pp. 1924–1931 (2001).

[31] N.A. Issa, M.A. van Eijkelenborg, G. Henry, M. Fellew, and M.C.J. Large, “Fa-brication and characterization of microstructured optical fibres with elliptical ho-les,”Optics Letters, Vol. 23, pp. 2255–2266 (2005).

[32] S.C. Xue, R.I. Tanner, G.W. Barton, R. Lwin, M.C.J. Large, and L. Poladian, “Fabri-cation of microstructured optical fibers - Part II: numerical modeling of steady-statedraw process,”Journal of Lightwave Tecnology, Vol. 23, pp. 2255–2266 (2005).

[33] J. I. Ramos, “Drawing of annular liquid jets at low Reynolds num-bers,”Computational Theoretical Polymer Science, Vol. 11, pp. 429–443 (2001).

[34] J. I. Ramos, “Non linear dynamics of hollow, compound jets at low Reynolds num-bers,”International Journal of Engineering Science, Vol. 39, pp. 1289–1314 (2001).

[35] A.L. Yarin, P. Gospodinov and V.I. Roussinov, “Stability loss and sensitivity inhollow,”Physics Fluids, Vol. 6, pp. 1454–1463 (1994).

[36] P. Gospodinov and A.L. Yarin, “Draw resonance of optical microcapillaries in hollowfiber drawing,”International Journal of Multiphase Flow, Vol. 23, pp. 967–976 (1997).

[37] H.M. Reeve, A.M. Mescher, and A.F. Emery, “Experimental and numerical investi-gation of polymer preform heating,”Journal of Materials Processing and Manufac-turing Science, Vol. 9, pp. 285–301 (2001).

73

[38] J.I. Ramos, “Convection and radiation effects in hollow, compound optical fi-bers,”International Journal of Thermal Sciences, Vol. 44, pp. 832–850 (2005).

[39] K. Lyytikainen, J. Zagari, G. Barton, and J. Canning, “Heat transfer within a mi-crostructured polymer optical fibre perform,”Modelling and Simulation in MaterialScience and Engineering, Vol. 12, pp. S255–S265 (2004).

[40] A.J. Chapman, Heat Transfer, 3rd edition, MacMillan Publishing Co., New York(1974).

[41] J.M. Caruthers, K.-C. Chao, V. Venkatasubramanian, R. Sy-Siong-Kiao, C.R. Nove-nario, and A. Sundaram, Handbook of Diffusion and Thermal Properties of Polymersand Polymer Solutions, DIPPR, New York (1998).

[42] P. Dashora and G. Gupta, “On the temperature dependence of the thermal conduc-tivity of linear amorphous polymers,”Polymer, Vol. 37, pp. 231–234 (1996).

[43] W. Hayden, W.G. Moffatt, and J. Wolff, The structure and properties of materials,Wiley Eastern, New Delhi (1968).

[44] Y.K. Godovsky, Thermophysical Properties of Polymers, Springer, New York (1992).

[45] Y. Agari, A. Ueda, Y. Omura, and S. Nagai, “Thermal diffusivity and conductivityof PMMA/PC blends,”Polymer, Vol. 38, pp. 801–807 (1997).

[46] J. Bicerano, Prediction of Polymer Properties, 2nd Edition, Marcel Dekker, NewYork (1996).

[47] D.W. van Krevelen, Properties of Polymers. Their Correlation with Chemical Struc-ture, Their Numerical Estimation and Prediction from Additive Group Contribu-tions, 3rd Edition, Elsevier, Amsterdam (1990).

[48] C. Zhong, Q. Yang, W. Wang, “Correlation and prediction of the thermal conducti-vity of amorphous polymers,”Fluid Phase Equilibria, Vol. 181, pp. 195–202 (2001).

[49] S. Rudtsch and U. Hammerschmidt, “Intercomparison of Measurements of the Ther-mophysical Properties of Polymethyl Methacrylate,”International Journal of Ther-mophysics, Vol. 25, pp. 1475–1482 (2004).

[50] K. Kurabayashi, “Anisotropic Thermal Properties of solid Polymers,”InternationalJournal of Thermophysics, Vol. 22, pp. 277–288 (2001).

[51] C. Stacey, “NPL vacuum guarded hot-plate for measuring thermal conductivity andtotal hemispherical emittance of insulation materials,”ASTM Special Technical Pu-blication, Vol. 1426, pp. 131–144 (2002).

74

[52] J.E.K. Schawe and S. Theobald, “Linearity limits of dynamic calorimetric responseat the glass transition of polystyrene ,”Journal of Non-Crystalline Solids, Vol. 235-237, pp. 496–503 (1998).

[53] S. Theobald, W, Pechhold, and B. Stoll, “The pressure and temperature dependenceof the relaxation processes in poly(methylmethacrylate),”Polymer, Vol. 42, pp. 289–295 (2001).

[54] B. Wunderlich, “The Athas Data Base on Heat Capacities of Polymers,”Pure andApplied Chemistry, Vol. 67, pp. 1019–1026 (1995).

[55] F.R. Villatoro and P. Udias de la Mora, “Numerical Analysis of the Heat TransferWithin Photonic Crystal Fibre Preforms,”Proceedings of SPIE, Vol. 5840, pp. 397–408 (2005).

[56] R. Eymard, T. Gallouet, and R. Herbin, The Finite Volume Method, to appear inHandbook of Numerical Analysis, Ciarlet and Lions eds, North Holland (2006).

[57] D.A. Field, “Qualitative measures for initial meshes,”International Journal for Nu-merical Methods in Engineering, Vol. 47, pp. 887–906 (2000).

Apendice A

Contenido del CD

A.1. Contenidos

El CD-ROM que se adjunta proporciona la informacion necesaria para que el lectorpueda repetir los resultados obtenidos con los codigos implementados. El CD esta estruc-turado en cuatro carpetas.

Carpeta MEMORIA: Proporciona la memoria del proyecto en formato .pdf, ejecutablecon el programa Acrobat Reader y su version en LaTeX.

Carpeta SIMULACION 3D: Presenta la respuesta transitoria del calentamiento parauna preforma de PCF mediante archivos de video. Se representa la evolucion de latemperatura en cada punto de las preformas mediante figuras 3D.

Carpeta RESULTADOS: Contiene varios archivos .mat a ejecutar en la lınea de co-mandos de Matlab. En ellos se puede encontrar la evolucion del error para todoslos nodos de la malla bajo la influencia del parametro de malla en un problemano lineal. Tambien se aportan soluciones para el problema concreto de PCF enel proceso de calentamiento del modelo no lineal mencionado en la memoria delproyecto. Asimismo, se adjuntan los resultados numericos del problema al variarlas condiciones externas (temperatura del horno, dimensiones preforma) como laspropiedades termicas (densidad, calor especifico, conductividad) del material de lafibra.

Carpeta CODIGOS: Subdivida en tres subcarpetas llamadas Exacta, PDE lineal yPDE no lineal. En ellas se hallan los codigos que resuelven los distintos tipos deproblemas tratados. La carpeta PDE lineal contiene dos subcarpetas FVMB y FVMC

con los codigos numericos de los volumenes finitos baricentrico y circuncentrico,respectivamente, implementados para la resolucion de problemas lineales. Se anadeun breve manual para el uso de las funciones implementadas.

75

APENDICE A 76

A.2. Codigos implementados

En esta seccion se describen brevemente los codigos implementados. Distinguimosentre las tecnicas utilizadas para resolver de forma analıtica y numerica los problemasregidos por PDEs lineales, y las tecnicas numericas para los problemas regidos por PDESno lineales.

PDEs lineales

Para el caso concreto del problema de valores iniciales dado por la ecuacion ∂u∂t−∆u =

0, en un dominio circular con condiciones de contorno Dirichlet, se puede obtener lasolucion analıtica mediante transformaciones de Laplace. La solucion expresada en formade un sumatorio de infinitos terminos requiere del computo de un determinado numero deceros de la funcion de Bessel de primera especie y orden nulo. Para ello se ha implementadoel programa solucion_exacta.m que devuelve la solucion analıtica. Este mismo problemase puede resolver por FEM haciendo uso de la funcion parabolic.m de la PDEToolboxde Matlab.

Para este sencillo problema lineal y otros mas complicados como el expresado por laecuacion (2.20) se aportan las funciones necesarias para resolverlos numericamente en lascarpetas FVMB y FVMC. Su uso ya fue explicado en la memoria.

PDEs no lineales

El problema del calentamiento de preformas de fibras de cristal fotonico estudia-do en este proyecto ha sido resuelto mediante las funciones incluidas en la carpetaPDE no lineal. Para la resolucion de este problema se hace una llamada a la funcionmatrizT1.m que nos devuelve de forma estructurada toda la informacion necesaria sobrela malla. Posteriormente, se hace uso de la funcion solucion_no_lineal.m encargada deplantear el sistema matricial iterativo y de devolver la solucion para todos los nodos de lamalla y en todos los instantes de tiempo. Los coeficientes de la PDE que representan laspropiedades termicas del material son no lineales, su valor y el de sus derivadas respectode la temperatura se calcula mediante funciones externas.

Apendice B

Codigos en Matlab

A continuacion se muestran algunos codigos sencillos que han sido mencionados en lamemoria. La totalidad de los codigos utilizados en la realizacion del proyecto se adjuntanen un CD-Rom, como es el caso de la implementacion del FVM circuncentrico pararesolver problemas PDEs no lineales .

B.1. Distancia mınima

En el desarrollo del metodo de volumenes finitos se requiere del calculo de la distanciamınima entre un nodo, i, dado (p.e. circuncentro de un VC) y una recta dada por dospuntos, 1 y 2, que definen una arista frontera). El codigo siguiente nos devuelve estadistancia.

distanciaminima.m

function ljs1=distanciaminima(xi,yi,x1,y1,x2,y2);

1 imagi=sqrt(-1);

2 vec_12=(x2-x1)+(y2-y1)*imagi;

3 vec_1i=(xi-x1)+(yi-y1)*imagi;

4 alfa=angle(vec_1i)-angle(vec_12);

5 ljs1=abs(abs(vec_1i)*sin(alfa));

Fin distanciaminima.m

77

APENDICE B 78

B.2. Circuncentro

Existen programas estandar que generan la posicion del circuncentro de un triangulodadas las coordenadas de los vertices. Sin embargo, se esta interesado en desarrollar porcompleto el metodo numerico lo que nos lleva a proponer una nueva version que halla lalocalizacion del circuncentro calculando la interseccion de las mediatrices de las aristas.El codigo es el siguiente.

Circuncentro.m

function [Xcirc,Ycirc]=circuncentro(x1,y1,x2,y2,x3,y3);

1 xma=(x1+x2)/2;

2 yma=(y1+y2)/2;

3 xmb=(x1+x3)/2;

4 ymb=(y1+y3)/2;

5 a=’’; b=’’;

6 if (y2-y1)==0

7 Xcirc=xma;

8 a=1;

9 else

10 ma=(x2-x1)/-(y2-y1);

11 pa=yma-ma*xma;

12 end

13 if (y3-y1)==0

14 Xcirc=xmb;

15 b=1;

16 else

17 mb=(x3-x1)/-(y3-y1);

18 pb=ymb-mb*xmb;

19 end

20

21 if ~isempty(a)

22 Ycirc=mb*Xcirc+pb;

23 elseif ~isempty(b)

24 Ycirc=ma*Xcirc+pa;

25 else

26 Xcirc=(pb-pa)/(ma-mb);

27 Ycirc=ma*Xcirc+pa;

28 end

Fin Circuncentro.m

APENDICE B 79

B.3. Sistema EDO baricentro

Con la informacion generada mediante la funcion MatrizT2 se puede proceder a cons-truir el sistema EDO que resuelve el problema gobernado por la PDE parabolica delcapıtulo 3, mediante el FVM baricentrico. El codigo es el siguiente.

sistemaABC2.m

function [A,B,C]=sistemaABC2(T2,a,b,c,d,uext);



1 for i=1:length(t1)

2 if t1(14,i)==0

Triangulos interiores, ecuacion (2.14)3 B(i,t1(15,i))=b*T2(8,i)*T2(20,i)/T2(11,i);

4 B(i,t1(16,i))=b*T2(9,i)*T2(21,i)/T2(12,i);

5 B(i,t1(17,i))=b*T2(10,i)*T2(22,i)/T2(13,i);

6 B(i,i)=c*T2(5,i)-b*(T2(8,i)*T2(20,i)/T2(11,i)+T2(9,i)*...

T2(21,i)/T2(12,i)+T2(10,i)*T2(22,i)/T2(13,i));

7 C(i,1)=d*T2(5,i);

Triangulos con una arista perteneciente al contorno, ecuacion (2.15)


9 ljs1=T2(18,i)

10 if T2(11,i)==0

11 C(i)=d*T2(5,1)-uext*b*T2(8,i)/ljs1;

12 B(i,T2(16,i))=b*T2(9,i)*T2(21,i)/T2(12,i)+c*T2(5,i);

13 B(i,T2(17,i))=b*T2(10,i)*T2(22,i)/T2(13,i)+c*T2(5,i);

14 B(i,i)=c*T2(5,i)-b*(T2(9,i)*T2(21,i)/T2(12,i)+...

T2(10,i)*T2(22,i)/T2(13,i)+T2(8,i)/ljs1);

15 elseif T2(12,i)==0

16 C(i)=d*T2(5,1)-uext*b*T2(9,i)/ljs1;

17 B(i,T2(15,i))=b*T2(8,i)*T2(20,i)/T2(11,i)+c*T2(5,i);

18 B(i,T2(17,i))=b*T2(10,i)*T2(22,i)/T2(13,i)+c*T2(5,i);

19 B(i,i)=c*T2(5,i)-b*(T2(8,i)*T2(20,i)/T2(11,i)+...

T2(10,i)*T2(22,i)/T2(13,i)+T2(9,i)/ljs1);

20 elseif t1(13,i)==0

21 C(i)=d*T2(5,1)-uext*b*T2(10,i)/ljs1;

22 B(i,T2(15,i))=b*T2(8,i)*T2(20,i)/T2(11,i)+c*T2(5,i);

Sigue . . .

APENDICE B 80

sistemaABC2.m (cont.)

23 B(i,T2(16,i))=b*T2(9,i)*T2(21,i)/T2(12,i)+c*T2(5,i);

24 B(i,i)=c*T2(5,i)-b*(T2(8,i)*T2(20,i)/T2(11,i)+...

T2(9,i)*T2(21,i)/T2(12,i)+T2(10,i)/ljs1);

25 end

Triangulos con dos aristas perteneciente al contorno, ecuacion (2.16)

26 elseif T2(14,i)==2

27 ljs1=T2(18,i)

28 ljs2=T2(19,i)

29 if T2(11,i)~=0


31 B(i,T2(15,i))=b*T2(8,i)*T2(20,i)/T2(11,i);

32 B(i,i)=c*T2(5,i)-b*(T2(9,i)/ljs1+T2(10,i)/ljs2+T2(8,i)*...

T2(20,i)/T2(11,i));

33 elseif T2(12,i)~=0


35 B(i,T2(16,i))=b*T2(9,i)*T2(21,i)/T2(12,i);


T2(21,i)/T2(12,i));

37 elseif T2(13,i)~=0


39 B(i,T2(17,i))=b*T2(10,i)*T2(22,i)/T2(13,i);


T2(22,i)/T2(13,i));

41 end

42 end

43 A(i,i)=a*t1(5,i);

44 end

Fin sistemaABC2.m

Apendice C

Calidad del mallado

La influencia de la calidad de la malla en las tecnicas numericas es crucial para laobtencion de resultados de calidad. Es de especial importancia el hacer un estudio de lamalla a utilizar previamente a evaluar el metodo. Con el termino alta calidad, mencionadoen los apartados 2.3.2 y apartado 2.6, nos referimos a la relacion entre longitudes de aristasdel triangulo, cuanto mas cercana sea de la unidad mejor calidad se tiene. Por tanto, setrata de conseguir triangulos lo mas equilaterales posible. El metodo de Delaunay usadopara la construccion del mallado asegura una alta calidad de los triangulos que lo formany por tanto una buena disposicion centrada de los nodos de calculo localizados en loscircuncentros. Sin embargo en ciertos casos se tienen triangulos que no cumplen estosrequisitos, en particular los triangulos obtusos con los que el circuncentro se halla fueradel triangulo en cuestion, figura 2.1.

(a) (b)

Figura C.1: Calidad de los triangulos (a) de una malla de Delaunay, y (b) simetrica

Por calidad de triangulos se entiende la relacion q = 4a√

3/(h21 + h2

2 + h23), siendo a el

area del triangulo y h1, h2, h3 las longitudes de sus aristas [57]. Los triangulos equilaterales

81

APENDICE C 82

tienen la propiedad de tener por calidad la unidad y son los idoneos para el calculonumerico. Sin embargo, los dominios de calculo complejos no permiten el uso de este tipode triangulos.

En el capıtulo 2, se hizo un estudio del comportamiento del FVM baricentrico endistintas mallas. Los resultados fueron poco satisfactorios al utilizar mallas simetricasobtenidas a partir de refinamientos sucesivos. Para mallas de Delaunay los resultadosmejoran ostensiblemente. Este comportamiento es debido a la calidad de ambas mallas.La calidad de los dos tipos de malla se muestra en la figura C.1. Se observa que paramallas de Delaunay existen un determinado numero de triangulos con calidad muy pordebajo de la calidad mınima de los triangulos de la malla simetrica. No obstante, lacalidad media de la malla de Delaunay es superior a la de la malla simetrica. Por tanto,es mas importante disponer de una malla con calidad media de triangulos elevada aunconteniendo ciertos triangulos de mala calidad.


Apendice D

Artıculo enviado para publicacion

Los resultados obtenidos durante el desarrollo de este proyecto fin de carrera hanpermitido elaborar el artıculo, que se adjunta seguidamente, enviado a la revista interna-cional Computer Physics Communications (ISI IF 1.5):

Pablo Udias de la Mora, Francisco R. Villatoro and Jose F. Aldana, “BarycentricVersus Circumcentric Cell-Centered Finite Volume Methods”.

Los siguientes 3 artıculos de investigacion estan en preparacion:

1. Pablo Udias de la Mora, Francisco R. Villatoro, and Jose F. Aldana, “Finite VolumeMethods Using Matlab,” to be submitted to Computers Chemical Engineering, Inpreparation.

2. Pablo Udias de la Mora, Francisco R. Villatoro, and Jose F. Aldana, “NonlinearHeat Transfer in Polymer Photonic Crystal Fibre Preforms,” to be submitted toModelling and Simulation in Materials Science and Engineering, In preparation.

3. Pablo Udias de la Mora, Francisco R. Villatoro, and Jose F. Aldana, “Finite VolumeMethods by Quasilinearization for Nonlinear Heat Transfer,” to be submitted toComputers and Mathematics with Applications, In preparation.


83

Barycentric Versus Circumcentric

Cell-Centered Finite Volume Methods

Pablo Udias de la Mora a, Francisco R. Villatoro a,∗, andJose F. Aldana b

aDepartamento de Lenguajes y Ciencias de la Computacion, E.T.S. IngenierosIndustriales, Universidad de Malaga, Plaza El Ejido, s/n, 29013-Malaga, SpainbDepartamento de Lenguajes y Ciencias de la Computacion, E.T.S. IngenierıaInformatica, Universidad de Malaga, Campus de Teatinos, s/n, 29071-Malaga,

Spain

Abstract

Cell-centered finite volume methods on triangular meshes based on piecewise con-stant functions over each triangle where the solution is located at either the barycen-ter or the circumcenter are compared between them and with a finite elementmethod for the solution of linear parabolic equations. The circumcenter of theboundary triangles must be inside them, in other case the solution presents spuriousoscillations. A method to attain such a procedure is presented. Circumcentric finitevolume methods with piecewise constant approximations of the unknown functionare as accurate as finite element methods with piecewise approximations, with theadditional advantage of presenting a discrete maximum principle which avoids unre-alistic solutions. Barycentric finite volume methods are less accurate and less robustto mesh parameters than circumcentric ones for irregular meshes obtained by con-strained Delaunay triangulation, so the last ones are recommended in applications.

Key words: Cell-centered finite volume methods, constrained Delaunay mesh,triangular mesh, barycenter, circumcenterPACS: 02.70.Bf, 47.11.-j

∗ Corresponding author.Email addresses: [email protected] (Pablo Udias de la Mora),

[email protected] (Francisco R. Villatoro), [email protected] (Jose F. Aldana).

Preprint submitted to Elsevier Science 12 May 2006

1 Introduction

Finite volume methods (FVM) are the most preferred numerical methodsamong physicists and engineers for the solution of thermal and fluid physicalproblems modelled with conservation laws represented by partial differentialequations (PDE) [1–3]. These methods allow the use of irregular meshes andare based on the integration of the PDE on certain control volumes, whicheither may belong to the mesh, in the so-called cell-centered FVMs [4,5], ormay belong to a dual mesh, in the vertex-centered FVMs [6,7], and then nu-merically approximate the conservative fluxes on the boundary of these con-trol volumes. In triangular meshes with the Delaunay property, cell-centeredmethods use the triangles as control volumes and vertex-centered methodsuse Voronoi polygons [8,9]. By the Euler formula, the number of triangles isalways larger than the number of vertices, therefore vertex-centered FVM re-sults in algebraic systems of equations with a smaller number of unknownsthan cell-centered FVMs, so they are preferred in the majority of applica-tions. However, the number of non-null elements in each row of the resultingmatrices is larger than that of the cell-centered FVM, yielding sparse matri-ces with best banded structure. The theoretical analysis [3] of both cell- andvertex-centered FVMs shows that both approaches are convergent with simi-lar theoretical error estimates, so the choice between them is a matter of thepractitioner’s preference, some ones using cell-centered FVM [4,5] and othersusing vertex-centered FVM [6,7].

This paper is focused on the cell-centered FVMs on Delaunay triangularmeshes using piecewise constant approximations of the unknown function.Several approaches are possible depending on the position of the point in-side each triangle where the solution is assumed to be located. The two mainpossibilities are the use of either the circumcenter [10] or the barycenter [11].The circumcentric FVM has the advantages that the fluxes are collinear tothe normal vector of each triangle side, which simplifies the computationalimplementation, but the disadvantage that the circumcenter is outside theobtuse triangles, which may appear in a Delaunay triangulation [12,13]. Thebarycentric FVM has the advantage that the barycenter is always located in-side the triangle, but the fluxes has to be projected in the normal vector ofthe triangle sides, so a cosine is required increasing the computational cost ofthe implementation.

To the authors’s knowledge there is no technical comparison between circum-centric and barycentric FVMs techniques in a practical context, this papertries to fill this gap, being its contents as follows. The next section presentsthe cell-centered finite volume methods in triangular meshes based on boththe circumcenter and the barycenter for a linear parabolic partial differentialequation with the trapezoidal rule for the integration in time. Section 3 shows

2

the main results of the comparison of the methods developed in this paper.Finally, the last section summarizes the main conclusions.

2 Problem formulation

Let us consider the initial-boundary value problem

α∂U

∂t+∇ · (β∇U) + γ U = q, t > 0, ~x ∈ Ω ⊂ R2, (1)

U(~x, 0) = u0(~x), ~x ∈ Ω, (2)

U(~x, t) = uext, t > 0, ~x ∈ ∂Ω, (3)

where U ≡ U(~x, t) is a scalar function, ∇ is the gradient operator, α, β, γ, q,and uext are constant, and Ω = Ω ∪ ∂Ω. Let us also assume that the Dirichletboundary condition is compatible with the initial one (u0(~x) = uext, ~x ∈ ∂Ω).

The linear parabolic equation (1) models, for example, heat conduction in abidimensional, homogeneous and isotropic medium when U(~x, t) is tempera-ture at position ~x and time t, α = ρ c, ρ is the density, c is the specific heat,−β is the thermal conductivity, γ is the convective heat transfer coefficient, qis a heat source, u0 is the initial distribution of temperature in the domain,and uext is a constant temperature at its boundary.

Equations (1)–(3) may be solved by means of cell-centered FVM using trian-gular meshes. Let us assume that the boundary of Ω is polygonal and thata constrained Delaunay triangulation TΩ of the domain has been obtained.This triangulation may be decomposed into TΩ = T∂Ω ∪TΩ, where T∂Ω are thetriangles with at least one side on the boundary, referred to as the boundarytriangles, and TΩ are the rest ones, referred to as the inner triangles.

Let us take as control volume a triangle Ri ∈ TΩ and integrate Eq. (1) in spaceover this control volume yielding

∫

Ri

(α∂U

∂t+∇ · (β∇U) + γ U − q) d~x = 0, ∀ Ri ∈ TΩ, (4)

and the divergence theorem may be applied giving

∫

Ri

(α∂U

∂t+ γ U − q) d~x +

∮

∂Ri

β∇U · ~n d~s = 0, (5)

3

where ∂Ri represents the three sides of triangle Ri, ~n is the outward-pointingunit normal to ∂Ri, and d~s is a triangle side differential element.

The solution U of Eq. (5) may be numerically approximated by a piecewiseconstant function u over the mesh TΩ, whose constant value ui over eachtriangle Ri is assumed to be located at its “center” ~xi; both the barycenterand the circumcenter will be used in this paper as triangle centers. Let Ni

be the index set of the neighboring triangles of Ri, those sharing a side withit. This set has three elements for inner triangles, but only either one or twofor boundary ones. Let Si be the area of the triangle Ri, eij be the length ofthe side shared by triangles Ri and Rj, with j ∈ Ni, dij = ‖~xj − ~xi‖ be thedistance between their centers, and αij be the angle between the vector ~xj−~xi

and the outward normal of the side shared by Ri and Rj. Let Ni be the indexset of the boundary sides of a triangle Ri ∈ T∂Ω, i.e., such that the cardinality#Ni < 3, and δik be the distance between ~xi and the side k ∈ Ni, i.e.,the distance between ~xi and the intersection point of the line passing through~xi perpendicular to the side of the triangle. Let ek be the length of the sidek ∈ Ni.

The integral over Ri in Eq. (5) may be evaluated by means of the rectanglequadrature rule and the integral over ∂Ri by means of a projection of thegradient ∇u into the normal vector ~n of the three sides of the triangle, wherethe gradient is approximated by a central finite difference formula, yielding

(αdui

dt+ γ ui − q) Si +

∑

j∈Ni

β

(uj − ui

dij

)cos(αij) eij = 0, (6)

for Ri ∈ TΩ, and

(αdui

dt+ γ ui − q) Si +

∑

j∈Ni

β

(uj − ui

dij

)cos(αij) eij

− ∑

k∈Ni

β(

ui

δik

)ek = − ∑

k∈Ni

β(

uext

δik

)ek, (7)

for Ri ∈ T∂Ω, where the Dirichlet boundary conditions have been introduced.

2.1 Circumcentric cell-centered finite volume method

The circumcenter is the center of a triangle’s circumcircle, i.e., a point suchthat its distance to the three vertices of a triangle is constant. It can be foundas the intersection of the perpendicular bisectors of its sides, as shown inFig. 1 (a). The use of the circumcenter as the “center” ~xi of a control volume

4

Ri ∈ TΩ makes the evaluation of Eqs. (6) and (7) cheapest since then all thecosines cos(αij) are exactly unity.

Circumcentric finite volume methods (C-FVM) have two drawbacks. In theone side, a division by zero occurs in either Eq. (6) or Eq. (7) when, respec-tively, either the circumcenter of two neighboring triangles coincides, or thecircumcenter for a boundary triangle is located in the boundary (this occurs forright boundary triangles). In the other side, in a constrained Delaunay trian-gulation, such as that used in this paper, obtuse triangles may appear. In sucha case, the circumcenter is located outside the corresponding triangle. For theinner triangles, this a small source of error since for a triangulation having theDelaunay property, i.e., if the circumcircle of any mesh triangle does not con-tain any other triangle vertex, the circumcenters are well distributed over thedomain, “representing” the circumcircles covering the triangles of the mesh.However, for the boundary triangles, the resulting error is very high, showingspurious oscillations, to be illustrated in Section 3. Fortunately, a constrainedDelaunay triangulation assures high quality triangles with a high degree ofequilaterally, hence only few triangles are (slightly) obtuse. Usually, none islocated in the boundary. Note that a Delaunay mesh without obtuse trianglesis an admissible mesh, as defined in Ref. [3], when using the circumcenters.

The two drawbacks of using the circumcenter as control volume center ina constrained Delaunay triangulation may be easily solved by “manual” re-touching of the mesh, since only a few exceptional triangles are conflicting.In such cases, the problem may be avoided by means of a displacement ofthe coordinates of one of the vertices of the conflicting triangle, which mustnot be located in the domain boundary, in order to relocate the correspond-ing circumcenter inside it. In the whole set of simulations summarized in thispaper where the mesh was generated by the Matlab PDE Toolbox, a vectordisplacement whose modulus is one-tenth of the smallest side of the triangleand whose direction is the mean of the unit normal vectors of the sides joiningin such a vertex have always solved the problem. Obviously, the resulting tri-angulation may loss the Delaunay property in a few triangles, but the qualityof the mesh is only slightly affected with such a small displacement.

2.2 Barycentric cell-centered finite volume method

The barycenter, also referred to as either the centroid or the center of gravity,of a triangle is the intersection of its three medians, the lines passing trougha vertex and the midpoint of the opposite side, being always inside of it, asshown in Fig. 1 (b).

The use of the barycenter as the “center” ~xi of a control volume Ri ∈ TΩ has

5

(a) (b)

Fig. 1. Illustration of the position of the circumcenter (a) and the barycenter (b) ofa triangle and their neighborhood ones, including the notation used in this paper.

the advantage of always being inside the triangle, but requires the explicit eval-uation of the cosines cos(αij) in Eqs. (6) and (7) increasing the computationalcost of the method. For high quality triangles, like those in constrained Delau-nay triangulation, the cosines are near unity, so some authors recommend thesubstitution of all the cosines by unity when implementing barycentric finitevolume methods (B-FVM).

2.3 Time integration

Equations (6) and (7) for, respectively, the inner and boundary triangles yielda system of ordinary differential equations for each control volume of the mesh,which may be assemble into matrix form as

Adu

dt+ B u = C, (8)

where N is the number of triangles in the mesh TΩ, u ∈ RN is a column vectorwith the values ui over all the triangles of the mesh, A is a N × N diagonalmatrix, B is a N ×N sparse matrix with four non-null elements in each rowfor the control volumes in TΩ, and two or three non-null elements depending ifthe control volume in T∂Ω has one or two sides in the boundary, respectively,and C is a null vector for the rows corresponding to triangles in TΩ, beingnon-null otherwise.

In this paper, the system of differential equations (8) is solved in time bymeans of the trapezoidal rule, which yields

un+1 =

(A +

4t

2B

)−1 (4t C +

(A− 4t

2B

)un

), (9)

6

where un corresponds to vector u at time tn, and the initial condition is u0 =u0.

The numerical method developed in this section is second-order accurate andunconditionally stable in time. However, physically realistic solutions are notobtained for arbitrary time steps, since solutions with high frequency compo-nents may develop spurious oscillations due to the trapezoidal rule is A stable,but not L stable [14,15]. In order to avoid these spurious oscillations the timestep may be properly selected [15]. Either being smaller than a critical onewhich depends on the highest frequency of the solution or being adaptiveand logarithmically distributed. In this paper this second solution is adoptedto discretize the time interval [0, T ] by means of t0 = 0, tn = 10(n−1) τ ∆t,τ = (log10(T )− log10(∆t))/(N − 1), and n = 1, 2, . . . , N .

3 Presentation of results

Let us consider the parabolic problem

∂u

∂t− 1

r

∂

∂r

(r

∂u

∂r

)= 0, 0 ≤ r < R, t > 0, (10)

with u(r, 0) = 0, 0 ≤ r < R, and u(R, t) = uext, whose analytical solution maybe obtained by means of the Laplace transform yielding [16, Chap. III, Eq. (7)]

u = uext − 2 uext

R

∞∑

n=1

J0(rαn)

αn J1(Rαn)e−α2

n t, (11)

where J0 and J1 are Bessel functions of the first kind and αn are the roots ofJ0(Rα) = 0. This analytical solution has been numerically evaluated by calcu-lating the roots αn with a relative error of 10−15 by using Halley’s method [17]and summing the series up to the relative error is smaller than 10−12.

The solution of the circumcentric and the barycentric FVM developed in thispaper has been compared with both the analytical solution and the numer-ical solution obtained by means of Matlab’s PDEToolbox which uses a fi-nite element method (FEM) based on linear piecewise functions. Figure 2shows the error between the analytical solution and the numerical solution,i.e., U(xi, t

n) − uni , for the FEM (dotted line), the circumcentric FVM (line

with + signs), the barycentric FVM using cosines (line with × signs) and thebarycentric FVM without cosines (line with signs) where the circular domainΩ of unit radius has been triangulated by means of Matlab’s PDEToolbox,which uses a constrained Delaunay method, taking the option that the largest

7

0 0.2 0.4 0.6 0.8 1−2

0

2

4

6

8

10

12x 10

−3

t

U(a

na

l.)−

u(n

um

.)

0 0.2 0.4 0.6 0.8 1−2

0

2

4

6

8

10

12x 10

−3

t

U(a

na

l.)−

u(n

um

.)

(a) (b)

Fig. 2. Difference between the analytical solution of Eq. (10) and the numericalone obtained by means of a FEM (dotted line), a circumcentric FVM (line with+ signs), a barycentric FVM using cosines (line with × signs), and a barycentricFVM without cosines (line with signs) for (a) the central node of the mesh (r = 0)and (b) a node located at r = 0.5.

side of a triangle of the mesh is 0.1, resulting in a mesh with 1058 trianglesand 562 vertices. Figure 2 (a) and (b) show the error at the centre of thecircle (r = 0) and at a point with r = R/2, respectively. In order to avoidpossible interpolation errors in the calculation of the solution at these points,the mesh has been manually fine tuned in order to locate the nearest vertex(in the FEM), or either the nearest circumcenter (in the circumcentric FVM)or the nearest barycenter (in the barycentric FVM) exactly in that position.Hence, the triangle mesh is not exactly the same in all the three methods, butonly differs in the position of a small number of triangles.

Figures 2 (a) and (b) show that the error for the FEM and the circumcen-tric FVM are very similar in magnitude being smaller than that for bothbarycentric FVMs. Surprisingly, the apparently more inaccurate barycentricFVM with all the cosines approximated by unity is, in fact, slightly more ac-curate. The evolution in time of the error shows the same trends for threemethods. The error initially grows in time until reaching a maximum, wherethe analytical solution is larger than the numerical one, and then decreasesasymptotically to zero, for both barycentric FVMs, and up to a negative min-imum, for both the FEM and the circumcentric FVM; in this last case, theerror grows again to reach a new positive maximum after which also decreasesasymptotically to zero. Note also that the FEM method is more efficient, incomputational time, than the FVMs since it requires the solution of a sparselinear system with only 562 unknowns instead of the 1058 ones required by theFVMs. However, the FEM do not assures the minimum-maximum principleof the solution and unrealistic negative values appear in the first time steps,completely avoided by the FVMs.

8

0.05 0.1 0.15 0.2

0

0.02

0.04

0.06

0.08

0.1

t

U(a

na

l.)−

u(n

um

.)

0 0.2 0.4 0.6 0.8 1

10−4

10−2

100

t

U(a

na

l.)−

u(n

um

.)

(a) (b)

Fig. 3. Maximum error between the analytical solution of Eq. (10) and the numericalone obtained by means of a FEM (dotted line), a circumcentric FVM (line with +signs), a barycentric FVM using cosines (line with × signs), and a barycentric FVMwithout cosines (line with signs) for t ∈ [0, 0.2] (a) and t ∈ [0, 1].

Figure 3 shows the evolution in time of the maximum error for the FEMand the three FVMs. The error is very large in the first time steps, as shownin Fig. 3 (a), being located in the boundary triangles and, therefore, due tothe inconsistency between the initial condition and the Dirichlet boundaryconditions. During this transient, the smallest error is that of the FEM, andthe largest of the barycentric FVM with cosines. Figure 3 (b) shows that theerror of the FEM and circumcentric FVM is similar in magnitude during thewhole simulation, but that of the barycentric FVM is one order of magnitudelarger, being practically the same for both with and without cosines. Afterreaching the steady state, near t = 1, the four methods yield similar maximumerror.

Figure 4 shows the relative error between the analytical solution and thesolution of the FEM (a), the circumcentric FVM (b), the barycentric FVMwith cosines (c) and without cosines (d) at time t = 0.5. The FEM andcircumcentric FVM methods yield a smaller error with also smaller dispersionthan that of the barycentric FVMs, both very similar among them. Since thesolution of Eq. (12) is radially symmetric, this indicates that the solution withthe barycentric FVMs is less symmetric (and hence worst) than that of theother two methods. This result may be unexpected since a more expensivemethod (B-FVM with cosines) yields worst results than an approximation(B-FVM without cosines) and even the other two methods.

Let us consider now the numerical solution of a problem without analyticalsolution given by the partial differential equation

100∂U

∂t−∆U + U = 10, t > 0, ~x ∈ Ω ⊂ R2, (12)

9

0 0.2 0.4 0.6 0.8 1−3

−2

−1

0

1x 10

−3

t

U(a

na

l.)−

u(F

EM

)

0 0.2 0.4 0.6 0.8 1−3

−2

−1

0

1x 10

−3

t

U(a

na

l.)−

u(C

−F

VM

)

(a) (b)

0 0.2 0.4 0.6 0.8 1−3

−2.5

−2

−1.5

−1

−0.5

0

0.5x 10

−3

t

U(a

na

l.)−

u(B

−F

VM

)

0 0.2 0.4 0.6 0.8 1−3

−2.5

−2

−1.5

−1

−0.5

0

0.5x 10

−3

t

U(a

na

l.)−

u(B

−F

VM

with

ou

t)

(c) (d)

Fig. 4. Relative error between the analytical solution and the solution of theFEM (a), the circumcentric FVM (b), the barycentric FVM using cosines (c), andthe barycentric FVM without cosines (d) at time t = 0.5.

with Dirichlet boundary conditions U(~x, t) = 1, for ~x ∈ ∂Ω and t > 0, andinitial condition U(~x, 0) = 0, for ~x ∈ Ω. Let us use exactly the same Delaunaymesh that in the previous problem with 1058 triangles. Figure 5 (a) shows azoom of part of the mesh illustrating the location of both the circumcentersand the barycenters. In this mesh, the circumcenters of 24 triangles are outsidethem, representing the 2.3% of the total. Among them, only one boundarytriangle has its circumcenter outside the domain, indicated by C1 in Fig. 5 (a).The mesh may be manually retouched as indicated in Section 2.1 in order torelocate this circumcenter inside the corresponding boundary triangle, butwith the possible loss of the Delaunay property in adjacent triangles. Theresulting mesh is shown in Fig. 5 (b), where B′

1 and C ′1 corresponds to the

new centers of the triangle having, respectively, B1 and C1 in Fig. 5 (a).

Figures 6 (a) and 6 (b) show the solution of Eq. (12) at several nodes of themesh in Fig. 5 (a) by means of cell-centered FVM using, respectively, thecircumcenters and the barycenters. The most noticeable aspect in Fig. 6 (a)

10

(a) (b)

Fig. 5. Zooming of a part of a constrained Delaunay mesh for a circular domainshowing the location of the circumcenters (+ signs) and barycenters (× signs) ofthe triangles, as originally generated by Matlab PDEToolbox (a) and after manualretouching, with the circumcenter of the boundary triangles inside them.

is that the solution at the circumcenter of the boundary triangle lying out-side the domain is not monotonic, with a big error for the first time stepsand an oscillatory behavior with a frequency diminishing as time progresses.Outside this region the solution evolves always monotonically, illustrated bythat located at C2 in Fig. 6 (a). The results for the barycentric FVM shownFig. 6 (b) are always monotonic. The largest difference between the solutionsof the circumcentric and barycentric FVM is smaller than the 0.15%. A man-ual retouching of the conflicting triangle is easily accomplished as indicated inSection 2.1, completely solving the problem with the oscillations of the solu-tion in such a triangle, yielding a monotonic result which is indistinguishableto that of Fig. 6 (b) at that plot resolution. This also happens for the solutionwith both barycentric FVMs, with and without cosines.

Figure 7 (a) shows a three-dimensional representation of the solution of thecircumcentric FVM, with the mesh in the bottom of the plot, for all thecircumcenters of mesh at time t = 10. Figure 7 (b) shows the solution over theplane π in Fig. 7 (a) for several times. This figure illustrates how the solutionevolves in time, shown that the solution at the center of domain grows in timeand that for intermediate times, cf. t = 10, the largest value of the solutionis neither in the boundary nor in the center. The final steady state solutionis practically reached at t = 100. Note that the corresponding figures for thebarycentric FVM are practically indistinguishable at the graphical resolutionused in the plots.

Let us study the influence of the quality of the mesh in the numerical results,considering two Delaunay meshes, one symmetrical, Fig. 8 (a), and the otherone asymmetrical, Fig. 8 (b), the last one obtained by the Matlab jiggle

function. There are several metrics for triangle quality, but the most widely

11

0 50 1000

1

2

3

t

u

0 50 1000

1

2

3

t

u

(a) (b)

Fig. 6. (a) Numerical solution with the circumcenter FVM for nodes C1 (continuousline) and C2 (dashed line) in Fig. 5 (a), a circumcenter at r = 0.5 (dotted line)and another one at r=0 (dashed-dotted line), and (b) Numerical solution with thebarycenter FVM for nodes B1 (continuous line) and B2 (dashed line) in Fig. 5 (b),a barycenter at r = 0.5 (dotted line) and another one at r=0 (dashed-dotted line).

0 0.2 0.4 0.6 0.8 10

0.5

1

1.5

2

2.5

3

r

u(C

−F

VM

)

π

(a) (b)

Fig. 7. (a) Three-dimensional representation of the solution of the circumcentricFVM for the mesh in Fig. 5 (a) showing the values of the solution at all the cir-cumcenters in the mesh at time t = 100, and (b) the solution at plane π for t = 2,5, 10, and 100.

used is q = 4 a√

3/(l21 + l22 + l23), where a is the triangle area, and l1, l2, andl3 are the side lengths of the triangle [18]. For an equilateral triangle q = 1.The triangle quality of the two meshes in Fig. 8 are indicated in the plots bymeans of a gray scale. The mean triangle quality of the asymmetrical mesh, seeFig. 8 (b), is better than that of the symmetrical one, see Fig. 8 (a). However,the worst quality triangles appear on the asymmetrical mesh.

The relative error between the analytical solution and the barycentric FVMwith cosines (×-signed points) and the circumcentric FVM (+-signed points)at steady state (t = 100) are shown in Fig. 9 (a) and (b) for, respectively,

12

0.88

0.9

0.92

0.94

0.75

0.8

0.85

0.9

0.95

(a) (b)

Fig. 8. Symmetrical (a) and asymmetrical (b) triangular constrained Delaunaymeshes indicating the triangle quality using a gray scale.

0 0.2 0.4 0.6 0.8 1−0.02

0

0.02

0.04

0.06

0.08

0.1

r

Er

0 0.2 0.4 0.6 0.8 1−0.015

−0.01

−0.005

0

0.005

0.01

0.015

r

Er

(a) (b)

Fig. 9. Relative error between the finite element solution and the barycentric FVM(×-signed points) and the circumcentric FVM (+-signed points) at t = 100 for thesymmetrical (a) and the asymmetrical (b) Delaunay meshes shown in Fig. 8 (a)and (b), respectively.

the asymmetrical and the symmetrical Delaunay meshes shown in Fig. 8 (a)and (b). The error for the circumcentric FVM is an order of magnitude smallerthan that for the barycentric FVM in both meshes. The symmetry of themesh only slightly affects the error for the circumcentric FVM. However, thebarycentric FVM solution on the symmetrical mesh has a larger error than thatof the asymmetrical one. Therefore, the symmetry of the mesh has a significanteffect on the barycentric FVM and only slightly affects the circumcentric FVM,showing is best robustness on changes on the mesh.

The global error between both the circumcentric FVM and the barycentricFVM (using cosines) with respect to the FEM solution for both asymmetricaland symmetrical Delaunay meshes is shown in Table 1, which also indicatesthe number of triangles, the mesh size parameter, and the mean and minimumtriangle quality of the mesh. For the circumcentric FVM, Table 1 shows that

13

Mesh Quality Method Max (%) Mean (%) STD

Asymmetrical qmean=0.95 C-FVM 2.58 1.73 0.19

(120 t.,h=0.3) qmin=0.78 B-FVM 8.13 1.47 0.23


(490 t.,h=0.15) qmin=0.73 B-FVM 3.56 0.44 0.14


(2184 t.,h=0.07) qmin=0.77 B-FVM 2.03 0.03 0.18

Symmetrical qmean=0.90 C-FVM 2.56 1.85 0.21

(128 t., h=0.3) qmin=0.87 B-FVM 9.74 6.70 0.8


(512 t., h=0.155) qmin=0.86 B-FVM 7.45 4.86 1.21


(2048 t., h=0.079) qmin=0.85 B-FVM 7.64 4.40 2.22

Table 1Relative error between the circumcentric FVM and barycentric FVM using cosineswith respect to the FEM for asymmetrical and symmetrical meshes, indicating thenumber of triangles and the mesh size (first column), the mean and minimum tri-angle quality of the mesh (second column), the corresponding FVM (third column),and the maximum, mean, and standard deviation of the relative error (last threecolumns).

the error diminishes, as expected, when the mesh size is reduced (the numberof triangles grows), being slightly better for the symmetrical mesh than forthe asymmetrical one due to the error in low quality triangles is large and theminimum triangle quality in such a case is the smallest. Moreover, the standarddeviation of the error is smaller for asymmetrical meshes than for symmetricalones, due to the fact that the standard deviation in triangle quality around itsmean is also small. For the barycentric FVM, Table 1 shows that the influenceof the low quality triangles is higher than for the circumcentric FVM, beingthe error smaller for the asymmetrical mesh than for the symmetrical one.The error of the barycentric FVM for asymmetrical meshes is smaller thanthat of the circumcentric FVM. In general, the standard deviation of the errorfor the barycentric FVM is larger than that for the circumcentric FVM.

14

4 Conclusions

A parabolic partial differential equation has been solved by means of cell-centered FVMs on a circular domain triangulated using a constrained Delau-nay mesh with either the circumcenter or the barycenter of the triangles ascenters of the control volumes and a piecewise constant approximation for thesolution. The methods were compared among them and with a FEM based onpiecewise linear elements. The circumcentric FVM is as accurate as the FEM,being more accurate than both barycentric FVMs, with and without cosinesin the formulation, whose error have the largest standard deviation. The prob-lem with the FEM is the lack of discrete maximum principle which may yieldunrealistic solutions for the first time steps, for example, with negative so-lutions. The circumcentric FVM has the problem that for obtuse boundarytriangles the circumcenter is outside the domain yielding a spurious oscillatorysolution, whose effect is local affecting the rest of the solution very small. Asolution for this problem has been proposed based on retouching of the mesh.Barycentric FVM using the projection of the fluxes on the normal vectors ofthe triangle sides, introducing some cosines in the formulation, and barycen-tric FVM approximating the cosines by unity have been compared showingthat, surprisingly, this second approach yields slightly more accurate results.The effect of the symmetry and triangle quality of the mesh were also stud-ied, showing that the presence of low quality triangles only slightly affects theerror for the circumcentric FVM, but instead the barycentric FVM solution islargely affected, being the best results obtained for symmetrical meshes andwhen the lowest triangle quality is high. In this sense, the circumcentric FVMis more robust on changes on the mesh than the barycentric FVMs

In summary, our results indicate that circumcentric finite volume methodswith a retouched mesh in order to avoid circumcenters outside the domainmust be preferred to barycentric finite volume methods, which always mustbe implemented approximating the cosines by unity.

The comparison of circumcentric and barycentric cell-centered methods onnonlinear partial differential equations and advection dominated problems areinteresting topics for further study. Moreover, the advantages introduced bythe use of new Delaunay refinement algorithms assuring the absence of obtusetriangles in the mesh [19] must also be studied.

Acknowledgments

The research reported in this paper was partially supported by Projects FIS2005-03191 and TIN2005-09098-C05-01 from the Ministerio de Educacion y Ciencia,

15

Spain.

References

[1] R.J. LeVeque, Finite Volume Methods for Hyperbolic Problems, CambridgeUniversity Press, Cambdridge, UK (2002).

[2] C. Hirsch, Numerical Computation of Internal and External Flows, Vol. 1,Wiley (1988).

[3] R. Eymard, T. Gallouet, and R. Herbin, “Finite volume methods,” in Handbookof Numerical Analysis, edited by P.G. Ciarlet and J.L. Lions, North Holland,pp. 715–1022 (2000).

[4] E. Bertolazzi and G. Manzini, “A cell-centered second-order accurate finitevolume method for convection-diffusion problems on unstructured meshes,”Mathematical Models and Methods in Applied Sciences, Vol. 14, pp. 1235–1260(2004).

[5] R.E. Ewing, O. Sævareid, and J. Shen, “Discretization schemes on triangulargrids,” Computer Methods in Applied Mechanics and Engineering, Vol. 152, pp.219–238 (1998).

[6] M. Ohlberger, “A posteriori error estimates for vertex centered finite volumeapproximations of convection-diffusion-reaction equations,” MathematicalModelling and Numerical Analysis, Vol. 35, pp. 355–387 (2001).

[7] P.A. Jayantha and I.W. Turner, “A comparison of gradient approximationsfor use in finite-volume computational models for two-dimensional diffusionequations,” Numerical Heat Transfer, Part B, Vol. 40, pp. 367–390 (2000).

[8] Q. Du and L. Ju, “Numerical simulations of the quantized vortices on a thinsuperconducting hollow sphere,” Journal of Computational Physics, Vol. 201,pp. 511–530 (2004).

[9] Q. Du and D. Wang, “The optimal centroidal Voronoi tessellations andthe Gersho’s conjecture in the three-dimensional space,” Computers andMathematics with Applications, Vol. 49, pp. 1355–1373 (2005).

[10] R. Herbin and O. Labergerie, “Finite volume schemes for elliptic and elliptic-hyperbolic problems on triangular meshes,” Computational Methods in AppliedMechanics and Engineering, Vol. 147, pp. 85–103 (1997).

[11] S. Chou and S. Tang, “Comparing Two Approaches of Analyzing Mixed FiniteVolume Methods,” Journal of the KSIAM, Vol. 5, pp. 55–78 (2001).

[12] F. Aurenhammer and R. Klein, Handbook of Computational Geometry, North-Holland, Amsterdam (2000).

16

[13] Z.C. Li and S. Wang, “The finite volume method and application incombinations,” Journal of Computational and Applied Mathematics , Vol. 106,pp. 21–53 (1999).

[14] S.V. Patankar, Numerical Heat Transfer and Fluid Flow, Taylor & Francis(1980).

[15] W.L. Wood, Practical Time-stepping Schemes, Clarendon Press, Oxford (1990).

[16] H.C. Carslaw and J.C. Jaeger, Conduction of Heat in Solids, Clarendon Press,Oxford (1986).

[17] D. Kincaid and W. Cheney, Numerical Analysis. Mathematics of ScientificComputing, (3rd. ed.) Brooks/Cole, Belmont (2002).

[18] D.A. Field, “Qualitative measures for initial meshes,” International Journal forNumerical Methods in Engineering, Vol. 47, pp. 887–906 (2000).

[19] J.R. Shewchuk, “Delaunay refinement algorithms for triangular meshgeneration,” Computational Geometry, Vol. 22, pp. 21–74 (2002).

17

pdetool-440

Documents

Transcript of pdetool-440