Explicit time integration with large time steps for the heat flow equations

Resumen

Este trabajo presenta un nuevo método de integración temporal puramente explícito que es estable para grandes pasos de tiempo. El mismo no requiere inversión de matrices ni resolución de un sistema de ecuaciones y por lo tanto es apto para ser paralelizado fácilmente y con alta eficiencia. En este primer trabajo se aplicará la teoría propuesta para resolver problemas de transmisión del calor por conducción, demostrando que es estable para cualquier paso de tiempo como ocurre con los métodos implícitos pero con un costo computacional mucho menor.

Abstract

We present a novel fully explicit time integration method that remains stable for large time steps, requires neither matrix inversions nor solving a system of equations and therefore allows for nearly effort-less parallelization. In this first paper the proposed approach is applied to solve conduction heat transfer problems, showing that it is stable for any time step as is the case with implicit methods but with a much lower computation time.

Palabras clave

Elementos finitos ; Método explícito ; Grandes pasos de tiempo ; Transmisión de calor transitoria

Keywords

Finite elements ; Explicit method ; Large time steps ; Transient heat transfer

1. Introducción

El análisis de la transmisión de calor no estacionaria presenta un indudable interés práctico no solamente debido a la importancia que los procesos de enfriamiento y calentamiento poseen en un gran número de aplicaciones industriales, sino también por su similitud con un sinnúmero de otras ecuaciones de la física-matemática que presentan las mismas dificultades para ser resueltas numéricamente. Entre los problemas relacionados directamente a problemas de transmisión del calor no estacionaria cabe citar, por ejemplo, la problemática asociada con el arranque y parada de máquinas, procesos que pueden conllevar, por un lado, la aparición de importantes cargas térmicas en el sólido con el consiguiente deterioro (álabes de turbinas, por ejemplo), y por otro lado, una mayor emisión de especies contaminantes (generadores de vapor, turbinas de gas, motores alternativos de combustión interna, etc.). Entre las ecuaciones diferenciales que tienen la misma problemática para ser resueltas numéricamente, se pueden citar las ecuaciones de conservación de momento, tanto en sólidos como en fluidos, donde en las mismas aparecen términos en derivadas parciales en el espacio de segundo orden y temporales de primer orden.

Por los motivos mencionados anteriormente, la resolución eficiente y precisa de las ecuaciones temporales de transmisión del calor, juega un papel fundamental, no solo por la cantidad de problemas ingenieriles que ellas representan, sino también como caso inicial o de prueba para un número importante de otros problemas más complejos, en particular las ecuaciones de la mecánica de fluidos, cuya solución correcta depende de la eficiencia mostrada por los diferentes métodos en la resolución de las ecuaciones de transmisión del calor [1] .

La integración temporal numérica de las ecuaciones se puede hacer en forma explícita o implícita. En la primera, su principal ventaja es la posibilidad de resolver cada paso de tiempo en forma local, sin necesidad de resolver un sistema de ecuaciones donde se involucran la totalidad de los grados de libertad. Su desventaja es la estabilidad condicional de sus resultados lo que en muchos casos conlleva a pasos de tiempo demasiados pequeños para ser eficientes. Por el contrario, los métodos implícitos permiten, bajo ciertas condiciones, pasos de tiempo siempre estables, cuyo tamaño debe ser regulado solo por criterios de precisión y no de estabilidad. A cambio, se debe resolver un sistema de ecuaciones que involucra a la totalidad de sus grados de libertad, lo cual hace que el método sea costoso y de difícil paralelización.

El creciente uso de ordenadores basados en multiprocesadores, ya sea usando programación en paralelo o GPU, ha incentivado últimamente el uso de los métodos explícitos, muchos mas aptos a ser paralelizados debido a su condición de tratamiento local de la información [2] . Podríamos decir que un método explícito, aunque sea mas caro que un método implícito usando un ordenador serial, tiene muchas mas posibilidades de ser más eficiente en un ordenador paralelo, incluso es muy posible que siempre lo llegue a superar, dependiendo de la cantidad de procesadores en paralelo que se use [3] .

Este trabajo trata de cómo desarrollar métodos explícitos que sean estables para pasos de tiempo mucho mayores que los normalmente utilizados en este tipo de integración temporal. Los pasos de tiempo que se logran son, para una precisión similar, del mismo orden que los utilizados en los métodos implícitos. De esta forma se mantienen todas las ventajas referentes a la facilidad del tratamiento paralelo propio de los explícitos con las ventajas de los implícitos referentes a su estabilidad para grandes pasos de tiempo.

En este trabajo nos concentramos en el problema de conducción pura, aunque en la Sección 4, se plantea como generalizarlo para problemas de conducción con convección.

2. Ecuaciones a resolver y su integración temporal

Las ecuaciones de transmisión del calor en un medio que se mueve con una velocidad conocida V(x''','''t''') se las puede escribir como:

${\begin{array}{c}C{\frac {DT}{Dt}}=\nabla (\kappa \nabla T)+Q{\mbox{ }}{\mbox{en}}\Omega {\mbox{ }}\\C{\frac {Dx}{Dt}}=V{\mbox{ }}{\mbox{en}}\Omega \end{array}}{\begin{array}{c}T={\overline {T}}{\mbox{ }}{\mbox{en}}{\mbox{   }}{\Gamma }_{D}\\\nabla T{\boldsymbol {n}}={\overline {q}}_{n}{\mbox{ }}{\mbox{en}}{\mbox{   }}{\Gamma }_{N}\end{array}}$

( 2.1)

donde

${\boldsymbol {\mbox{x}}}={\boldsymbol {\mbox{x}}}(t)={\boldsymbol {\mbox{x}}}^{t};{\mbox{ }}T=$$T({\boldsymbol {\mbox{x}}}^{t},t)=T^{t}({\boldsymbol {\mbox{x}}}^{t});{\mbox{ }}{\boldsymbol {\mbox{V}}}=$${\boldsymbol {\mbox{V}}}({\boldsymbol {\mbox{x}}}^{t},t)={\boldsymbol {\mbox{V}}}^{t}({\boldsymbol {\mbox{x}}}^{t});{\mbox{ }}Q=$$Q({\boldsymbol {\mbox{x}}}^{t},t)=Q^{t}({\boldsymbol {\mbox{x}}}^{t})$

$\kappa =\left[{\frac {Kg{\mbox{ }}m}{s^{4}{\mbox{ }}K}}\right];{\mbox{   }}C=$$\left[{\frac {Kg}{s^{2}{\mbox{ }}m^{0}{\mbox{ }}K}}\right];{\mbox{   }}{\frac {\kappa }{C}}=$$\left[{\frac {m^{2}}{s}}\right]$

representan, la posición de una partícula en el espacio en función del tiempo, la temperatura, el campo conocido de velocidad del material, una fuente externa de calor, el coeficiente de conductividad térmica y el coeficiente de capacidad calorífica. El término ${\textstyle D\varphi /Dt}$ representa la derivada material de cualquier función ϕ   , ${\textstyle {\overline {T}}}$ representa la temperatura impuesta en los contornos Γ_D donde hay condiciones de Dirichlet y ${\textstyle {\overline {q}}_{n}}$ el flujo térmico en los contornos Γ_N donde hay condiciones de Neumman

La ecuación (2.1) está escrita en forma Lagrangiana. Su forma Euleriana se obtiene reemplazando la derivada material por la derivada espacial más el término convectivo (Dφ /Dt ) = ∂φ /∂t  + V ∇ φ .

La integración temporal de primer orden estándar de la ecuación anterior entre 2 instantes de tiempo tⁿ y t^{n +1}  = tⁿ  + Δt es:

${\begin{array}{c}T^{n+1}({\boldsymbol {\mbox{x}}}^{n+1})=T^{n}({\boldsymbol {\mbox{x}}}^{n})+{\left\{\nabla \left[{\frac {\kappa }{C}}\nabla T({\boldsymbol {\mbox{x}}})\right]+{\frac {Q({\boldsymbol {\mbox{x}}})}{C}}\right\}}^{n+\theta }{\mbox{ }}\Delta t\\{\boldsymbol {\mbox{x}}}^{n+1}={\boldsymbol {\mbox{x}}}^{n}+{\left[{\boldsymbol {\mbox{V}}}({\boldsymbol {\mbox{x}}})\right]}^{n+\theta }{\mbox{ }}\Delta t\end{array}}$

( 2.2)

donde el coeficiente θ puede variar entre 0 y 1. Para θ=0, el método es explícito y para θ≠0 el método es implícito.

Para θ<1/2 el método es condicionalmente estable. En particular para el caso explícito (θ=0), los límites de estabilidad vienen dados por 2 coeficientes adimensionales, el número de Courant ${\textstyle \left(Co={\frac {\left|{\boldsymbol {\mbox{V}}}\right|\Delta t}{2\delta }}\right)}$ y el numero de Fourier ${\textstyle \left(Fo={\frac {{\frac {\kappa }{C}}\Delta t}{{\delta }^{2}}}\right)}$ .

El parámetro δ representa el tamaño mínimo en el que se puede aproximar la variación de la función incógnita (en nuestro caso la temperatura) en la discretización. Por ejemplo, en un problema unidimensional, δ representa la menor distancia entre 2 nodos, pero en un problema tridimensional discretizado con tetraedros, δ es la menor altura de todos los tetraedros.

El límite de estabilidad para el caso explícito viene dado por ${\textstyle (Co+Fo)\leq {\frac {1}{2}}}$ .

En un trabajo anterior [4] los autores explicaron la forma de resolver problemas con grandes pasos de tiempo utilizando una integración temporal explícita para los casos dominados por convección (grandes números de Courant y pequeños números de Fourier). Por este motivo en este trabajo solo se tratará el caso dominado por difusión (Co =0), enviando al lector interesado a la referencia [4] para los casos dominados por convección. En la mencionada referencia, la integración explícita para obtener métodos estables para cualquier paso de tiempo en el caso dominado por convección, se la denominó integración X-IVAS. Será repasada en este trabajo para facilidad del lector. En el capítulo cuarto introduciremos nuestra propuesta para combinar ambas metodologías para tratar problemas de conducción-convección, integrados temporalmente en forma explícita, estables y con grandes pasos de tiempo.

3. Caso con conducción sin convección

Definiremos la función aceleración o tasa de crecimiento térmico g (x , t ) = g^t (x ) como:

$g^{t}({\boldsymbol {\mbox{x}}})=\nabla \left[{\frac {\kappa }{C}}\nabla T^{t}({\boldsymbol {\mbox{x}}})\right]+$${\frac {Q^{t}({\boldsymbol {\mbox{x}}})}{C}}$

( 3.1)

Para aumentar la estabilidad de la solución explícita, el cálculo de la función g^t (x ) en cualquier punto x_i del espacio lo calcularemos aplicando el siguiente kernel:

$g^{t}({\boldsymbol {\mbox{x}}}_{i})\cong {\overset {\frown }{g}}^{t}({\boldsymbol {\mbox{x}}}_{i})=$${\frac {1}{\gamma }}\int _{\Omega }\varphi ({\boldsymbol {\mbox{x}}}-$${\boldsymbol {\mbox{x}}}_{i})g^{t}(x)d\Omega$

( 3.2)

donde φ(x) es una función de peso arbitraria pero que cumple al menos que:

${\begin{array}{c}\varphi ({\boldsymbol {\mbox{x}}})=0{\mbox{   }}\forall {\boldsymbol {\mbox{x}}}/{\mbox{   }}\left|{\boldsymbol {\mbox{x}}}\right|>R_{ball}\\\int _{\Omega }\varphi ({\boldsymbol {\mbox{x}}})d\Omega =\gamma \end{array}}$

( 3.3)

En la ecuación (3.3), R_ball representa el radio del volumen Ω que tiene relación con la distancia R^* hasta donde llega la información suponiendo que la misma se desplaza en forma radial con una velocidad ficticia V^* :

$V^{\ast }={\sqrt {\frac {2\kappa }{C\Delta t}}}$

( 3.4)

Para lograr integraciones explicitas estables, el radio R^* debe ser:

$R^{\ast }\geq V^{\ast }\Delta t$

( 3.5)

Por tanto

$\Delta t\leq {\frac {{\left(R^{\ast }\right)}^{2}C}{2\kappa }}{\mbox{ }}{\mbox{o}}{\mbox{ }}{\mbox{también}}{\mbox{ }}Fo^{\ast }\leq {\frac {1}{2}}$

( 3.6)

donde hemos definido un número de Fourier ficticio como ${\textstyle Fo^{\ast }={\frac {\kappa }{C}}{\frac {\Delta t}{{\left(R^{\ast }\right)}^{2}}}}$

En la Sección 5 de este artículo desarrollaremos la teoría sobre la condición de estabilidad (3.5) o (3.6) y la relación existente entre R_ball y R^* .

Entonces, la integración temporal explicita que en la forma estándar es:

$T^{n+1}({\boldsymbol {\mbox{x}}}_{i})=T^{n}({\boldsymbol {\mbox{x}}}_{i})+$$\Delta t\nabla \left\{\left[{\frac {\kappa }{C}}\nabla T^{n}({\boldsymbol {\mbox{x}}})\right]+\right.$$\left.{\frac {Q^{n}({\boldsymbol {\mbox{x}}})}{C}}\right\}$

( 3.7)

queda ahora usando la aproximación (3.2):

$T^{n+1}({\boldsymbol {\mbox{x}}}_{i})=T^{n}({\boldsymbol {\mbox{x}}}_{i})+$$\Delta t{\frac {1}{\gamma }}\int _{\Omega }\varphi ({\boldsymbol {\mbox{x}}}-$${\boldsymbol {\mbox{x}}}_{i})\left\{\nabla \left[{\frac {\kappa }{C}}\nabla T^{n}({\boldsymbol {\mbox{x}}})\right]+\right.$$\left.{\frac {Q^{n}({\boldsymbol {\mbox{x}}})}{C}}\right\}d\Omega$

( 3.8)

Existen en la literatura otras propuestas para ampliar el límite de estabilidad de la integración explícita [5] , [6] , [7]  and [8] . Por ejemplo en la referencia [5] , se propone hacer pasar un polinomio de grado M en cada punto a integrar y se toman tantos puntos en su vecindad como los necesarios para definir biunívocamente un polinomio de grado M . El método solo se lo utiliza para casos con distribución de puntos regulares dado el costo computacional que significaría utilizarlo para distribuciones arbitrarias de nodos. En una forma mucho más elemental la referencia [6] propone simplemente usar para la definición de la derivada segunda, puntos más alejados que los vecinos. Lógicamente solo se presentan ejemplos unidimensionales debido a la imposibilidad de generalizar esta metodología a espacios mayores. La propuesta realizada en este trabajo es aplicable a problemas tanto en 1, 2 como en 3-dimensiones y con distribuciones arbitrarias de nodos, como se podrá apreciar en los ejemplos realizados.

3.1. Dizcretización espacial

Utilizaremos para la discretización espacial de la temperatura el método de los elementos finitos (FEM):

$T^{t}({\boldsymbol {\mbox{x}}})={\boldsymbol {\mbox{N}}}^{T}({\boldsymbol {\mbox{x}}}){\boldsymbol {\mbox{T}}}^{t}$

( 3.9)

donde el vector N(x) representa las clásicas funciones de forma del FEM y el vector T los valores locales de temperatura en los nodos de la malla.

De igual forma, para el cálculo de la integrales dentro de (3.8), aproximaremos la función que recién denomináramos como aceleración ${\textstyle \nabla \left[{\frac {\kappa }{C}}\nabla T^{n}({\boldsymbol {\mbox{x}}})\right]+}$${\frac {Q^{n}({\boldsymbol {\mbox{x}}})}{C}}$ por una función continua con las mismas funciones de forma que las usadas para la temperatura:

$g^{n}({\boldsymbol {\mbox{x}}})={\boldsymbol {\mbox{N}}}^{T}({\boldsymbol {\mbox{x}}}){\mbox{ }}{\boldsymbol {\mbox{g}}}^{n}$

( 3.10)

El valor de gⁿ (x ) lo obtenemos planteando la ecuación en residuos ponderados:

$\int _{\Omega }{\boldsymbol {\mbox{N}}}\left\{\nabla \left[{\frac {\kappa }{C}}\nabla T^{n}({\boldsymbol {\mbox{x}}})\right]+\right.$$\left.{\frac {Q^{n}({\boldsymbol {\mbox{x}}})}{C}}-{\boldsymbol {\mbox{N}}}^{T}({\boldsymbol {\mbox{x}}}){\mbox{ }}{\boldsymbol {\mbox{g}}}^{n}\right\}d\Omega +$$\int _{{\Gamma }_{N}}{\boldsymbol {\mbox{N}}}\left({\frac {{\overline {q}}_{N}}{C}}-\right.$$\left.{\frac {\kappa }{C}}\nabla T^{n}({\boldsymbol {\mbox{x}}}){\boldsymbol {\mbox{n}}}\right)d\Gamma =$$0$

( 3.11)

que luego de integrar por partes el primer termino queda:

$-\int _{\Omega }\nabla {\boldsymbol {\mbox{N}}}\left\{{\frac {\kappa }{C}}\nabla T^{n}({\boldsymbol {\mbox{x}}})\right\}d\Omega +$$\int _{\Omega }{\boldsymbol {\mbox{N}}}{\frac {Q^{n}({\boldsymbol {\mbox{x}}})}{C}}d\Omega +$${\underset {{\Gamma }_{N}}{\oint }}{\boldsymbol {\mbox{N}}}{\frac {{\overline {q}}_{N}}{C}}{\mbox{ }}d\Gamma -$$\int _{\Omega }{\boldsymbol {\mbox{N}}}{\boldsymbol {\mbox{N}}}^{T}d\Omega {\mbox{ }}{\boldsymbol {g}}^{n}=$$0$

( 3.12)

o en forma matricial:

$-{\boldsymbol {\mbox{K}}}{\boldsymbol {\mbox{T}}}^{n}+{\boldsymbol {\mbox{F}}}^{n}=$${\boldsymbol {\mbox{M}}}{\boldsymbol {\mbox{g}}}^{n}$

( 3.13)

o sea:

${\boldsymbol {\mbox{g}}}^{n}={\boldsymbol {\mbox{M}}}^{-1}({\boldsymbol {\mbox{F}}}^{n}-$${\boldsymbol {\mbox{K}}}{\boldsymbol {\mbox{T}}}^{n})$

( 3.14)

donde:

${\begin{array}{c}{\boldsymbol {\mbox{F}}}^{n}=\int _{\Omega }{\boldsymbol {\mbox{N}}}{\frac {Q^{n}({\boldsymbol {\mbox{x}}})}{C}}d\Omega +{\underset {{\Gamma }_{N}}{\oint }}{\boldsymbol {\mbox{N}}}{\frac {{\overline {q}}_{N}}{C}}{\mbox{ }}d\Gamma ;\\{\boldsymbol {\mbox{K}}}=\int _{\Omega }\nabla {\boldsymbol {\mbox{N}}}{\frac {\kappa }{C}}\nabla T^{n}({\boldsymbol {\mbox{x}}})d\Omega {\mbox{   }};{\mbox{   }}{\boldsymbol {\mbox{M}}}=\int _{\Omega }{\boldsymbol {\mbox{N}}}{\boldsymbol {\mbox{N}}}^{T}d\Omega \end{array}}$

( 3.15)

Por lo tanto la expresión final para evaluar en forma explícita el valor de la temperatura en el tiempo n+1 queda:

${\begin{array}{c}{\boldsymbol {\mbox{g}}}^{n}={\boldsymbol {\mbox{M}}}^{-1}({\boldsymbol {\mbox{F}}}^{n}-{\boldsymbol {\mbox{K}}}{\boldsymbol {\mbox{T}}}^{n})\\{\boldsymbol {\mbox{T}}}^{n+1}({\boldsymbol {\mbox{x}}}_{i})={\boldsymbol {\mbox{T}}}^{n}({\boldsymbol {\mbox{x}}}_{i})+\Delta t{\frac {1}{\gamma }}\int _{\Omega }\varphi ({\boldsymbol {\mbox{x}}}-{\boldsymbol {\mbox{x}}}_{i}){\boldsymbol {\mbox{N}}}^{T}({\boldsymbol {\mbox{x}}})d\Omega {\mbox{ }}{\boldsymbol {\mbox{g}}}^{n}\\\gamma =\int _{\Omega }\varphi ({\boldsymbol {\mbox{x}}})d\Omega \end{array}}$

( 3.16)

3.2. Evaluación de las integrales espaciales

Ambas integrales espaciales en (3.16) se las pueden calcular de diversas maneras, ya sea en forma analítica o numérica. Evidentemente, el costo computacional del método propuesto dependerá de la eficiencia en calcular dichas integrales. En este trabajo proponemos hacerlo en una forma numérica muy sencilla que tiene la ventaja de ser muy eficiente desde el punto de vista computacional. Para esto subdividimos el dominio Ω en varios subdominios δ Ω de forma que:

${\begin{array}{c}\gamma =\int _{\Omega }\varphi ({\boldsymbol {\mbox{x}}}-{\boldsymbol {\mbox{x}}}_{i})d\Omega \approx \sum _{l=1}^{m}\left[\varphi ({\boldsymbol {\mbox{x}}}_{l}-{\boldsymbol {\mbox{x}}}_{i}){\mbox{ }}\delta {\Omega }_{l}\right]\\\int _{\Omega }\varphi ({\boldsymbol {\mbox{x}}}-{\boldsymbol {\mbox{x}}}_{i}){\mbox{ }}g^{n}({\boldsymbol {\mbox{x}}}){\mbox{ }}d\Omega \approx \sum _{l=1}^{m}\left[\varphi ({\boldsymbol {\mbox{x}}}_{l}-{\boldsymbol {\mbox{x}}}_{i}\right)g^{n}({\boldsymbol {\mbox{x}}}_{l})]{\mbox{ }}\delta {\Omega }_{l}\end{array}}$

( 3.17)

donde x_l son los diferentes puntos de integración.

Si consideramos como puntos de integración los mismos nodos de la malla, para cada nodo podemos definir una integral:

${\begin{array}{c}{\gamma }_{i}=\int _{\Omega }\varphi ({\boldsymbol {\mbox{x}}}-{\boldsymbol {\mbox{x}}}_{i})d\Omega \approx \sum _{j=1}^{m}\left[\varphi ({\boldsymbol {\mbox{x}}}_{j}-{\boldsymbol {\mbox{x}}}_{i}){\mbox{ }}\delta {\Omega }_{j}\right]\\\int _{\Omega }{\varphi }_{i}({\boldsymbol {\mbox{x}}}-{\boldsymbol {\mbox{x}}}_{i})g^{n}({\boldsymbol {\mbox{x}}})d\Omega \approx \sum _{j=1}^{m}[\varphi ({\boldsymbol {\mbox{x}}}_{j}-{\boldsymbol {\mbox{x}}}_{i}){\mbox{ }}g^{n}({\boldsymbol {\mbox{x}}}_{j})]\delta {\Omega }_{j}\end{array}}$

( 3.18)

Ahora los x_j de la ecuación 3.18 son los propios nodos de la malla que están dentro del espacio donde ϕ_i ≠0.

Sea M^L la matriz agrupada (del inglés lumped ) de la matriz M . Teniendo en cuenta que cada uno de sus términos diagonales es el volumen asociado a cada nodo ${\textstyle \left(\delta {\Omega }_{i}=M_{ii}^{L}\right)}$ , entonces:

${\gamma }_{i}=\int _{\Omega }{\varphi }_{i}({\boldsymbol {\mbox{x}}})d\Omega \approx traza\left[{\Phi }_{i}^{T}{\boldsymbol {\mbox{M}}}^{L}\right]{\mbox{   }}{\mbox{y}}{\mbox{   }}\int _{\Omega }{\varphi }_{i}({\boldsymbol {\mbox{x}}}){\mbox{ }}g^{n}({\boldsymbol {\mbox{x}}}){\mbox{ }}d\Omega \approx {\Phi }_{i}^{T}{\boldsymbol {\mbox{M}}}^{L}{\boldsymbol {\mbox{g}}}^{n}$

( 3.19)

donde ${\textstyle {\Phi }_{i}^{T}=[{\phi }_{i}({\boldsymbol {\mbox{x}}}_{1}),{\phi }_{i}({\boldsymbol {\mbox{x}}}_{2}),{\phi }_{i}({\boldsymbol {\mbox{x}}}_{3}),...{\phi }_{i}({\boldsymbol {\mbox{x}}}_{i}),...{\phi }_{i}({\boldsymbol {\mbox{x}}}_{m})]}$ .

Entonces (3.16) se puede escribir como:

${\begin{array}{c}{\boldsymbol {\mbox{g}}}^{n}={\boldsymbol {\mbox{M}}}^{-1}({\boldsymbol {\mbox{F}}}^{n}-{\boldsymbol {\mbox{K}}}{\boldsymbol {\mbox{T}}}^{n})\\{\boldsymbol {\mbox{T}}}^{n+1}={\boldsymbol {\mbox{T}}}^{n}+\Delta t{\mbox{ }}{\Phi }^{T}{\boldsymbol {\mbox{M}}}^{L}{\boldsymbol {\mbox{g}}}^{n}\end{array}}$

( 3.20)

en donde hemos agrupado dentro de la matriz Φ^T los diferentes vectores ${\textstyle {\Phi }_{i}^{T}}$ divididos por su correspondiente γ_i

${\Phi }^{T}=\left[{\begin{array}{c}{\Phi }_{1}^{T}/{\gamma }_{1}\\{\Phi }_{2}^{T}/{\gamma }_{2}\\{\Phi }_{3}^{T}/{\gamma }_{3}\\\ldots \end{array}}\right]$

( 3.21)

Para el caso en que se utilizan matrices de masa agrupada (aproximación estándar en los métodos explícitos), el algoritmo a resolver queda simplemente:

${\boldsymbol {\mbox{T}}}^{n+1}={\boldsymbol {\mbox{T}}}^{n}+$$\Delta t{\mbox{   }}{\Phi }^{T}({\boldsymbol {\mbox{F}}}^{n}-$${\boldsymbol {\mbox{K}}}{\boldsymbol {\mbox{T}}}^{n})$

( 3.22)

Como resumen el método explicito estándar es:

${\boldsymbol {\mbox{T}}}^{n+1}={\boldsymbol {\mbox{T}}}^{n}+$$\Delta t{\mbox{   }}{\boldsymbol {\mbox{M}}}^{-1}({\boldsymbol {\mbox{F}}}^{n}-$${\boldsymbol {\mbox{K}}}{\boldsymbol {\mbox{T}}}^{n})$

( 3.23)

y el nuevo método estabilizado para grandes pasos de tiempo queda:

${\boldsymbol {\mbox{T}}}^{n+1}={\boldsymbol {\mbox{T}}}^{n}+$$\Delta t{\mbox{   }}[{\Phi }^{T}{\boldsymbol {\mbox{M}}}^{L}]{\mbox{   }}{\boldsymbol {\mbox{M}}}^{-1}({\boldsymbol {\mbox{F}}}^{n}-$${\boldsymbol {\mbox{K}}}{\boldsymbol {\mbox{T}}}^{n})$

( 3.24)

Es de hacer notar que la matriz Φ depende exclusivamente de la malla y del radio r elegido para las funciones ϕ_i . Mientras se cumpla la Ecuación 3.6, la resolución será estable.

3.3. Condiciones de contorno

Supondremos 2 tipos de condiciones de contorno:

• Γ_T , contorno donde se conoce la temperatura ${\textstyle {\overline {T}}({\boldsymbol {\mbox{x}}}_{{\Gamma }_{T}})}$ y donde supondremos que la misma no varía en el tiempo ${\textstyle \left({\frac {DT}{Dt}}_{{\Gamma }_{T}}=\right.}$$\left.0\right)$ .
• Γ_q , contorno adiabático ${\textstyle \nabla T({\boldsymbol {\mbox{x}}}_{{\Gamma }_{q}}){\mbox{ }}{\boldsymbol {\mbox{n}}}({\boldsymbol {\mbox{x}}}_{{\Gamma }_{q}})=}$$0$

Esta última es la condición natural, por lo tanto en las fronteras donde no se impone nada, serán consideradas fronteras adiabáticas.

Como es estándar en FEM en el contorno Γ_T las funciones de forma ${\textstyle {\boldsymbol {\mbox{N}}}({\boldsymbol {\mbox{x}}}_{{\Gamma }_{T}})}$ de aproximación de la temperatura son nulas, pero como además aquí estamos representando también la función aceleración si esta condición de contorno no cambia en el tiempo entonces hay que imponer también las funciones de forma ${\textstyle {\boldsymbol {\mbox{N}}}({\boldsymbol {\mbox{x}}}_{{\Gamma }_{T}})}$ de aproximación de la ${\textstyle {\frac {DT}{Dt}}_{{\Gamma }_{T}}}$ a un valor nulo. En forma equivalente se puede hacer cero los valores de gⁿ calculados en los nodos pertenecientes a esas fronteras.

4. Caso con conducción y convección

Si bien en este artículo no se mostrarán ejemplos del caso combinado de problemas de conducción-convección resueltos en forma explicita, se presenta a continuación una propuesta de cómo resolverlo. Utilizaremos una modificación del método X-IVAS propuesto en [4] .

El método X-IVAS consiste en integrar temporalmente dentro de cada paso de tiempo la posición de las partículas de material siguiendo las líneas de corriente:

${\begin{array}{c}T^{n+1}({\boldsymbol {\mbox{x}}}^{n+1})=T^{n}({\boldsymbol {\mbox{x}}}^{n})+\int _{n}^{n+1}\left\{\nabla \left[{\frac {\kappa }{C}}\nabla T^{n}({\boldsymbol {\mbox{x}}}^{\tau })\right]+{\frac {Q^{n}({\boldsymbol {\mbox{x}}}^{\tau })}{C}}\right\}d\tau \\{\boldsymbol {\mbox{x}}}^{n+1}={\boldsymbol {\mbox{x}}}^{n}+\int _{n}^{n+1}{\boldsymbol {\mbox{V}}}^{n}({\boldsymbol {\mbox{x}}}^{\tau })d\tau \end{array}}$

( 4.1)

Este método tiene la ventaja que es explícito y estable para ${\textstyle Fo\leq 1/2}$ independientemente del número de Courant. A continuación propondremos una modificación del método X-IVAS, basada en la forma de integrar desarrollada en la sección 3, para salvar también el inconveniente de la necesidad de mantener el paso de tiempo suficientemente pequeño y poder lograr un método de integración explícito con ${\textstyle Co>1{\mbox{   }}{\mbox{y}}{\mbox{   }}Fo>1/2}$ .

Si llamamos ${\textstyle {\overset {\frown }{\boldsymbol {\mbox{g}}}}^{n}}$ al vector de aceleraciones modificado según (3.20)

${\overset {\mbox{ˆ}}{\boldsymbol {\mbox{g}}}}^{n}={\Phi }^{T}{\boldsymbol {\mbox{M}}}^{L}{\boldsymbol {\mbox{g}}}^{n}$

( 4.2)

entonces el método descrito anteriormente para conducción pura se escribirá:

${\begin{array}{c}{\overset {\mbox{ˆ}}{\boldsymbol {\mbox{g}}}}^{n}={\Phi }^{T}({\boldsymbol {\mbox{F}}}^{n}-{\boldsymbol {\mbox{K}}}{\boldsymbol {\mbox{T}}}^{n})\\{\boldsymbol {\mbox{T}}}^{n+1}={\boldsymbol {\mbox{T}}}^{n}+\Delta t{\overset {\mbox{ˆ}}{\boldsymbol {\mbox{g}}}}^{n}\end{array}}$

( 4.3)

Cuando también exista convección, la aplicación de X-IVAS a este esquema queda:

${\begin{array}{c}{\overset {\mbox{ˆ}}{\boldsymbol {\mbox{g}}}}^{n}={\Phi }^{T}({\boldsymbol {\mbox{F}}}^{n}-{\boldsymbol {\mbox{K}}}{\boldsymbol {\mbox{T}}}^{n})\\T^{n+1}({\boldsymbol {\mbox{x}}}^{n+1})=T^{n}({\boldsymbol {\mbox{x}}}^{n})+\int _{n}^{n+1}{\overset {\mbox{ˆ}}{g}}^{n}({\boldsymbol {\mbox{x}}}^{\tau })dt=T^{n}({\boldsymbol {\mbox{x}}}^{n})+\int _{n}^{n+1}{\boldsymbol {\mbox{N}}}^{T}({\boldsymbol {\mbox{x}}}^{\tau })d\tau {\overset {\mbox{ˆ}}{\boldsymbol {\mbox{g}}}}^{n}\\{\boldsymbol {\mbox{x}}}^{t}={\boldsymbol {\mbox{x}}}^{n}+\int _{n}^{n+t}{\boldsymbol {\mbox{V}}}^{n}({\boldsymbol {\mbox{x}}}^{\tau })d\tau \end{array}}$

( 4.4)

Es de hacer notar que la matriz Φ se la puede calcular una sola vez antes de empezar la integración. Lo mismo la matriz K , si los coeficientes κ /C no cambian con el tiempo.

Por tanto la evaluación del vector ${\textstyle {\overset {\frown }{\boldsymbol {\mbox{g}}}}^{n}}$ es explicita y se puede hacer antes de comenzar cada paso de tiempo. Por último la evaluación de T^{n +1} (x^{n +1} ) de cualquier partícula por separado, se la puede calcular integrando las funciones de forma ${\textstyle \int _{n}^{n+1}{\boldsymbol {\mbox{N}}}^{T}({\boldsymbol {\mbox{x}}}^{\tau })d\tau }$ y las velocidades ${\textstyle \int _{n}^{n+t}{\boldsymbol {\mbox{V}}}^{n}({\boldsymbol {\mbox{x}}}^{\tau })d\tau }$ para lo cual se pueden utilizar métodos numéricos (dividiendo el paso de tiempo Δt   en varios sub pasos ${\textstyle \left(\delta t=\Delta t/nsubpasos\right)}$ o realizando una integración analítica [9]  and [10] .

Por último es conveniente remarcar que por tratarse de una formulación Lagrangiana y al haber convección, la posición x^{n +1} de cualquier partícula no es la misma que la posición en xⁿ . Por tanto, el método se puede utilizar en una formulación con malla móvil, moviendo los nodos de la malla a cada paso de tiempo de xⁿ a x^{n +1} , o en una formulación con malla fija, previa proyección de las temperaturas sobre los nodos de la malla Ref [9]  and [10] .

5. Consistencia y condiciones de estabilidad

5.1. Análisis de la consistencia

Una de las primeras ideas fue verificar que el método desarrollado fuera de primer orden en el tiempo y segundo orden en el espacio como lo es el esquema centrado en el espacio y explícito en el tiempo denominado comúnmente Forward Euler, en adelante, el método estándar. Recordemos que el método explícito estándar tiene como ecuación discreta en el caso unidimensional la siguiente expresión, que para mallas regulares coincide con el método de las diferencias finitas:

${\begin{array}{c}T^{n+1}(x_{i})=T^{n}(x_{i})+\Delta tg^{n}(x_{i})=T^{n}(x_{i})+{\frac {\kappa }{C}}\Delta t{\frac {T^{n}(x_{i-1})-2T^{n}(x_{i})+T^{n}(x_{i+1})}{\Delta x^{2}}}\\T^{n+1}(x_{i})=T^{n}(x_{i})+Fo(T^{n}(x_{i-1})-2T^{n}(x_{i})+T^{n}(x_{i+1}))\end{array}}$

( 5.1)

5.2. Si desarrollamos por Taylor el esténcil anterior usando

$T^{n+1}(x_{i})=T^{n}(x_{i})+{\frac {\partial T^{n}}{\partial t}}(x_{i})\Delta t+$$O{\left(\Delta t\right)}^{2}$

( 5.2)

$T^{n}(x_{i\pm 1})={T_{i}}^{n}\pm {\frac {\partial T^{n}}{\partial x}}(x_{i})\Delta x+$${\frac {{\partial }^{2}T^{n}}{\partial x^{2}}}(x_{i}){\frac {{\left(\Delta x\right)}^{2}}{2}}\pm {\frac {{\partial }^{3}T^{n}}{\partial x^{3}}}(x_{i}){\frac {{\left(\Delta x\right)}^{3}}{6}}+$$O{\left(\Delta x\right)}^{4}$

( 5.3)

y reemplazando en la anterior vemos que

${\begin{array}{c}T^{n}(x_{i})+{\frac {\partial T^{n}}{\partial t}}(x_{i})\Delta t+O{\left(\Delta t\right)}^{2}=T^{n}(x_{i})+Fo\left({\frac {{\partial }^{2}T^{n}}{\partial x^{2}}}(x_{i}){\left(\Delta x\right)}^{2}+O{\left(\Delta x\right)}^{4}\right)\\{\frac {\partial T^{n}}{\partial t}}(x_{i})\Delta t+O{\left(\Delta t\right)}^{2}=Fo\left({\frac {{\partial }^{2}T^{n}}{\partial x^{2}}}(x_{i}){\left(\Delta x\right)}^{2}+O{\left(\Delta x\right)}^{4}\right)\\{\frac {\partial T^{n}}{\partial t}}(x_{i})+O(\Delta t)={\frac {\kappa }{C}}{\frac {{\partial }^{2}T^{n}}{\partial x^{2}}}(x_{i})+O{\left(\Delta x\right)}^{2}\end{array}}$

( 5.4)

Utilizando la ecuación del continuo se establece que el esquema es consistente (la diferencia entre la ecuación a resolver y la versión discreta de la misma va a cero cuando ambos pasos discretos, el del tiempo y el de la malla tienden a cero) y nos queda que el error del método depende linealmente con el paso de tiempo y cuadráticamente con el espaciamiento de la malla.

Del mismo modo se puede analizar la consistencia del esquema explícito modificado. En ese caso hemos usado un polinomio φ(r ), o función de peso, definido como:

$\varphi (r)=1-2r+r^{2}$ $r(x)={\frac {\left\|x-x_{i}\right\|}{R_{ball}}}$

( 5.5)

con x_i la coordenada del nodo i en cuestión y R_ball un radio de influencia que marca aproximadamente el semi-ancho del esténcil. Siendo que el esquema modificado se escribe como:

$T^{n+1}(x_{i})=T^{n}(x_{i})+{\frac {\Delta t}{{\sum }_{k=-N}^{k=+N}{\varphi }_{i}(x_{i+k})d{\Omega }_{i+k}}}\sum _{j=-N}^{j=+N}{\varphi }_{i}(x_{i+j})g^{n}(x_{i+j})d{\Omega }_{i+j}$

( 5.6)

Como el cómputo es tedioso de hacerlo en forma manual se recurrió a un análisis simbólico del cual, escribiendo la ecuación que queda después de reemplazar cada término por su correspondiente expansión de Taylor surgió lo siguiente:

${\begin{array}{c}T^{n}(x_{i})+{\frac {\partial T^{n}}{\partial t}}(x_{i})\Delta t+O{\left(\Delta t\right)}^{2}=T^{n}(x_{i})+Fo\left({\beta }_{0}T^{n}(x_{i})+{\beta }_{2}{\frac {{\partial }^{2}T^{n}}{\partial x^{2}}}(x_{i}){\left(\Delta x\right)}^{2}+{\beta }_{4}{\frac {{\partial }^{4}T^{n}}{\partial x^{4}}}(x_{i}){\left(\Delta x\right)}^{4}\right)\\{\frac {\partial T^{n}}{\partial t}}(x_{i})\Delta t+O{\left(\Delta t\right)}^{2}=Fo\left({\beta }_{0}T^{n}(x_{i})+{\beta }_{2}{\frac {{\partial }^{2}T^{n}}{\partial x^{2}}}(x_{i}){\left(\Delta x\right)}^{2}+{\beta }_{4}{\frac {{\partial }^{4}T^{n}}{\partial x^{4}}}(x_{i}){\left(\Delta x\right)}^{4}\right)\end{array}}$

( 5.7)

de lo cual para poder demostrar que el esquema es consistente se debe cumplir que β₀  ≅ 0 y β₂  ≅ 1.

Para el caso de N =2 surgió lo siguiente:

${\beta }_{0}={\frac {1}{18014398509481984}}\cong 0;{\mbox{   }}{\beta }_{2}=$${\frac {36028797018963969}{36028797018963968}}\cong 1;{\mbox{   }}{\beta }_{4}=$${\frac {6569149026422195}{144115188075855872}}\cong 0,45583$

Para el caso de N =3 surgió lo siguiente:

${\beta }_{0}={\frac {1}{288230376151711744}}\cong 0;{\mbox{   }}{\beta }_{2}=$${\frac {576460752303423501}{576460752303423488}}\cong 1;{\mbox{   }}{\beta }_{4}=$${\frac {2069732352903506911}{2305843009213693952}}\cong 0,89760$

Finalmente para el caso de N =4 tenemos:

${\beta }_{0}={\frac {1}{576460752303423488}}\cong 0;{\mbox{   }}{\beta }_{2}=$${\frac {288230376151711757}{288230376151711744}}\cong 1;{\mbox{   }}{\beta }_{4}=$${\frac {5174850966523186049}{3458764513820540928}}\cong 1,4962$

Por todo lo mostrado podemos asegurar que el esquema modificado goza de ser consistente.

5.3. Análisis de estabilidad

A continuación planteamos un análisis de estabilidad para el problema 1D. Para ello utilizamos el esténcil que genera el método explícito aquí presentado, a saber:

$T^{n+1}(x_{i})=T^{n}(x_{i})+{\frac {\Delta t}{{\sum }_{k=-N}^{k=+N}{\varphi }_{i}(x_{i+k})}}\sum _{j=-N}^{j=+N}{\varphi }_{i}(x_{i+j})g^{n}(x_{i+j})$

( 5.8)

Por cuestiones de simplicidad hemos asumido que la malla es de tamaño constante motivo por el cual hemos obviado en la expresión anterior tanto en el numerador como en el denominador el tamaño de la misma. Reemplazando la expresión de la función gⁿ para el caso de un esténcil 1D de conducción pura con paso de malla constante por:

$g^{n}(x_{i+j})={\frac {\kappa /C}{\Delta x^{2}}}[T^{n}(x_{i-1})-$$2T^{n}(x_{i})+T^{n}(x_{i+1})]$

( 5.9)

Usando la definición del número de Fourier (5.8) queda:

$T^{n+1}(x_{i})=T^{n}(x_{i})+{\frac {Fo}{{\sum }_{k=-N}^{k=+N}{\varphi }_{i}(x_{i+k})}}\sum _{j=-N}^{j=+N}{\varphi }_{i}(x_{i+j})[T^{n}(x_{i+j-1})-$$2T^{n}(x_{i+j})+T^{n}(x_{i+j+1})]$

( 5.10)

que podemos reescribir como:

$T^{n+1}(x_{i})=T^{n}(x_{i})+{\frac {Fo}{\gamma }}\sum _{j=-N}^{j=+N}{\varphi }_{i}(x_{i+j})\left(T^{n}(x_{i+j-1})-\right.$$\left.2T^{n}(x_{i+j})+T^{n}(x_{i+j+1})\right)con\gamma =$${\sum }_{k=-N}^{k=+N}{\varphi }_{i}(x_{i+k})$

Si tomamos por ejemplo un radio de integración de 3Δx obtenemos la función de ponderación φ graficada en la figura 1 . Se puede apreciar que se anula en los extremos del esténcil ${\textstyle \left({\phi }_{r=1}=0\right)}$ por lo que las gⁿ (x_{i ±3} ) no son necesarias.

Figura 1. Función φ tomando Rball =3h .

Utilizando el método de integración propuesto la temperatura en un nodo i para el paso de tiempo n+1 toma la siguiente forma:

${\begin{array}{c}T^{n+1}(x_{i})=T^{n}(x_{i})+Fo/\gamma {\mbox{ }}[T^{n}(x_{i-1})-2T^{n}(x_{i})+T^{n}(x_{i+1})]+\\F_{o}/\gamma {\phi }_{1}{\mbox{ }}[T^{n}(x_{i-2}-2T^{n}(x_{i-1}+2T^{n}(x_{i})-2T(x_{i+1})+T^{n}(x_{i+2})]+\\F_{o}/\gamma {\phi }_{2}[T^{n}(x_{i-3})-2T^{n}(x_{i-2})+T^{n}(x_{i-1}){\mbox{ }}+T^{n}(x_{i+1})-2T^{n}(x_{i+2})+T^{n}(x_{i+3})]\end{array}}$

Considerando ahora la evolución de un único modo de Fourier para la onda k , podemos expresar la temperatura como:

$T^{t}(x)={\overset {\mbox{ˆ}}{T}}(t)e^{ikx}{\mbox{ }}entonces{\mbox{ }}T^{t}(x_{i+1})=$${\overset {\mbox{ˆ}}{T}}(t)e^{ik(x_{i}+h)}$

( 5.11)

Utilizando la definición cos(kh ) = 1/2(e^ikh  + e^−ikh ) en la ecuación 5.11 obtenemos:

${\overset {\mbox{ˆ}}{T}}^{n+1}={\overset {\mbox{ˆ}}{T}}^{n}A$

( 5.12)

con

$A=1-2Fo/\gamma ((1-{\varphi }_{1})-(1+{\varphi }_{2}-2{\varphi }_{1})cos(kh)-$$({\varphi }_{1}-2{\varphi }_{2})cos(2kh)-{\varphi }_{2}cos(3kh))$

( 5.13)

donde A es un factor de amplificación y para que el sistema sea estable su valor absoluto debe ser menor a 1 para todas las ondas.

Usando ahora la identidad trigonométrica [1-cos(x )]=2sin² (x /2), la fórmula 5.13 se convierte en:

$A=1-4Fo/\gamma ((1+{\varphi }_{2}-2{\varphi }_{1}){sin}^{2}(1/2kh)+$$({\varphi }_{1}-2{\varphi }_{2}){sin}^{2}(kh)+{\varphi }_{2}{sin}^{2}(3/2kh))$

( 5.14)

Buscando el valor de la fase kh que maximice los términos multiplicados por (Fo /γ ) podemos encontrar el límite de estabilidad. Además esto nos permite determinar la longitud de la onda que causa la divergencia. Esta tarea se realizó utilizando Matlab [11] para diferentes radios de integración. Los valores obtenidos mediante este análisis y los Δt_críticos para una malla de 2000 elementos con κ /C =1 se pueden observar en la tabla 1 .

Tabla 1. Δt_crít obtenidos mediante diferentes evaluaciones

Rball	Δtcrít obtenido mediante el análisis de V. Neumann	Δtcrít obtenido en los experim. numéricos	${\textstyle \Delta t_{crit}={\frac {1}{2}}{\frac {{\left(R^{\ast }\right)}^{2}}{\kappa /C}}}$	Norma L2 error	Tiempo de ejecución [s]
<1	0,500	0,50	0,50	3,82E-04	37,69
2	1,333	1,33	1,12	4,59 E-04	18,57
3	2,693	2,69	2,22	5,76 E-04	11,34
4	4,606	4,61	3,78	8,01 E-04	7,93
5	7,069	7,07	5,78	1,10 E-03	5,98
6	10,08	10,09	8,22	1,52 E-03	4,79
7	13,64	13,65	11,11	2,03 E-03	4,02
8	17,75	17,77	14,44	2,41 E-03	3,42
9	22,40	22,42	18,23	3,18 E-03	3,02

5.4. Relación entre el ancho del kernel y la estabilidad

El análisis por Von-Neumann presentado anteriormente resulta adecuado para evaluar la estabilidad del método propuesto y obtener una relación entre el ancho del kernel R_ball y el radio R^* para el cual se cumple la ecuación 3.5.

Para obtener una fórmula rápida, que permita deducir en forma inmediata el análisis del paso de tiempo crítico en función de la ϕ utilizada se ha realizado el siguiente análisis:

Si la función ϕ es lineal, de forma tal que valga uno para x =0 y cero para x =R_ball , entonces es fácil ver que en este caso, la aproximación de g (x ) realizada en (3.2) es equivalente a tomar un elemento de tamaño h =R_ball . Como es bien sabido, el Δt_crit para un elemento de este tamaño es ${\textstyle \Delta t_{crit}={\frac {{R_{ball}}^{2}}{2\kappa /C}}}$ . Este es el paso máximo de estabilidad para una ϕ de radio R_ball . El inconveniente de esta ϕ es que tiene la precisión de un elemento de h =R_ball . Si queremos mejorar la precisión, sin disminuir excesivamente el paso de tiempo, se pueden utilizar funciones ϕ que tengan mayor peso en la cercanías de x =0 y menor peso en x =R_ball . Si llamamos ϕ_L a la función lineal podemos establecer la siguiente relación:

${\frac {\int _{\Omega }\varphi {\mbox{ }}d\Omega }{\int _{\Omega }{\varphi }_{L}{\mbox{ }}d\Omega }}=$${\frac {R^{\ast }}{R_{ball}}}{\mbox{ }}{\mbox{por}}{\mbox{ }}{\mbox{tanto}}{\mbox{ }}R^{\ast }=$$R_{ball}{\frac {\int _{\Omega }\varphi {\mbox{ }}d\Omega }{\int _{\Omega }{\varphi }_{L}{\mbox{ }}d\Omega }}$

( 5.15)

Utilizando la función propuesta en (5.5e) integrándola en un dominio unidimensional de radio R_ball obtenemos:

$R_{1D}^{\ast }=R_{ball}{\frac {\int _{-R_{ball}}^{+R_{ball}}\left(1-2\left|{\frac {x}{R_{ball}}}\right|+{\left|{\frac {x}{R_{ball}}}\right|}^{2}\right)d\Omega }{\int _{-R_{ball}}^{+R_{ball}}\left(1-\left|{\frac {x}{R_{ball}}}\right|\right){\mbox{   }}d\Omega }}=$$R_{ball}{\frac {2/3{\mbox{ }}R_{ball}}{R_{ball}}}=2/3{\mbox{ }}R_{ball}$

( 5.16)

Extrapolando el concepto para problemas 3D y tomando como coordenadas del nodo: [0,0,0] la función φ_L debe cumplir:

${\begin{array}{c}\varphi ([0,0,0])=1\\\varphi (\pm [R_{ball},0,0])=\varphi (\pm [0,R_{ball},0])=\varphi (\pm [0,0,R_{ball}])=0\end{array}}$

( 5.17)

por lo que:

${\varphi }_{L}=1-\left|{\frac {x}{R_{ball}}}\right|-\left|{\frac {y}{R_{ball}}}\right|-$$\left|{\frac {z}{R_{ball}}}\right|$

Tomando como dominio de integración para la φ_L el octaedro definido por los puntos dados en 5.17 y una esfera de radio R_ball para la φ evaluada finalmente se obtiene el R^* para el caso tridimensional como:

$R_{3D}^{\ast }=R_{ball}{\frac {1/3}{2/15{\mbox{ }}\pi }}\approx 0,796{\mbox{ }}R_{ball}$

( 5.18)

Nota: aunque quede fuera del alcance de este trabajo, siguiendo la misma idea anterior, se podrían obtener aproximaciones explícitas de mayores órdenes en el tiempo. Por ejemplo, en vez de pensar en una aceleración g^t (x ) = gⁿ (x ) constante en el tiempo, podemos pensar en una aceleración linealmente variable en el tiempo:

$g^{t}({\boldsymbol {\mbox{x}}})=g^{n}({\boldsymbol {\mbox{x}}})+$${\frac {g^{n}({\boldsymbol {\mbox{x}}})-g^{n-1}({\boldsymbol {\mbox{x}}})}{\Delta t}}{\mbox{ }}t$

( 5.19)

Y por tanto, el método explícito modificado de segundo orden en el tiempo nos queda (Adams–Bashforth):

$T^{n+1}({\boldsymbol {\mbox{x}}})=T^{n}({\boldsymbol {\mbox{x}}})+$$\left[{\frac {3}{2}}{\overset {\mbox{ˆ}}{g}}^{n}({\boldsymbol {\mbox{x}}})-\right.$$\left.{\frac {1}{2}}{\overset {\mbox{ˆ}}{g}}^{n-1}({\boldsymbol {\mbox{x}}})\right]{\mbox{ }}\Delta t$

( 5.20)

6. Ejemplos numéricos

6.1. Ejemplo 1D

Para verificar la consistencia del método propuesto se simuló con Matlab un caso de malla homogénea con Δx =1 y κ /C =1, con 2.000 nodos. Las condiciones iniciales de temperatura fueron dadas por una función Gaussiana. Dado que se conoce la solución analítica para este problema se computó la norma de error L² para cada paso de tiempo y se guardó su valor más alto siempre y cuando la temperatura en los contornos se mantenga cercana a cero. Los resultados se pueden apreciar en la tabla 1 . El tiempo de ejecución es para 100.000 segundos de simulación.

Observando los resultados obtenidos podemos concluir que:

• el costo del método modificado decrece linealmente con el ancho del esténcil empleado.
• esto se explica porque si bien el cálculo es más costoso, el paso de tiempo empleado crece más que linealmente con el ancho del esténcil, diríamos que cuadráticamente, con lo cual el costo final se reduce linealmente. Esto quiere decir que si se utiliza un esténcil de semi-ancho 4, el nodo en cuestión interactúa con 4 vecinos a cada lado, el costo se reduce 4 veces en una malla homogénea.
• la precisión cae al aumentar el radio de integración.
• comparando la columna 2 con la columna 3 se aprecia que el análisis por von Neumann resulta adecuado para evaluar la estabilidad del método propuesto.
• en la columna 4 de la tabla 1 se observa los valores de Δt_crít obtenidos con la aproximación propuesta en (5.15) son aceptables y se puede ver que son conservativos.

6.2. Ejemplos 3D

La generación de una malla 3D por medio de elementos finitos de forma de tetraedros es, posiblemente, la mas utilizada para resolver problemas prácticos. Sin embargo, hasta ahora ha resultado muy difícil obtener mallas de tetraedros que sean luego eficientes para obtener soluciones utilizando métodos explícitos. El motivo es que es casi imposible que toda la malla tenga elementos tetraédricos regulares. Efectivamente, la mejor de las mallas 3D tetraédricas, luego de aplicarle algún algoritmo de eliminación de los slivers (tetraedro achatado con volumen casi 0), presentan algunos elementos con un espesor δ que, en el mejor de los casos son del orden h /10, siendo h el espaciamiento promedio entre los nodos de la malla [12] . Como el número de Fourier necesario para la estabilidad de los métodos de integración explícitos depende de δ² , esto hace que para una malla de tamaño h , para la cual uno espera un límite de integración estable para

$\Delta t<{\frac {h^{2}}{2\kappa /C}}$

( 6.1)

obtenga estabilidad solo para

$\Delta t<{\frac {{\delta }^{2}}{2\kappa /C}}\approx {\frac {h^{2}}{200\kappa /C}}$

( 6.2)

Por este motivo la utilización de un método explicito estándar resulta prohibitivo en la mayoría de los casos. Como veremos en los ejemplos 3D, la modificación propuesta en este artículo permite salvar este inconveniente y obtener pasos de tiempos proporcionales a h² (en vez de ser proporcionales a δ² ). Incluso se pueden obtener pasos de tiempo mayores al correspondiente al Fourier de la malla a precisiones de cálculo muy aceptables.

Las versiones del código en 2 y 3 dimensiones fueron desarrolladas dentro del ambiente de KRATOS [14] . De esta manera se evitó programar las tareas comunes de todos los programas de elementos finitos.

6.2.1. Prisma con malla homogénea y campana de Gauss de temperatura en su interior

Se trata de un prisma de dimensiones 50x50x10 con una malla de tetraedros con una distancia entre nodos más o menos constante igual a ${\textstyle h=1/{\sqrt {3}}}$ , con ${\textstyle \kappa /C=0.95}$ y con contornos adiabáticos. En el interior se introduce un campo de temperaturas inicial en forma de una campana de Gauss radial que se deja difundir en el tiempo (ver figura 2 ). La solución analítica del problema para un dominio de dimensiones infinitas se obtiene mediante la fórmula 6.3. Se tomó como tiempo inicial t =5 y se dejó difundir hasta t =55, de forma tal que los límites del dominio no alteren la solución. (Es decir que la temperatura no ha sido sensiblemente modificada en los contornos).

$T^{t}(r)={\frac {e^{\left({\frac {-r^{2}}{4\kappa /c{\mbox{   }}t}}\right)}}{4\pi {\mbox{   }}t}}$

( 6.3)

Figura 2. Prisma con h constante. Condiciones de temperatura iniciales.

El ejemplo fue procesado usando 3 métodos, a) un método explicito estándar con el máximo paso de tiempo posible para que sea estable, b) usando un método implícito con ${\textstyle \theta =1/2}$ para varios pasos de tiempo diferentes y c) con el método explicito modificado presentado en este artículo para distintos valores del paso de tiempo y usando la función de ponderación definida en (5.5):

$\varphi =1-2\left({\frac {d}{R_{ball}}}\right)+{\left({\frac {d}{R_{ball}}}\right)}^{2}$

( 6.4)

donde R_ball es el radio de integración y d la distancia al nodo central.

La malla utilizada fue de 500.000 tetraedros con una distribución uniforme. La misma fue generada con GID [13] . A pesar de tener una distribución de nodos aproximadamente constante, la malla tiene algunos elementos distorsionados (con ángulos diedros menores a 10°). Aún así es considerada una buena malla para un dominio 3D.

La tabla 2 muestra los resultados obtenidos con los distintos métodos. La columna 1 muestra los radios de integración empleados. La segunda columna el paso de tiempo utilizado. En la tercera y cuarta columna aparecen el error y el tiempo de cálculo para el método explícito modificado (MEM) presentado en este trabajo y en las columnas quinta y sexta los mismos valores para un método implícito (MI) con ${\textstyle \theta =1/2}$ . Para cuantificar la exactitud se calculó la norma de error L² en el dominio para cada paso de tiempo y se retuvo el valor máximo. Errores por debajo de 10^–04 son errores de la discretización espacial para esa malla y por lo tanto no deben tenerse en cuenta para el análisis del error de la integración temporal.

Tabla 2. Prisma con ${\textstyle h=1/{\sqrt {3}}}$ constante. Resultados para diferentes Δt

r	Δt	error MEM	tejec MEM	error MI	tejec MI
<h	0,013	1,142E-04	356	1,14E-04	1.430
0,635	0,055	1,136E-04	94	1,14E-04	358
1,154	0,26	1,161E-04	19,72	1,49E-04	90
1,732	1,0	5,052E-04	5,85	2,31E-04	30
2,309	1,95	1,723E-03	3,61	9,77E-04	18
2,887	2,95	6,701E-03	3,25	1,48E-03	14

El máximo paso de tiempo (Δt ) estable para el caso explícito estándar fue de 0,013 segundos. Este valor es coherente para la malla descripta anteriormente con esos valores de tetraedros deformados.

Se calcularon usando el método implícito y el método explicito modificado pasos de tiempo sucesivamente cada vez más grandes hasta llegar a pasos de tiempo casi 200 veces mas grandes y se analizó para cada caso, el tiempo total de cálculo para los 55 segundos de tiempo de simulación y la precisión alcanzada. Es de hacer notar que todos los cálculos se hicieron en una maquina serial, por tanto se supone que en caso de paralelización, el método explícito debería tener una performance mucho mayor que el método implícito.

En los resultados de la tabla 2 se pueden observar las siguientes características:

• Para el método explícito modificado (MEM) propuesto en este artículo, la solución fue siempre estable, aun para pasos de tiempo casi 200 veces mayores que el paso límite de estabilidad para el explícito estándar (caso con r < h ).
• El tiempo total de cálculo, a igualdad de paso de tiempo, es siempre menor en el MEM que en el método implícito (MI). Esto es menos de 4 segundos para el MEM contra más de 14 s para el MI en el caso del máximo paso de tiempo estudiado.
• A igualdad de paso de tiempo, la precisión del MI es mas alta que el MEM. No obstante los resultados en cuanto a precisión del MEM son razonables y perfectamente aceptables para los límites normales de precisión en cálculos ingenieriles.
• Comparado con el método explícito estandar, la modificación propuesta permite un método explícito con pasos de tiempo 200 veces más grandes en una malla 3D con una distribución de nodos regular.
• El paso de tiempo estable del método explícito estándar se haría aun más pequeño en caso de que la malla tenga muchos elementos defectuosos. Para el MEM los elementos defectuosos no modifican el paso de tiempo estable si se utiliza un radio de integración suficiente grande, como se explica en el punto siguiente.
• La aproximación numérica de las integrales de ϕ y gⁿ propuesta en (3.18) son válidas para mallas regulares sin elementos distorsionados. En el caso de mallas irregulares, con grandes variaciones de la distancia entre nodos o en el caso de mallas regulares, pero con elementos distorsionados, la aproximación (3.18) deja de ser una integración sobre un radio constante R_ball para pasar a ser una integración sobre un dominio irregular con un radio muy variable dependiendo de la distorsión de los elementos. Por este motivo, la ecuación (3.6) que demuestra ser apropiada para problemas 1D y 2D con mallas regulares, se puede interpretar como imprecisa para mallas 3D en presencia de elementos distorsionados. Esta imprecisión es menor a medida que aumentamos el radio de integración. El motivo es que para radios pequeños (de tamaño h o 2h ) se toman pocos nodos de la vecindad del nodo x_i y por lo tanto la influencia de los elementos distorsionados es grande. Para radios grandes la influencia de los elementos distorsionados se hace menos sensible ya que se comportan como macro-elementos generados por la nube de nodos dentro del radio de integración. Este comportamiento se puede apreciar en la figura 3 donde se ha trazado el Δt estimado en una malla no distorsionada frente al Δt realmente obtenido para la malla 3D en estudio.

Figura 3. Mínimo R* requerido para que el método sea estable versus Δt . a) para una malla regular (gris) y b) para la malla del ejemplo (negro).

• En la figura 4 se puede observar que para obtener un error máximo determinado independientemente del Δt , el MEM siempre resuelve el problema para todo el intervalo en menos tiempo.

Figura 4. Tiempo de cálculo requerido para obtener un determinado error: a) Método explícito (negro) y b) método implícito (gris).

6.2.2. Prisma con malla variable

Para evaluar la capacidad del MEM de resolver problemas que poseen mallas de h no constante, se remalló el dominio del problema anterior utilizando 2 tamaños diferentes como se aprecia en la figura 5 . Con esto se busca captar con mayor precisión la zona donde el gradiente de temperatura es más elevado.

Figura 5. Malla con 2 tamaños de elementos: h1 =0,3 y h2 =0,9.

La tabla 3 muestra los resultados obtenidos para diferentes Δt por diferentes métodos.

Tabla 3. Prisma con h variable. Resultados para diferentes Δt por distintos métodos

r	dt	tejec MEM	error MEM	tejec MI	error MI
< h1	0,0032	1441	9,256E-05	21.616	9,241E-05
0,36	0,0067	711	1,471E-04	10.098	9,241E-05
0,6	0,021	231	2,879E-04	2.891	9,243E-05
0,9	0,033	194	6,284E-04	1.721	9,972E-05
1,2	0,072	113	9,078E-04	676	1,130E-04
1,5	0,15	68,1	9,617E-04	322	1,153E-04
1,8	0,29	52,7	1,038E-03	190	1,693E-04
2,1	0,51	41,1	1,656E-03	128	1,837E-04
2,4	0,82	34,9	2,215E-03	95	2,338E-04
2,7	1,25	30,1	3,148E-03	68	5,391E-04

Al igual que en la tabla 2 , la primera fila corresponde al caso explicito estándar (r <0), que para este ejemplo tiene un Δt máximo de estabilidad de 0,0032 s y tarda en analizar los primeros 55 s 1.441 s de cálculo.

La primer columna muestra los diferentes radios de integración utilizados y la columna 2 los diferentes pasos de tiempo estables para ese radio de integración. Se puede apreciar que se logran pasos de tiempo hasta 390 veces más grandes que el límite de estabilidad estándar con una ganancia en el tiempo total de cálculo de 47 veces. El error se deteriora, pero se mantiene siempre dentro de los límites ingenieriles aceptables. Los resultados del método implícito (MI) en las columnas 5 y 6 de la tabla, muestran que MI da siempre errores inferiores a los errores espaciales, pero con tiempos de cálculo siempre superiores a los obtenidos con MEM.

En este problema es interesante observar que para obtener un mismo Δt se debieron utilizar radios de integración mayores que en el caso anterior de h =constante. Esto se debe a que en la zona donde la malla es más gruesa (h₂ =0,9) el coeficiente entre el radio de integración y h es más bajo (ver figura 3 ) y por lo tanto la integración pierde precisión, como se explicó en el punto 6 del apartado anterior.

6.2.3. Cubos con salto brusco del coeficiente de conductividad en su interior

Con este ejemplo (ver figura 6 ) se buscó determinar la degradación de la solución debido a saltos bruscos en el coeficiente de conductividad entre un elemento y otro El modelo son 2 cubos contiguos de 1x1x1 con una malla de h= 0,07. El primero con un κ /C =1,0 y el segundo de κ /C =0,1. Se fijan como condiciones de contorno T=1 en un extremo y T=0 en el otro. El resto de las caras se consideran adiabáticas.

Figura 6. Dominio del problema: el cubo derecho posee elementos con κ /C =1,0 mientras que el izquierdo κ /C =0,1.

Por analogía eléctrica la temperatura en la sección central debe ser de 1/1,1= 0,909 y hacia ambos lados variar con un gradiente constante hasta llegar a la temperatura de los contornos. La figura 7 muestra los resultados obtenidos con el MEM y el MI, usando en ambos casos un Δt =0,0147 y r =5h en el MEM. Es de hacer notar que Δt =0,0147 corresponde a un paso de tiempo 100 veces mayor que el Δt estable en un método explícito estándar.

Figura 7. Prueba con 2 materiales diferentes: Distribución de temperatura al llegar a la solución estacionaria.

El método explicito propuesto en este trabajo, no presentó ninguna variación en cuanto a la estabilidad ya expuesta anteriormente para el caso de mallas regulares. No obstante, al tomar información en un cierto radio, suaviza levemente la transición entre ambos materiales como se puede observar en la figura 7 .

6.2.4. Cubos con salto brusco de temperatura en su interior

Este ejemplo se utilizó para ver el comportamiento del MEM para casos con grandes gradientes de la función incógnita. La geometría analizada es la misma del ejemplo anterior utilizando ahora un coeficiente de conductividad constante en todo el dominio y condiciones de contorno y condiciones iniciales:

${\Gamma }_{D}={\begin{array}{c}T=1{\mbox{ }}para{\mbox{ }}x=0\\T=0{\mbox{ }}para{\mbox{ }}x=1\end{array}}{\mbox{ }}T_{t=0}={\begin{array}{c}T=1{\mbox{ }}para{\mbox{ }}x<1\\T=0{\mbox{ }}para{\mbox{ }}x>=1\end{array}}$

El problema fue resuelto en forma explicita mediante el MEM con un paso de tiempo mucho mayor al límite de estabilidad para el método explícito estándar y los resultados no presentaron ningún tipo de inestabilidad. Para apreciar la precisión del método, el mismo problema fue resuelto en forma implícita (MI). Como es sabido el método implícito con θ =1/2 es incondicionalmente estable. No obstante en este problema con grandes gradientes de temperatura, los métodos implícitos presentan al comienzo del período de tiempo, importantes oscilaciones espurias cuando se utilizan grandes pasos de tiempo. Si bien estas oscilaciones pueden amortiguarse utilizando θ> 0,8, esto conlleva a un aumento del error (por acercarnos al menos exacto Backward-Euler) y por lo tanto se necesitaría adaptar este factor en cada problema de forma tal de obtener la mayor precisión posible sin oscilaciones.

La figura 8 muestra la evolución de la temperatura en un nodo que se encuentra en x=1 comparando los 2 métodos utilizando un Δt tal que Fo =5 Además se traza como referencia la solución brindada por el MI utilizando un Δt pequeño (Fo<<1 ). Se puede apreciar que la convergencia del MEM a la solución estacionaria es más lenta que el MI, pero con la ventaja de no presentar oscilaciones espurias y ser completamente explícito. La fuerte discontinuidad que presenta el problema hace que el método implícito necesite al menos un θ ≥0,8 en lugar del valor teórico de θ ≥0,5

Figura 8. Prueba con grandes gradientes de temperatura. Evolución de la temperatura para un nodo a una distancia x=1 para el MEM y el MI a Fo =5 comparada con una solución con Fo<<1.

6.2.5. Disipador 3D

Por último se calcularon los resultados sobre una aplicación ingenieril del problema de transmisión del calor. Se trata de un modelo de disipador del calor cuyas dimensiones son 10 x 10 x 10 y posee las siguientes propiedades: conductividad κ /C =1.000. El mismo es sometido a una fuente de calor variable en la base, que oscila en intensidad con una función sinusoidal. El dominio incluye además una región con conductividad menor (κ /C =10) como se aprecia en la figura 9 . En los contornos de este dominio se fijó la temperatura en T=20°C .

Figura 9. Disipador 3D. Dominio del problema.

Dado que no se tiene una solución analítica para este problema, se tomaron como resultados de comparación los obtenidos con un método implícito (MI) con un paso de tiempo correspondiente a Fo= 0,07, que es el máximo paso de tiempo estable para el esquema explícito estándar (ME). Se compararon estos resultados con los brindados por el MEM para diferentes radios de integración y con el MI para los mismos Δt . Los tiempos de ejecución para un segundo de simulación y los errores utilizando la norma L₂ , se pueden ver en la tabla 4 . El error indicado en la tabla es el máximo error obtenido en un paso de tiempo durante el intervalo de cálculo.

Tabla 4. Disipador 3D. Resultados para diferentes Δt por distintos métodos

r	Δt	Error MEM	tejec MEM	Error MI	tejec MI
<1 (ME)	0,0000055	3,15E-06	8759	-	-
2,2	0,0000111	4,9E-04	4530	3,65E-06	19379
3	0,000050	1,12E-03	1101	3,26E-05	6451
4	0,000189	2,35E-03	281	1,12E-04	2747
5	0,000444	4,03E-03	178	3,35E-04	1684
6	0,000689	6,0E-03	160	4,5E-04	1141

Como se aprecia en la tabla 4 , con el MEM se obtuvieron errores siempre menores al 1%, requiriendo un tiempo de ejecución total mucho menor al del método explícito estándar utilizando el Δt_crít o el MI utilizando el mismo paso de tiempo. Esta ventaja del MEM en cuanto a tiempo de cálculo sería aún mayor si se sacara provecho del mejor speed-up que gozan los métodos explícitos frente a los implícitos en cuanto a la eficiencia de la paralelización.

7. Conclusiones

Se ha presentado una modificación a introducir en la integración temporal explícita de las ecuaciones de transmisión del calor por conducción con las siguientes características:

• Permite obtener soluciones transitorias sin necesidad de invertir ninguna matriz
• Es estable para pasos de tiempo mucho mayores que los pasos de tiempo usados en los métodos explícitos estándar.
• Consume mucho menos tiempo de cálculo que un método explícito estándar y a igual tamaño de paso de tiempo también consume menos tiempo de cálculo que un método implícito.
• Por ser un método explícito es fácil de programar (autómata celulares) en maquinas con múltiples procesadores y por la misma razón se espera una alta eficiencia en la paralelización.
• Para grandes pasos de tiempo, tiene menos precisión que un método implícito, pero es más eficiente en cuanto a tiempo de cálculo.
• El método es simple de programar y/o de adaptar a programas existentes.

El método es fundamentalmente útil en problemas 3D donde la utilización de una integración explícita estándar es muy costosa debido a que la estabilidad depende del elemento más distorsionado. Con este método se recupera la posibilidad de regular el paso de tiempo deseado en función del error requerido y no en función de la estabilidad del método.

Sin embargo debe tenerse en cuenta que el tamaño de la matriz Φ aumenta según el cubo del radio de integración, por lo que grandes radios de integración (del orden de 10h o más) tornan al método menos eficiente debido a la cantidad de operaciones que requiere y sobre todo a la memoria requerida para almacenar dicha matriz; a partir de ciertos valores de r el tiempo de cálculo en el MEM es tal que no compensa la ganancia por aumentar el paso de tiempo.

En conjunto con el método X-IVAS propuesto en la referencia [4] permitiría tratar en forma explícita problemas donde hay también convección, sin límites de estabilidad.

Su generalización para la resolución de otras ecuaciones similares como las ecuaciones de Navier-Stokes parece evidente y será probada en futuros trabajos.

Agradecimientos

Este trabajo fue parcialmente financiado por el Advanced Grant: ERC-2009 -AdG Real Time Computational Mechanics Techniques for Multi-Fluid Problems otorgado por la European Research Council . Norberto Nigro desea agradecer al CONICET, la Universidad Nacional del Litoral y a la Agencia Nacional de Promoción Científica y Tecnológica (ANPCyT) por su apoyo económico a través de los proyectos PICT 1645 BID (2008) y CAI+D 65-333 (2009)). Un especial agradecimiento al Dr.Pedro Morín por sus invalorables comentarios en los aspecto matemáticos del método aquí desarrollado.

Bibliografía

1. [1] S.R. Idelsohn, E. Oñate, F. Del Pin; The particle finite element method: a powerful tool to solve incompressible flows with free-surfaces and breaking waves; Int J Numer Meth Eng., 61 (2004), pp. 964–984
2. [2] F. Mossaiby, R. Rossi, P. Dadvand, S.R. Idelsohn; OpenCL-based implementation of an untructured edge-based Finite Element convectio-diffusion solver on graphic hardware; Int J Numer Meth Eng. (2011) http://dx.doi.org/10.1002/nme.3302
3. [3] R. Rossi, J. Cotela, N.M. Lafontaine, P. Dadvand, S.R. Idelsohn; Parallel adaptive mesh refinement for incompressible flow problems In Press, Corrected Proof, Available on line 8 February 2012; Comput Fluids (2012)
4. [4] S. Idelsohn, N. Nigro, A. Limache, E. Oñate; Large time-step explicit integration method for solving problem with dominant convection; Comput. Method. Appl. M., 217–220 (2012), pp. 168–185 http://dx.doi.org/10.1016/j.cma.2011.12.008
5. [5] N. Satofuka; A new explicit method for the numérical solution of parabolic differntial equations; Numerical Properties and Methodologies in Heat Transfer (A84-40651 19-34), Hemisphere Publishing Corp, Washington, DC (1983) p. 97-108
6. [6] J. Teixeira; Stable schemes for partial differential equations: the one-dimensional diffusion equation; J. Comput. Phys., 153 (1999), pp. 401–417
7. [7] C. Girard, Y. Delage; Stable schemes for the vertical diffusion in atmospheric circulation models; Mon Weather Rev, 118 (1990), p. 737
8. [8] E. Hairer; Unconditionally stable explicit methods for parabolic equations; Numer. Math, 35 (1980), pp. 57–68
9. [9] Norberto Nigro, Juan Gimenez, Alejandro Limache, Sergio Idelsohn, Eugenio Oñate, Nestor Calvo; A new approach to solve Incompressible Navier-Stokes Equations using a Particle Method; Mec. Comp., XXX (6) (2011), pp. 451–483
10. [10] N. Nigro, S. Idelsohn, J. Gimenez, A. Limache, P. Novara, N. Calvo, et al.; Searching efficient Navier-Stokes solver using a particle method; TCCM 2011, Padova (Italia), 12-14 setiembre (2011)
11. [11] Matlab. Numerical computing environment. Disponible en: http://www.mathworks.com/ [consultado 4 Jul 2012].
12. [12] S.R. Idelsohn, E. Oñate; To mesh or not to mesh. That is the question…; Comput. Method. Appl. M., 195 (37–40) (2006), pp. 4681–4696
13. [13] GID. The Personal Pre and Post Processor. Disponible en: http://gid.cimne.upc.es/ [consultado 4 Jul 2012].
14. [14] P. Dadvand, R. Rossi, S.E. Oñate; An object-oriented environment for developing finite element codes for multi-disciplinary applications; Arch Comput Methods Eng, 17 (2010), pp. 253–297