## Abstract

In this paper we have developed the numerical solution of two problems of diffusion-convection (DC), using the finite element method of Streamline Upwind Petrov-Galerkin (SUPG). The parameters that define the behavior of the equations are modeled as stochastic fields, specifically, are used: the convective velocity, diffusion and heat capacity as values of random type. Therefore, we have included SUPG method to DC, with dominant convection, with the stochastic spectral finite element method. Each parameter was described by Karhunen-Loève expansion, while the unknown is represented by the polynomial expansion of the chaos. The objectives of the paper are: 1. To study the influence of stochastic fields in solving problems with SUPG DC and 2. Get the solution of each expansion unknown variable. The results show the versatility of the method for solving different physical problems due to the generality of the statistical description of the solution and the richness in the representation of the areas where there is the greater variability in response. The patterns shown in the unknown uncertainty depends on the dynamics of diffusion, convective velocity and the type of solution used.

## Palabras clave

Petrov-Galerkin ; Elementos finitos estocásticos ; Streamline Upwind Petrov Galerkin ; Advección ; Difusión

## Keywords

Petrov-Galerkin ; Stochastic finite element ; Streamline Upwind Petrov Galerkin ; Advection ; Diffusion

## 1. Introducción

Muchos problemas de la física, la química, la economía, la biología, la bioingeniería, e incluso la ecología, entre otros campos, pueden ser modelados a través del balance de 3 fenómenos: la difusión, la convección y la reacción [1] , [2] , [3] , [4] , [5] , [6] , [7] , [8] , [9] , [10] , [11] , [12]  and [13] . Dicho balance se plantea en la ecuación diferencial (1) de difusión convección-reacción, y se complementa con las condiciones de contorno descritas en (2), (3) y (4).

 ${\displaystyle {\frac {\partial \varphi }{\partial t}}-\nabla \cdot ({\boldsymbol {\mbox{k}}}\nabla \varphi )+}$${\displaystyle {\boldsymbol {\mbox{u}}}\cdot \nabla \varphi -s\varphi =}$${\displaystyle f({\boldsymbol {\mbox{x}}}){\mbox{ }}en{\mbox{ }}\Omega }$
( 1)
 ${\displaystyle \varphi (t=0)={\varphi }_{o}}$
( 2)
 ${\displaystyle \varphi ({\boldsymbol {\mbox{x}}})=g({\boldsymbol {\mbox{x}}}){\mbox{ }}{\mbox{en}}{\mbox{ }}{\Gamma }_{\varphi }}$
( 3)
 ${\displaystyle {\overset {\rightarrow }{\nabla }}(\varphi )\bullet {\boldsymbol {\mbox{n}}}=}$${\displaystyle h({\boldsymbol {\mbox{x}}}){\mbox{ }}{\mbox{en}}{\mbox{ }}{\Gamma }_{\nabla }}$
( 4)

donde k es la matriz de coeficientes difusivos, u es el campo de velocidad asociado al proceso convectivo, s es el coeficiente fuente (s  > 0 implica producción y s  < 0 significa disipación), ${\textstyle f({\boldsymbol {\mbox{x}}})}$ es la función de generación, ${\textstyle g({\boldsymbol {\mbox{x}}})}$ es la función que define el valor del campo escalar ${\textstyle \varphi }$ sobre la frontera ${\textstyle {\Gamma }_{\varphi }}$ y ${\textstyle h({\boldsymbol {\mbox{x}}})}$ es la función que define el valor del flujo sobre la frontera ${\textstyle {\Gamma }_{\nabla }}$ , cuya normal es n .

El comportamiento y los patrones de la solución de los ejemplos señalados [1] , [2] , [3] , [4] , [5] , [6] , [7] , [8] , [9] , [10] , [11] , [12]  and [13] dependen de los parámetros físicos del problema. En (1) los parámetros que definen el comportamiento de la ecuación son k , u y ${\textstyle s}$ . En casos más generales, en ecuaciones de reacción-difusión-convección, donde ${\textstyle s}$ puede depender de más de una variable, del espacio y del tiempo, se puede obtener un conjunto de relaciones entre los parámetros que definen el comportamiento de la solución, por ejemplo, se pueden obtener patrones de Turing o de Hopf [2]  and [5] . Por tanto, los parámetros definen el tipo de solución de la ecuacion y, además, definen el tipo de método numérico que se debe utilizar. En efecto, aunque la formulación convencional de elementos finitos, o método de Bubnov-Galerkin, resulta útil y adecuada para el tratamiento de muchos problemas de la ingeniería, especialmente en el campo de la mecánica de sólidos, presenta problemas de estabilidad cuando en la ecuación diferencial aparecen operadores no autoadjuntos, tal como el término convectivo en la ecuación (1). Este término, bajo la formulación convencional (Bubnov-Galerkin), tiene un efecto desestabilizador de la matriz de rigidez, ya que introduce asimetría en esta y produce oscilaciones falsas en la aproximación. La eliminación de estas oscilaciones, o estabilización de la solución, se logra regresando el carácter simétrico a la matriz de rigidez, lo cual se puede alcanzar mediante diversos métodos, entre los cuales cabe citar: el método de las líneas de características [14] , el método de cálculo finito [15] , el método de paso fraccional ${\textstyle \theta }$[16] , el método de mínimos cuadrados de Galerkin [17] y el método Petrov-Galerkin de contracorriente [18] , [19] , [20] , [21] , [22] , [23] , [24] , [25] , [26] , [27] , [28]  and [29] .

De otro lado, los parámetros mencionados pueden ser deterministas o estocásticos. Cuando son deterministas el método de Petrov-Galerkin es adecuado para encontrar la solución aproximada [18] . En el caso de que los parámetros sean estocásticos (variables de campo aleatorio), no existen, hasta donde los autores conocen, reportes sobre la solución de ecuaciones de difusión-convección (con convección dominante) que permitan incluir la variabilidad de campo en los parámetros de la ecuación (1). Para solucionar el problema estocástico se han utilizado varios métodos, entre los que se cuentan el método de la perturbación y los métodos de elementos finitos estocásticos espectrales [30] .

Alternativamente al método de las perturbaciones, en los años 90 se introdujo el método de los elementos finitos estocásticos espectrales, los cuales solucionan ecuaciones diferenciales parciales y ordinarias con parámetros y condiciones de contorno aleatorias [30] . El método fue originalmente desarrollado para problemas de elasticidad lineal en los que las propiedades del material presentaban campos aleatorios gaussianos. Luego, utilizando el mismo principio se extendió a problemas de transferencia de calor [41]  and [42] y campos lognormales [41] . Desarrollos recientes han mostrado la versatilidad del método para diversas aplicaciones, entre las que se incluyen: mecánica de fluidos [43] , análisis modal [44] , interacción fluido estructura [45] , análisis de placas laminadas [46] y optimización [47] . El método consiste en representar el campo aleatorio (que, en general, son las propiedades del material) en términos de la descomposición de Karhunen-Loeve. Posteriormente, la variable que se debe encontrar, que también tiene características estocásticas, se representa mediante los polinomios de caos. Ambas representaciones se discretizan y se implementan en el método de los elementos finitos, para lo cual se utiliza la proyección de Galerkin [30] , [31] , [32] , [33] , [34] , [35] , [36] , [37] , [38] , [39] , [40] , [41] , [42] , [43] , [44] , [45] , [46]  and [47] . Este procedimiento se puede programar en un paquete estándar de elementos finitos, lo cual representa una ventaja económica en términos computacionales con respecto a otras técnicas como los métodos de Monte Carlo [30] y el método de la perturbación.

Pueden encontrarse algunas técnicas alternativas a las ya explicadas en la literatura. Por ejemplo, el método de los elementos finitos estocásticos difusos se ha utilizado para modelar problemas en mecánica de fluidos y sólidos [48] y transferencia de calor [49] , entre otros. Este método introduce en las variables estocásticas el uso de conjuntos difusos que describen su comportamiento. De igual manera, existen otras tantas alternativas para la descripción estocástica, como el método de los elementos finitos extendidos estocásticos [50] , el método de los elementos finitos estocásticos dinámicos de Neumann [51] , o el método de elementos finitos estocásticos multiescala [52] .

Todos los métodos aquí mencionados se han comparado típicamente con el método de Monte Carlo [30] . Este método requiere un exhaustivo cálculo con una gran población de valores aleatorios que simulan el espectro de eventos de un parámetro o una variable. La principal desventaja, motivo por el que se han desarrollado otros métodos, radica en el alto coste computacional, por tanto, los métodos alternos, principalmente el método de los elementos finitos estocásticos espectrales, han sido una buena alternativa.

Los métodos estocásticos han sido aplicados al diseño y al análisis en ingeniería con el objeto de representar fenómenos físicos de forma más precisa y fiable [30] . Esta metodología, en conjunto con las nuevas tecnologías computacionales, provee información de los fenómenos físicos en cuanto a los límites de desempeño real, las variaciones y sus implicaciones en la respuesta de cada fenómeno. En diseño en ingeniería, por ejemplo, el método de los elementos finitos estocásticos permite obtener modelos robustos sobre el comportamiento estructural y el desempeño de máquinas. El proceso de análisis estocástico provee ventajas sobre los métodos de diseño y análisis porque se pueden obtener valores medios, varianzas e intervalos de confianza que aportan mayor información sobre el comportamiento de una estructura, máquina o fenómeno físico. Adicionalmente, el análisis estocástico puede ayudar a desarrollar una guía inicial para identificar aquellas zonas en las cuales se debe hacer especial énfasis en el estudio de un fenómeno físico [53] .

Desde esta perspectiva, el presente artículo desarrolla numéricamente 2 problemas de difusión-convección (DC), empleando el método de Petrov-Galerkin en contracorriente (SUPG). Los parámetros de la ecuación se modelan como variables de campo aleatorio. Así pues, se combina el método SUPG y el método de los elementos finitos estocásticos espectrales. Los objetivos del artículo son: por un lado, estudiar la influencia de los campos estocásticos en la solución de problemas de DC con SUPG y, por otro, obtener los patrones de cada coeficiente del polinomio de caos. Aunque la propuesta teórica se ha desarrollado incluyendo el término reactivo y de fuente, los ejemplos aquí solucionados no consideran estos términos. Sin embargo, los resultados muestran que la incógnita se ve afectada por los campos estocásticos y que el patrón obtenido en los términos de la expansión con el polinomio de caos muestra la incertidumbre en la incógnita, que depende de la dinámica de la difusión, la velocidad convectiva y el tipo de solución utilizado.

## 2. Materiales y métodos

### 2.1. El método de estabilización de Petrov-Galerkin en contracorriente

El SUPG se basa en la estabilización de la matriz de rigidez mediante el empleo de funciones de ponderación modificadas, de modo que se otorgue un mayor peso a la información de los nodos ubicados aguas arriba. Este método de estabilización fue planteado por primera vez por Hughes y Brooks [21] .

En problemas unidimensionales la estabilización con el método de Petov-Galerkin se puede alcanzar añadiendo un término perturbador a la función de peso estándar (${\textstyle W_{i}}$ ), tal como se muestra en (5) [6] , donde h es el tamaño característico del elemento y α es un parámetro de perturbación positivo calculado por medio de la expresión (6) [7] :

 ${\displaystyle W_{i}^{\ast }=N_{i}+{\frac {\alpha h}{2}}{\frac {dN_{i}}{dx}}}$
( 5)

donde:

 ${\displaystyle \alpha =min\left(Coth\left|P_{e}\right|-{\frac {1}{\left|P_{e}\right|}}{\mbox{ }},{\mbox{ }}1-\right.}$${\displaystyle \left.{\frac {1}{\left|P_{e}\right|}}\right)}$
( 6)

En esta última expresión ${\textstyle P_{e}}$ es el número adimensional de Peclet definido como ${\textstyle P_{e}={\frac {\left\|{\boldsymbol {u}}\right\|h}{2k}}}$ . Nótese que la función de peso solo debe ser perturbada si ${\textstyle P_{e}>1}$ . En otras palabras, para valores inferiores a la unidad del número de Peclet, la formulación SUPG lleva implícito el método de Bubnov-Galerkin.

Para una ecuación de difusión-convección, los términos de la matriz de rigidez ${\textstyle \left[K\right]}$ de un elemento e , calculados de acuerdo con el método SUPG unidimensional, están definidos por medio de la expresión (7) [14] :

 ${\displaystyle K_{lm}^{e}=\int _{e}{\frac {dW_{l}^{\ast }}{dx}}k_{im}{\frac {dN_{m}}{dx}}dx+}$${\displaystyle \int _{e}W_{l}^{\ast }u{\frac {dN_{m}}{dx}}dx}$
( 7)

En esta última expresión, Nm son las funciones de forma empleadas para aproximar la solución [14] , tal como se observa en (8):

 ${\displaystyle {\varphi }^{e}=\sum _{m}N_{m}{\varphi }_{m}}$
( 8)

De forma análoga, los términos del vector ${\textstyle \left[F\right]}$ para un elemento ${\textstyle e}$ se definen de acuerdo con (9) [14] :

 ${\displaystyle F_{l}^{e}=\int _{e}W_{l}^{\ast }f(x)dx}$
( 9)

Adicionalmente, este método resulta especialmente útil para problemas en 2 o 3 dimensiones. Por tanto, de forma análoga, se puede perturbar la función de ponderación para estabilizar la convección para el caso bidimensional, como se observa en (10) [21] :

 ${\displaystyle W_{i}^{\ast }=N_{i}+{\frac {\alpha h}{2}}\left[{\frac {u_{x}}{\left\|{\boldsymbol {\mbox{u}}}\right\|}}\left({\frac {\partial N_{i}}{\partial x}}\right)+\right.}$${\displaystyle \left.{\frac {u_{y}}{\left\|{\boldsymbol {\mbox{u}}}\right\|}}\left({\frac {\partial N_{i}}{\partial y}}\right)\right]}$
( 10)

donde ${\textstyle u_{x}}$ y ${\textstyle u_{y}}$ son las componentes de la velocidad en las direcciones globales ${\textstyle x-y}$ , en tanto que α es el coeficiente de perturbación calculado de acuerdo con (6), y ${\textstyle h}$ es la longitud característica del elemento, como se muestra en la figura 1 .

 Figura 1. Línea de corriente al interior de un elemento bidimensional.

### 2.2. Implementación del método de los elementos finitos con el método de estabilización de Petrov-Galerkin en contracorriente

Al aplicar el método de los residuos ponderados a la ecuación diferencial de difusión-convección (1) (haciendo s  = 0), y después de debilitar el término difusivo, se obtiene la expresión (11):

 ${\displaystyle \int _{\Omega }{W_{l}}^{\ast }{\frac {\partial \phi }{\partial t}}d\Omega +}$${\displaystyle \int _{\Omega }{\overset {\rightarrow }{\nabla }}{W_{l}}^{\ast }\cdot {\boldsymbol {\mbox{k}}}{\overset {\rightarrow }{\nabla }}\phi d\Omega +}$${\displaystyle \int _{\Omega }{W_{l}}^{\ast }{\boldsymbol {\mbox{u}}}\cdot {\overset {\rightarrow }{\nabla }}\phi d\Omega =}$${\displaystyle \int _{\Omega }{W_{l}}^{\ast }fd\Omega +\int _{\Gamma }{W_{l}}^{\ast }\left({\boldsymbol {\mbox{k}}}{\overset {\rightarrow }{\nabla }}\phi \right)\bullet {\overset {\rightarrow }{n}}d\Gamma }$
( 11)

donde ${\textstyle l=1,...,nnod}$ , ${\textstyle {\overset {\rightarrow }{n}}}$ es el vector normal al borde de flujo o borde de Neumman y ${\textstyle nnod}$ es el número de nodos del dominio discretizado. Empleando ahora aproximaciones por tramos del tipo ${\textstyle {\varphi }^{e}=\sum _{m=1}^{nnod}N_{m}{\varphi }_{m}}$ , así como la función de peso modificada (10), se puede conducir de la ecuación global a una expresión a nivel elemental [14] , como la mostrada en (12):

 ${\displaystyle {\boldsymbol {\mbox{M}}}{\frac {\Delta {\boldsymbol {\phi }}}{\Delta t}}+}$${\displaystyle {\boldsymbol {\mbox{K}}}{\boldsymbol {\phi }}+{\boldsymbol {\mbox{C}}}{\boldsymbol {\phi }}=}$${\displaystyle {\boldsymbol {\mbox{F}}}}$
( 12)

donde:

 ${\displaystyle {\begin{array}{l}{\boldsymbol {\mbox{M}}}=\int _{{\Omega }^{e}}{\left({\boldsymbol {\mbox{W}}}^{\boldsymbol {\ast }}\right)}^{T}{\boldsymbol {\mbox{N}}}d{\Omega }^{e}\\{\boldsymbol {\mbox{K}}}=\int _{{\Omega }^{\boldsymbol {e}}}{\left({\overset {\rightarrow }{\nabla }}{\boldsymbol {\mbox{W}}}^{\boldsymbol {\ast }}\right)}^{T}\cdot {\boldsymbol {\mbox{k}}}{\overset {\rightarrow }{\nabla }}{\boldsymbol {\mbox{N}}}d{\Omega }^{e}\\{\boldsymbol {\mbox{C}}}=\int _{{\Omega }^{\boldsymbol {e}}}{\left({\boldsymbol {\mbox{W}}}^{\boldsymbol {\ast }}\right)}^{T}{\boldsymbol {\mbox{u}}}\cdot {\overset {\rightarrow }{\nabla }}{\boldsymbol {\mbox{N}}}d{\Omega }^{e}\\{\boldsymbol {\mbox{F}}}=\int _{{\Omega }^{e}}{\left({\boldsymbol {\mbox{W}}}^{\boldsymbol {\ast }}\right)}^{T}fd{\Omega }^{e}+c.c.\\\Delta {\boldsymbol {\phi }}={\boldsymbol {\phi }}^{t+\Delta t}-{\boldsymbol {\phi }}^{t}\\{\boldsymbol {\phi }}={\boldsymbol {\phi }}^{t+\Delta t}\end{array}}}$

donde se ha utilizado una integración temporal del tipo Backward-Euler [15] . Además, a nivel elemental se tienen los siguientes arreglos:

 ${\displaystyle {\begin{array}{l}{\boldsymbol {\mbox{W}}}^{\boldsymbol {\ast }}=\left[{\begin{array}{cccc}{W_{1}}^{\ast }&{W_{2}}^{\ast }&\cdots &{W_{nnodel}}^{\ast }\end{array}}\right]\\{\boldsymbol {\mbox{N}}}=\left[{\begin{array}{cccc}N_{1}&N_{2}&\cdots &N_{nnodel}\end{array}}\right]\\{\boldsymbol {\mbox{W}}}^{\boldsymbol {\ast }}=\left[{\begin{array}{c}{\begin{array}{cccc}{\frac {\partial {W_{1}}^{\ast }}{\partial x}}&{\frac {\partial {W_{2}}^{\ast }}{\partial x}}&\cdots &{\frac {\partial {W_{nnodel}}^{\ast }}{\partial x}}\end{array}}\\{\begin{array}{cccc}{\frac {\partial {W_{1}}^{\ast }}{\partial y}}&{\frac {\partial {W_{2}}^{\ast }}{\partial y}}&\cdots &{\frac {\partial {W_{nnodel}}^{\ast }}{\partial y}}\end{array}}\end{array}}\right]\\{\boldsymbol {\mbox{N}}}^{\boldsymbol {\ast }}=\left[{\begin{array}{c}{\begin{array}{cccc}{\frac {\partial N_{1}}{\partial x}}&{\frac {\partial N_{2}}{\partial x}}&\cdots &{\frac {\partial N_{nnodel}}{\partial x}}\end{array}}\\{\begin{array}{cccc}{\frac {\partial N_{1}}{\partial y}}&{\frac {\partial N_{2}}{\partial y}}&\cdots &{\frac {\partial N_{nnodel}}{\partial y}}\end{array}}\end{array}}\right]\end{array}}}$

donde nnodel es el número de nodos por elemento. Luego se lleva a cabo el ensamble global del conjunto de ecuaciones (12) para la solución numérica del problema completo.

A continuación se desarrolla, alternamente, el método de los elementos finitos espectrales en conjunto con el método SUPG.

### 2.3. El método de los elementos finitos espectrales

En el caso de que los parámetros de la ecuación se puedan representar como campos aleatorios, la ecuación (11) se convierte en:

 ${\displaystyle \left[\int _{{\Omega }^{e}}{\left({\boldsymbol {\mbox{W}}}^{\boldsymbol {\ast }}\right)}^{T}1({\boldsymbol {\mbox{x}}},\omega ){\boldsymbol {\mbox{N}}}d\Omega \right]{\frac {\Delta {\boldsymbol {\phi }}}{\Delta t}}+}$${\displaystyle \left[\int _{{\Omega }^{e}}{\left(\nabla {\boldsymbol {\mbox{W}}}^{\boldsymbol {\ast }}\right)}^{T}k({\boldsymbol {\mbox{x}}},\omega )\nabla {\boldsymbol {\mbox{N}}}d\Omega \right]{\boldsymbol {\phi }}+}$${\displaystyle \left[\int _{{\Omega }^{e}}\left({\boldsymbol {\mbox{W}}}\right)^{T}{\boldsymbol {\mbox{u}}}({\boldsymbol {\mbox{x}}},\omega )\bullet \nabla {\boldsymbol {\mbox{N}}}d\Omega \right]{\boldsymbol {\phi }}=}$${\displaystyle {\boldsymbol {\mbox{F}}}(x,\omega )}$
( 13)

donde los parámetros del modelo (1, ${\textstyle k}$ , ${\textstyle {\overset {\rightarrow }{u}}}$ y F ) dependen del espacio y poseen naturaleza aleatoria, la cual está indicada por el argumento ${\textstyle \omega }$[41] . A continuación, se muestra la metodología que se debe seguir para resolver el problema con elementos finitos estocásticos.

### 2.4. Expansión de Karhunen-Loève

La expansión de Karhunen-Loève [54] de un proceso estocástico, por ejemplo ${\textstyle k({\boldsymbol {\mbox{x}}},\omega )}$ , se basa en la expansión espectral de su función de covarianza ${\textstyle R_{EE}({\boldsymbol {\mbox{x}}}_{\boldsymbol {1}},{\boldsymbol {\mbox{x}}}_{\boldsymbol {2}})}$ , donde ${\textstyle {\boldsymbol {\mbox{x}}}_{\boldsymbol {1}}}$ y ${\textstyle {\boldsymbol {\mbox{x}}}_{\boldsymbol {2}}}$ denotan las coordenadas espaciales y ${\textstyle \omega }$ denota la aleatoriedad. Por definición, la función de covarianza es simétrica y definida positiva [30]  and [54] , por tanto, las autofunciones son mutuamente ortogonales y los autovalores son reales. De esta forma, el conjunto de autofunciones y autovalores que permiten obtener la representación de ${\textstyle k({\boldsymbol {\mbox{x}}},\omega )}$ toman la forma (14):

 ${\displaystyle k({\boldsymbol {\mbox{x}}},\omega )={\overline {k}}({\boldsymbol {\mbox{x}}})+}$${\displaystyle \sum _{j=1}^{\infty }{\sqrt {{\lambda }_{j}}}k_{j}({\boldsymbol {\mbox{x}}}){\xi }_{j}(\omega )}$
( 14)

donde ${\textstyle {\overline {k}}({\boldsymbol {\mbox{x}}})}$ es el valor promedio del proceso estocástico, en este caso, del coeficiente de difusión, los términos ${\textstyle \left\{{\xi }_{j}(\omega )\right\}}$ forman un conjunto ortogonal de variables aleatorias, ${\textstyle \left\{{\lambda }_{j}\right\}}$ es el conjunto de autovalores y ${\textstyle \left\{k_{j}({\boldsymbol {\mbox{x}}})\right\}}$ es el conjunto de autofunciones de la función de covarianza, los cuales pueden ser evaluados como la solución de la siguiente integral:

 ${\displaystyle \int _{\Omega }R_{EE}({\boldsymbol {\mbox{x}}},{\boldsymbol {\mbox{y}}})k_{j}({\boldsymbol {\mbox{y}}})d{\boldsymbol {\mbox{y}}}=}$${\displaystyle {\lambda }_{j}k_{j}({\boldsymbol {\mbox{x}}})}$
( 15)

donde ${\textstyle \Omega }$ denota el dominio espacial sobre el cual se presenta el proceso ${\textstyle k({\boldsymbol {\mbox{x}}},\omega )}$ . El punto más importante en la formulación radica en que las fluctuaciones aleatorias espaciales se pueden descomponer en un conjunto de funciones deterministas en el espacio que están multiplicadas por los coeficientes aleatorios que, a su vez, son independientes de dichas funciones [46] . En este punto se deben mencionar 2 características importantes. La primera se refiere al tipo de proceso: si el proceso es gaussiano   , las variables aleatorias ${\textstyle \left\{{\xi }_{j}(\omega )\right\}}$ forman un vector gaussiano ortonormal. En segundo lugar, si el proceso tiene un alto ruido se deben incluir un gran número de términos en la sumatoria. Por el contrario, se puede incluir únicamente un término para describir el proceso. En general, la mayoría de los fenómenos físicos tienen materiales con propiedades que varían suavemente en las escalas de interés en ingeniería y, por tanto, se requerirán pocos términos en la expansión de Karhunen Loève [30] .

En el problema aquí desarrollado las variables estocásticas (1, ${\textstyle k}$ , ${\textstyle {\overset {\rightarrow }{\boldsymbol {\mbox{u}}}}}$ y F ) dadas en (13) se pueden expresar según (16), (17) y (18). En el caso de (16) se hace una expansión del coeficiente de la matriz de masa M —ver ecuación (13)— de forma similar al desarrollo hecho en Ghanem [41] :

 ${\displaystyle 1({\boldsymbol {\mbox{x}}},\omega )=1+\sum _{j_{1}=1}^{\infty }{\sqrt {{\lambda }_{j_{1}}^{1}}}1_{j_{1}}({\boldsymbol {\mbox{x}}}){\xi }_{j_{1}}(\omega )}$
( 16)

El coeficiente 1(x ,t) tiene en cuenta la variabilidad del término temporal de la ecuación. Este coeficiente es homólogo al valor de k(x ,t) que se encuentra multiplicando al laplaciano de la concentración de la variable. Un acercamiento similar se ha hecho en Ghanem y Spanos [30] , donde el parámetro asociado a la derivada temporal se divide en toda la ecuación donde se pierde de vista su significado físico que, en ecuaciones diferenciales parabólicas del tipo de transferencia de calor, tiene en cuenta el calor específico de la sustancia en estudio. Por tanto, el coeficiente 1(x ,t) permite estudiar la variación de este coeficiente, aunque su valor promedio se encuentre en el término k(x ,t).

Para el caso del coeficiente de difusión se tiene (17):

 ${\displaystyle {\boldsymbol {\mbox{k}}}({\boldsymbol {\mbox{x}}},\omega )=}$${\displaystyle {\overline {\boldsymbol {\mbox{k}}}}({\boldsymbol {x}})+\sum _{j_{k}=1}^{\infty }{\sqrt {{\lambda }_{j_{k}}^{k}}}{\boldsymbol {\mbox{k}}}_{{\boldsymbol {j}}_{\boldsymbol {k}}}({\boldsymbol {\mbox{x}}}){\xi }_{j_{k}}(\omega )}$
( 17)

y en el caso de la velocidad, se puede hacer una expansión para cada una de las componentes de la misma, esto es (18):

 ${\displaystyle {\begin{array}{l}v_{x}({\boldsymbol {\mbox{x}}},\omega )={\overline {v_{x}}}({\boldsymbol {\mbox{x}}})+\sum _{j_{x}=1}^{\infty }{\sqrt {{\lambda }_{j_{x}}^{v_{x}}}}v_{xj_{x}}({\boldsymbol {\mbox{x}}}){\xi }_{j_{x}}(\omega )\\v_{y}({\boldsymbol {\mbox{x}}},\omega )={\overline {v_{y}}}({\boldsymbol {\mbox{x}}})+\sum _{j_{y}=1}^{\infty }{\sqrt {{\lambda }_{j_{y}}^{v_{y}}}}v_{yj_{y}}({\boldsymbol {\mbox{x}}}){\xi }_{j_{y}}(\omega )\end{array}}}$
( 18)

Las funciones ${\textstyle 1_{j_{1}}({\boldsymbol {\mbox{x}}})}$ , ${\textstyle k_{j_{k}}({\boldsymbol {\mbox{x}}})}$ , ${\textstyle v_{xj_{x}}}$ y ${\textstyle v_{yj_{y}}({\boldsymbol {\mbox{x}}})}$ corresponden a las autofunciones de la función de covarianza de la variable de campo. Igualmente, los valores ${\textstyle {\lambda }_{j_{1}}^{1}}$ , ${\textstyle {\lambda }_{j_{k}}^{k}}$ , ${\textstyle {\lambda }_{j_{x}}^{v_{x}}}$ y ${\textstyle {\lambda }_{j_{y}}^{v_{y}}}$ son los autovalores. De esta forma, cada conjunto de autofunciones y autovalores conforman el denominado autopar [30] .

Reemplazando (16), (17) y (18) en (13) se tiene (19):

 ${\displaystyle \left\{\left[\int _{{\Omega }^{e}}{\left({\boldsymbol {\mbox{W}}}^{\boldsymbol {\ast }}\right)}^{T}1{\boldsymbol {\mbox{N}}}d\Omega \right]+\right.}$${\displaystyle \left.\sum _{j_{1}=1}^{\infty }\left[\int _{{\Omega }^{e}}{\left({\boldsymbol {\mbox{W}}}^{\boldsymbol {\ast }}\right)}^{T}\left({\sqrt {{\lambda }_{j_{1}}^{1}}}1_{j_{1}}({\boldsymbol {\mbox{x}}})\right){\boldsymbol {\mbox{N}}}d\Omega \right]{\xi }_{k_{1}}\right\}{\frac {\Delta {\boldsymbol {\phi }}}{\Delta t}}+}$${\displaystyle \left\{\left[\int _{{\Omega }^{e}}{\left(\nabla {\boldsymbol {\mbox{W}}}^{\boldsymbol {\ast }}\right)}^{T}{\overline {k}}({\boldsymbol {\mbox{x}}})\nabla {\boldsymbol {\mbox{N}}}d\Omega \right]+\right.}$${\displaystyle \left.\sum _{j_{k}=1}^{\infty }\left[\int _{{\Omega }^{e}}{\left(\nabla {\boldsymbol {\mbox{W}}}^{\boldsymbol {\ast }}\right)}^{T}\left({\sqrt {{\lambda }_{j_{k}}^{k}}}k_{j_{K}}({\boldsymbol {\mbox{x}}})\right)\nabla {\boldsymbol {\mbox{N}}}d\Omega \right]{\xi }_{j_{k}}\right\}{\boldsymbol {\phi }}+}$${\displaystyle \left\{\left[\int _{{\Omega }^{e}}{\left({\boldsymbol {\mbox{W}}}^{\boldsymbol {\ast }}\right)}^{T}{\overline {{\boldsymbol {\mbox{u}}}({\boldsymbol {\mbox{x}}})}}\bullet \nabla {\boldsymbol {\mbox{N}}}d\Omega \right]+\right.}$${\displaystyle \left.\sum _{j_{u}=1}^{\infty }\left[\int _{{\Omega }^{e}}{\left({\boldsymbol {\mbox{W}}}^{\boldsymbol {\ast }}\right)}^{T}\left({\sqrt {{\boldsymbol {\lambda }}_{j_{u}}^{u}}}{\boldsymbol {\mbox{u}}}_{{\boldsymbol {j}}_{\boldsymbol {u}}}{\boldsymbol {(}}{\boldsymbol {\mbox{x}}}{\boldsymbol {)}}\right)\bullet \nabla {\boldsymbol {\mbox{N}}}d\Omega \right]{\xi }_{j_{u}}\right\}{\boldsymbol {\phi }}=}$${\displaystyle {\boldsymbol {\mbox{F}}}(x,\omega )}$
( 19)

Al definir ${\textstyle {\xi }_{0}=1}$ , ${\textstyle {\lambda }_{0}^{1}=1}$ , ${\textstyle {\lambda }_{0}^{k}=1}$ , ${\textstyle {\lambda }_{0}^{v_{x}}=1}$ , ${\textstyle {\lambda }_{0}^{v_{y}}=1}$ , ${\textstyle 1_{0}({\boldsymbol {\mbox{x}}})=}$${\displaystyle 1}$ , ${\textstyle k_{0}({\boldsymbol {\mbox{x}}})=}$${\displaystyle {\overline {k}}({\boldsymbol {\mbox{x}}})}$ , ${\textstyle v_{x_{0}}({\boldsymbol {\mbox{x}}})=}$${\displaystyle {\overline {v_{x}}}({\boldsymbol {\mbox{x}}})}$ y ${\textstyle v_{y_{0}}({\boldsymbol {\mbox{x}}})=}$${\displaystyle {\overline {v_{y}}}({\boldsymbol {\mbox{x}}})}$ y definir vectorialmente la velocidad (${\textstyle {\boldsymbol {\mbox{u}}}({\boldsymbol {\mbox{x}}})=}$${\displaystyle {\left[v_{x}({\boldsymbol {\mbox{x}}}),v_{y}({\boldsymbol {\mbox{x}}})\right]}^{T}}$ ) se llega a (20):

 ${\displaystyle \left\{\sum _{j_{1}=0}^{\infty }\left[\int _{{\Omega }^{e}}{\left({\boldsymbol {\mbox{W}}}^{\boldsymbol {\ast }}\right)}^{T}\left({\sqrt {{\lambda }_{j_{1}}^{1}}}1_{j_{1}}({\boldsymbol {\mbox{x}}})\right){\boldsymbol {\mbox{N}}}d\Omega \right]{\xi }_{k_{1}}\right\}{\frac {\Delta {\boldsymbol {\phi }}}{\Delta t}}+}$${\displaystyle \left\{\sum _{j_{k}=0}^{\infty }\left[\int _{{\Omega }^{e}}{\left(\nabla {\boldsymbol {\mbox{W}}}^{\boldsymbol {\ast }}\right)}^{T}\left({\sqrt {{\lambda }_{j_{k}}^{k}}}k_{j_{K}}({\boldsymbol {\mbox{x}}})\right)\nabla {\boldsymbol {\mbox{N}}}d\Omega \right]{\xi }_{j_{k}}\right\}{\boldsymbol {\phi }}+}$${\displaystyle \left\{\sum _{j_{u}=0}^{\infty }\left[\int _{{\Omega }^{e}}{\left({\boldsymbol {\mbox{W}}}^{\boldsymbol {\ast }}\right)}^{T}\left({\sqrt {{\boldsymbol {\lambda }}_{j_{u}}^{u}}}{\boldsymbol {\mbox{u}}}_{j_{u}}{\boldsymbol {(}}{\boldsymbol {\mbox{x}}}{\boldsymbol {)}}\right)\bullet \nabla {\boldsymbol {\mbox{N}}}d\Omega \right]{\xi }_{j_{u}}\right\}{\boldsymbol {\phi }}=}$${\displaystyle {\boldsymbol {\mbox{F}}}(x,\omega )}$
( 20)

La ecuación (20) se puede escribir como (21):

 ${\displaystyle \left\{\sum _{j_{1}=0}^{\infty }{\boldsymbol {\mbox{M}}}_{{\boldsymbol {\mbox{j}}}_{\boldsymbol {1}}}{\xi }_{k_{1}}\right\}{\frac {\Delta {\boldsymbol {\phi }}}{\Delta t}}+}$${\displaystyle \left\{\sum _{j_{k}=0}^{\infty }{\boldsymbol {\mbox{K}}}_{{\boldsymbol {\mbox{j}}}_{\boldsymbol {\mbox{k}}}}{\xi }_{j_{k}}\right\}{\boldsymbol {\phi }}+}$${\displaystyle \left\{\sum _{j_{u}=0}^{\infty }{\boldsymbol {\mbox{C}}}_{{\boldsymbol {\mbox{j}}}_{\boldsymbol {\mbox{u}}}}{\xi }_{j_{u}}\right\}{\boldsymbol {\phi }}=}$${\displaystyle {\boldsymbol {\mbox{F}}}(x,\omega )}$
( 21)

donde:

 ${\displaystyle {\begin{array}{l}{\boldsymbol {\mbox{M}}}_{j_{1}}=\int _{{\Omega }^{e}}{\left({\boldsymbol {\mbox{W}}}^{\boldsymbol {\ast }}\right)}^{T}\left({\sqrt {{\lambda }_{j_{1}}^{1}}}1_{j_{1}}({\boldsymbol {\mbox{x}}})\right){\boldsymbol {\mbox{N}}}d\Omega \\{\boldsymbol {\mbox{K}}}_{j_{k}}=\int _{{\Omega }^{e}}{\left(\nabla {\boldsymbol {\mbox{W}}}^{\boldsymbol {\ast }}\right)}^{T}\left({\sqrt {{\lambda }_{j_{k}}^{k}}}k_{j_{k}}({\boldsymbol {\mbox{x}}})\right)\nabla {\boldsymbol {\mbox{N}}}d\Omega \\{\boldsymbol {\mbox{C}}}_{{\boldsymbol {\mbox{j}}}_{\boldsymbol {U}}}=\int _{{\Omega }^{e}}{\left({\boldsymbol {\mbox{W}}}^{\boldsymbol {\ast }}\right)}^{T}\left({\sqrt {{\boldsymbol {\lambda }}_{j_{u}}^{u}}}{\boldsymbol {\mbox{u}}}_{{\boldsymbol {\mbox{j}}}_{\boldsymbol {\mbox{u}}}}({\boldsymbol {\mbox{x}}})\right)\bullet \nabla {\boldsymbol {\mbox{N}}}d\Omega \end{array}}}$

Utilizando el método de Backward Euler para la integración temporal, la ecuación (21) se convierte en (22):

 ${\displaystyle \left\{\sum _{j_{1}=0}^{\infty }{\boldsymbol {\mbox{M}}}_{{\boldsymbol {\mbox{j}}}_{\boldsymbol {1}}}{\xi }_{k_{1}}\right\}{\frac {{\boldsymbol {\phi }}^{t+\Delta t}-{\boldsymbol {\phi }}^{t}}{\Delta t}}+}$${\displaystyle \left\{\sum _{j_{k}=0}^{\infty }{\boldsymbol {\mbox{K}}}_{{\boldsymbol {\mbox{j}}}_{\boldsymbol {\mbox{k}}}}{\xi }_{j_{k}}\right\}{\boldsymbol {\phi }}^{t+\Delta t}+}$${\displaystyle \left\{\sum _{j_{u}=0}^{\infty }{\boldsymbol {\mbox{C}}}_{{\boldsymbol {j}}_{\boldsymbol {u}}}{\xi }_{j_{u}}\right\}{\boldsymbol {\phi }}^{t+\Delta t}=}$${\displaystyle {\boldsymbol {\mbox{F}}}(x,\omega )}$
( 22)

Además, mediante manipulación algebraica se llega a (23).

 ${\displaystyle \left[{\frac {1}{\Delta t}}\left\{\sum _{j_{1}=0}^{\infty }{\boldsymbol {\mbox{M}}}_{{\boldsymbol {\mbox{j}}}_{\boldsymbol {1}}}{\xi }_{k_{1}}\right\}+\right.}$${\displaystyle \left.\left\{\sum _{j_{k}=0}^{\infty }{\boldsymbol {\mbox{K}}}_{{\boldsymbol {\mbox{j}}}_{\boldsymbol {\mbox{k}}}}{\xi }_{j_{k}}\right\}+\right.}$${\displaystyle \left.\left\{\sum _{j_{u}=0}^{\infty }{\boldsymbol {\mbox{C}}}_{{\boldsymbol {\mbox{j}}}_{\boldsymbol {\mbox{u}}}}{\xi }_{j_{u}}\right\}\right]{\boldsymbol {\phi }}^{t+\Delta t}=}$${\displaystyle {\boldsymbol {\mbox{F}}}(x,\omega )+{\frac {1}{\Delta t}}\left\{\sum _{j_{1}=0}^{\infty }{\boldsymbol {\mbox{M}}}_{{\boldsymbol {\mbox{j}}}_{\boldsymbol {1}}}{\xi }_{k_{1}}\right\}{\boldsymbol {\phi }}^{t}}$
( 23)

### 2.5. Expansión de los polinomios de caos

La función de covarianza de la solución, en general, no es conocida de antemano [44] . Por tanto, no se puede llevar a cabo una expansión de Karhunen Loève para describir la incógnita que hay que encontrar. En este caso la solución se puede representar por una serie de términos en una combinación no lineal que es función de las variables aleatorias ${\textstyle \left\{{\xi }_{j}(\omega )\right\}}$ . Se ha demostrado [55] que la dependencia funcional de las variables aleatorias puede ser expresada en términos de polinomios de caos, descritos como (24):

 ${\displaystyle \phi (\omega )=a_{0}{\Gamma }_{0}+\sum _{i=1}^{\infty }a_{i}{\Gamma }_{1}({\xi }_{i}(\omega ))+}$${\displaystyle \sum _{j=1}^{\infty }\sum _{i=1}^{\infty }a_{ij}{\Gamma }_{2}({\xi }_{i}(\omega ),{\xi }_{j}(\omega ))+}$${\displaystyle ...{\mbox{.}}}$
( 24)

en donde ${\textstyle {\Gamma }_{n}({\xi }_{i_{1}}(\omega ),{\xi }_{i_{2}}(\omega ),...,{\xi }_{i_{n}}(\omega ))}$ son los polinomios de caos de orden n   en las variables ${\textstyle \left\{{\xi }_{j}(\omega )\right\}}$ . Estos polinomios son la generalización de los polinomios multidimensionales de Hermite [30] . La validez de la convergencia de esta expansión es independiente del comportamiento del material. En general, cualquier variable aleatoria, con comportamiento estocástico desconocido, puede ser expresada como un polinomio en variables gaussianas según (24) [30] . Si se desarrollan los polinomios de caos, se puede obtener un conjunto de polinomios bien ordenados indicialmente descritos por ${\textstyle {\Psi }_{i}(\omega )}$[30]  and [41] , de tal forma que (24) se puede reescribir como (25):

 ${\displaystyle \varphi (\omega )=\sum _{i=1}^{\infty }{\varphi }_{i}{\Psi }_{i}(\omega )}$
( 25)

donde la variable ${\textstyle \varphi (\omega )}$ representa un valor en «un» punto nodal. Estos polinomios (${\textstyle {\Psi }_{i}(\omega )}$ ) son mutuamente ortogonales en el sentido del producto interno ${\textstyle \left\langle {\Psi }_{i}(\omega ),{\Psi }_{j}(\omega )\right\rangle }$ , el cual está dado por (26):

 ${\displaystyle {\overline {{\Psi }_{i}(\omega ){\Psi }_{j}(\omega )}}=\left\langle {\Psi }_{i}(\omega ),{\Psi }_{j}(\omega )\right\rangle }$
( 26)

donde puede apreciarse que el producto interno está definido como el promedio (esperanza estadística) del producto de los 2 polinomios [41] . Este producto interno será cero en el caso en que ${\textstyle i\not =j}$ . La caracterización completa del proceso aleatorio tiene lugar cuando se encuentran los coeficientes deterministas, dados por ${\textstyle {\varphi }_{i}}$[30] .

Por tanto, reemplazando (25) en (23) se llega a (27):

 ${\displaystyle \left[{\frac {1}{\Delta t}}\left\{\sum _{j_{1}=0}^{\infty }{\boldsymbol {\mbox{M}}}_{{\boldsymbol {\mbox{j}}}_{\boldsymbol {1}}}{\xi }_{k_{1}}\right\}+\right.}$${\displaystyle \left.\left\{\sum _{j_{k}=0}^{\infty }{\boldsymbol {\mbox{K}}}_{{\boldsymbol {\mbox{j}}}_{\boldsymbol {\mbox{k}}}}{\xi }_{j_{k}}\right\}+\right.}$${\displaystyle \left.\left\{\sum _{j_{u}=0}^{\infty }{\boldsymbol {\mbox{C}}}_{{\boldsymbol {\mbox{j}}}_{\boldsymbol {\mbox{u}}}}{\xi }_{j_{u}}\right\}\right]\left(\sum _{i=1}^{\infty }{{\boldsymbol {\phi }}_{i}}^{t+\Delta t}{\Psi }_{i}(\omega )\right)=}$${\displaystyle {\boldsymbol {\mbox{F}}}(x,\omega )+{\frac {1}{\Delta t}}\left\{\sum _{j_{1}=0}^{\infty }{\boldsymbol {\mbox{M}}}_{{\boldsymbol {\mbox{j}}}_{\boldsymbol {1}}}{\xi }_{k_{1}}\right\}\left(\sum _{i=1}^{\infty }{{\boldsymbol {\phi }}_{\boldsymbol {i}}}^{t}{\Psi }_{i}(\omega )\right)}$
( 27)

donde el vector ${\textstyle {\boldsymbol {\phi }}_{\boldsymbol {\mbox{i}}}=}$${\displaystyle [{\varphi }_{i1},{\varphi }_{i2},....,{\varphi }_{innodel}]}$ representa la i-ésima incógnita estocástica.

Con un poco de álgebra se puede definir el error r dado por (28):

 ${\displaystyle {\boldsymbol {\mbox{r}}}=\left[{\frac {1}{\Delta t}}\left\{\sum _{j_{1}=0}^{\infty }{\boldsymbol {\mbox{M}}}_{{\boldsymbol {\mbox{j}}}_{\boldsymbol {1}}}{\xi }_{k_{1}}\right\}+\right.}$${\displaystyle \left.\left\{\sum _{j_{k}=0}^{\infty }{\boldsymbol {\mbox{K}}}_{{\boldsymbol {\mbox{j}}}_{\boldsymbol {\mbox{k}}}}{\xi }_{j_{k}}\right\}+\right.}$${\displaystyle \left.\left\{\sum _{j_{u}=0}^{\infty }{\boldsymbol {\mbox{C}}}_{{\boldsymbol {\mbox{j}}}_{\boldsymbol {\mbox{u}}}}{\xi }_{j_{u}}\right\}\right]\left(\sum _{i=1}^{\infty }{{\boldsymbol {\phi }}_{i}}^{t+\Delta t}{\Psi }_{i}(\omega )\right)-}$${\displaystyle {\boldsymbol {\mbox{F}}}(x,\omega )-{\frac {1}{\Delta t}}\left\{\sum _{j_{1}=0}^{\infty }{\boldsymbol {\mbox{M}}}_{{\boldsymbol {\mbox{j}}}_{\boldsymbol {1}}}{\xi }_{k_{1}}\right\}\left(\sum _{i=1}^{\infty }{{\boldsymbol {\phi }}_{\boldsymbol {\mbox{i}}}}^{t}{\Psi }_{i}(\omega )\right)=}$${\displaystyle 0}$
( 28)

### 2.6. Rumbo a la implementación computacional

Debido a la imposibilidad de calcular la sumatoria infinita en la expansión de Karhunen-Loève y en la expansión de los polinomios de caos, se hace un truncamiento en cada serie, como puede verse en (29):

 ${\displaystyle {\boldsymbol {\mbox{r}}}=\left[{\frac {1}{\Delta t}}\left\{\sum _{j_{1}=0}^{M_{1}}{\boldsymbol {\mbox{M}}}_{{\boldsymbol {\mbox{j}}}_{\boldsymbol {1}}}{\xi }_{k_{1}}\right\}+\right.}$${\displaystyle \left.\left\{\sum _{j_{k}=0}^{M_{k}}{\boldsymbol {\mbox{K}}}_{{\boldsymbol {\mbox{j}}}_{\boldsymbol {\mbox{k}}}}{\xi }_{j_{k}}\right\}+\right.}$${\displaystyle \left.\left\{\sum _{j_{u}=0}^{M_{u}}{\boldsymbol {\mbox{C}}}_{{\boldsymbol {\mbox{j}}}_{\boldsymbol {\mbox{u}}}}{\xi }_{j_{u}}\right\}\right]\left(\sum _{i=1}^{P}{{\boldsymbol {\phi }}_{i}}^{t+\Delta t}{\Psi }_{i}(\omega )\right)-}$${\displaystyle {\boldsymbol {\mbox{F}}}(x,\omega )-{\frac {1}{\Delta t}}\left\{\sum _{j_{1}=0}^{M_{1}}{\boldsymbol {\mbox{M}}}_{{\boldsymbol {\mbox{j}}}_{\boldsymbol {1}}}{\xi }_{k_{1}}\right\}\left(\sum _{i=1}^{P}{{\boldsymbol {\phi }}_{\boldsymbol {\mbox{i}}}}^{t}{\Psi }_{i}(\omega )\right)\approx 0}$
( 29)

donde ${\textstyle M_{1}}$ , ${\textstyle M_{K}}$ y ${\textstyle M_{u}}$ determinan el truncamiento de la expansión de Karhunen-Loève para cada uno de los campos aleatorios. De igual forma, ${\textstyle P}$ determina el truncamiento de la expansión de los polinomios de caos [30] . Nótese que debido al truncamiento el error no es exactamente cero, y se convierte en un valor apreciable. Para reducir el error se hace una proyección sobre los polinomios de caos, de tal forma que el producto interno queda definido por [30] :

 ${\displaystyle \left\langle {\boldsymbol {\mbox{r}}},{\Psi }_{l}(\omega )\right\rangle =}$${\displaystyle 0}$
( 30)

Para ${\textstyle l=0,1,2,...,P}$ . Mediante manipulación algebraica se obtiene la ecuación (31):

 ${\displaystyle \left[{\frac {1}{\Delta t}}\sum _{i=1}^{P}\left\{\sum _{j_{1}=0}^{M_{1}}{\boldsymbol {\mbox{M}}}_{{\boldsymbol {\mbox{j}}}_{\boldsymbol {1}}}\left\langle {\xi }_{k_{1}}{\Psi }_{i}(\omega ){\Psi }_{l}(\omega )\right\rangle {\boldsymbol {\phi }}_{i}^{t+\Delta t}\right\}+\right.}$${\displaystyle \left.\sum _{i=1}^{P}\left\{\sum _{j_{k}=0}^{M_{k}}{\boldsymbol {\mbox{K}}}_{{\boldsymbol {\mbox{j}}}_{\boldsymbol {\mbox{k}}}}\left\langle {\xi }_{j_{k}}{\Psi }_{i}(\omega ){\Psi }_{l}(\omega )\right\rangle {\boldsymbol {\phi }}_{i}^{t+\Delta t}\right\}+\right.}$${\displaystyle \left.\sum _{i=1}^{P}\left\{\sum _{j_{u}=0}^{M_{u}}{\boldsymbol {\mbox{C}}}_{{\boldsymbol {\mbox{j}}}_{\boldsymbol {\mbox{u}}}}\left\langle {\xi }_{j_{u}}{\Psi }_{i}(\omega ){\Psi }_{l}(\omega )\right\rangle {\boldsymbol {\phi }}_{i}^{t+\Delta t}\right\}-\right.}$${\displaystyle \left.{\frac {1}{\Delta t}}\sum _{i=1}^{P}\left\{\sum _{j_{1}=0}^{M_{1}}{\boldsymbol {\mbox{M}}}_{{\boldsymbol {\mbox{j}}}_{\boldsymbol {1}}}\left\langle {\xi }_{k_{1}}{\Psi }_{i}(\omega ){\Psi }_{l}(\omega )\right\rangle {\boldsymbol {\phi }}_{i}^{t}\right\}-\right.}$${\displaystyle \left.\left\langle {\boldsymbol {\mbox{F}}}(x,\omega ){\Psi }_{l}(\omega )\right\rangle =\right.}$${\displaystyle \left.0\right]}$
( 31)

De Ghanem y Spanos [30] y de Choi et al. [56] se pueden obtener los productos internos indicados en (31). Por ejemplo, en el segundo término de (31) se obtiene (32) [57] :

 ${\displaystyle \left\langle {\xi }_{k_{1}}{\Psi }_{i}(\omega ){\Psi }_{l}(\omega )\right\rangle =}$${\displaystyle {\overline {{\xi }_{k_{1}}{\Psi }_{i}(\omega ){\Psi }_{l}(\omega )}}=}$${\displaystyle d_{j_{k}il}}$
( 32)

donde el término ${\textstyle d_{j_{k}il}}$ está dado por (33) [56]  and [57] :

 ${\displaystyle d_{j_{k}il}={\begin{array}{ll}{\frac {j_{k}!i!l!}{\left({\frac {j_{k}+i-l}{2}}\right)!\left({\frac {i+l-j_{k}}{2}}\right)!\left({\frac {j_{k}+l-i}{2}}\right)!}}&{\mbox{si}}{\mbox{ }}(j_{k}+i+l){\mbox{ }}{\mbox{es}}{\mbox{ }}{\mbox{par}},{\mbox{ }}y{\mbox{ }}l\in \left[\left|j_{k}-i\right|,j_{k}-i\right]\\0&{\mbox{ }}\end{array}}}$
( 33)

De igual manera se obtienen los demás coeficientes. Por tanto, el sistema toma la forma (34):

 ${\displaystyle \left[{\frac {1}{\Delta t}}\sum _{i=1}^{P}\left\{\sum _{j_{1}=0}^{M_{1}}{\boldsymbol {\mbox{M}}}_{{\boldsymbol {\mbox{j}}}_{\boldsymbol {1}}}d_{j_{1}il}{\boldsymbol {\phi }}_{i}^{t+\Delta t}\right\}+\right.}$${\displaystyle \left.\sum _{i=1}^{P}\left\{\sum _{j_{k}=0}^{M_{k}}{\boldsymbol {\mbox{K}}}_{{\boldsymbol {\mbox{j}}}_{\boldsymbol {\mbox{k}}}}d_{j_{k}il}{\boldsymbol {\phi }}_{i}^{t+\Delta t}\right\}+\right.}$${\displaystyle \left.\sum _{i=1}^{P}\left\{\sum _{j_{u}=0}^{M_{u}}{\boldsymbol {\mbox{C}}}_{{\boldsymbol {\mbox{j}}}_{\boldsymbol {\mbox{u}}}}d_{j_{u}il}{\boldsymbol {\phi }}_{i}^{t+\Delta t}\right\}-\right.}$${\displaystyle \left.{\frac {1}{\Delta t}}\sum _{i=1}^{P}\left\{\sum _{j_{1}=0}^{M_{1}}{\boldsymbol {\mbox{M}}}_{{\boldsymbol {\mbox{j}}}_{\boldsymbol {1}}}d_{j_{1}il}{\boldsymbol {\phi }}_{i}^{t}\right\}-\right.}$${\displaystyle \left.\left\langle {\boldsymbol {\mbox{F}}}(x,\omega ){\Psi }_{l}(\omega )\right\rangle =\right.}$${\displaystyle \left.0\right]}$
( 34)

donde ${\textstyle l=0,1,2,...,P}$ .

Expandiendo (34) y considerando que el término ${\textstyle {\boldsymbol {F}}(x,\omega )}$ es determinista (esto es: ${\textstyle \left\langle {\boldsymbol {\mbox{F}}}(x,\omega ){\Psi }_{l}(\omega )\right\rangle =}$${\displaystyle {\boldsymbol {\mbox{F}}}_{\boldsymbol {0}}}$ ), la ecuación se convierte en el sistema (35):

 ${\displaystyle {\begin{array}{c}\sum _{i=1}^{P}\left[{\frac {1}{\Delta t}}\sum _{j_{1}=0}^{M_{1}}{\boldsymbol {\mbox{M}}}_{j_{1}}d_{j_{1}i0}+\sum _{j_{k}=0}^{M_{k}}{\boldsymbol {\mbox{K}}}_{j_{k}}d_{j_{k}i0}+\sum _{j_{u}=0}^{M_{u}}{\boldsymbol {\mbox{C}}}_{j_{u}}d_{j_{u}i0}\right]{\boldsymbol {\phi }}_{i}^{t+\Delta t}={\frac {1}{\Delta t}}\sum _{i=1}^{P}\left\{\sum _{j_{1}=0}^{M_{1}}{\boldsymbol {\mbox{M}}}_{j_{1}}d_{j_{1}i0}{\boldsymbol {\phi }}_{i}^{t}\right\}-{\boldsymbol {\mbox{F}}}_{0}\\\sum _{i=1}^{P}\left[{\frac {1}{\Delta t}}\sum _{j_{1}=0}^{M_{1}}{\boldsymbol {\mbox{M}}}_{j_{1}}d_{j_{1}i1}+\sum _{j_{k}=0}^{M_{k}}{\boldsymbol {\mbox{K}}}_{j_{k}}d_{j_{k}i1}+\sum _{j_{u}=0}^{M_{u}}{\boldsymbol {\mbox{C}}}_{j_{u}}d_{j_{u}i1}\right]{\boldsymbol {\phi }}_{i}^{t+\Delta t}={\frac {1}{\Delta t}}\sum _{i=1}^{P}\left\{\sum _{j_{1}=0}^{M_{1}}{\boldsymbol {\mbox{M}}}_{j_{1}}d_{j_{1}i1}{\boldsymbol {\phi }}_{i}^{t}\right\}\\\vdots \\\sum _{i=1}^{P}\left[{\frac {1}{\Delta t}}\sum _{j_{1}=0}^{M_{1}}{\boldsymbol {\mbox{M}}}_{j_{1}}d_{jis}+\sum _{j_{k}=0}^{M_{k}}{\boldsymbol {\mbox{K}}}_{j_{k}}d_{j_{k}is}+\sum _{j_{u}=0}^{M_{u}}{\boldsymbol {\mbox{C}}}_{j_{u}}d_{j_{u}is}\right]{\boldsymbol {\phi }}_{i}^{t+\Delta t}={\frac {1}{\Delta t}}\sum _{i=1}^{P}\left\{\sum _{j_{1}=0}^{M_{1}}{\boldsymbol {\mbox{M}}}_{j_{1}}d_{j_{1}is}{\boldsymbol {\phi }}_{i}^{t}\right\}\\\vdots \\\sum _{i=1}^{P}\left[{\frac {1}{\Delta t}}\sum _{j_{1}=0}^{M_{1}}{\boldsymbol {\mbox{M}}}_{j_{1}}d_{jiP}+\sum _{j_{k}=0}^{M_{k}}{\boldsymbol {\mbox{K}}}_{j_{k}}d_{j_{k}iP}+\sum _{j_{u}=0}^{M_{u}}{\boldsymbol {\mbox{C}}}_{j_{u}}d_{j_{u}iP}\right]{\boldsymbol {\phi }}_{i}^{t+\Delta t}={\frac {1}{\Delta t}}\sum _{i=1}^{P}\left\{\sum _{j_{1}=0}^{M_{1}}{\boldsymbol {\mbox{M}}}_{j_{1}}d_{j_{1}iP}{\boldsymbol {\phi }}_{i}^{t}\right\}\end{array}}}$
( 35)

Por tanto, para la s-ésima fila se puede definir (36):

 ${\displaystyle {\begin{array}{l}A_{is}={\frac {1}{\Delta t}}\sum _{j_{1}=0}^{M_{1}}{\boldsymbol {\mbox{M}}}_{{\boldsymbol {\mbox{j}}}_{\boldsymbol {1}}}d_{jis}+\sum _{j_{k}=0}^{M_{k}}{\boldsymbol {\mbox{K}}}_{{\boldsymbol {\mbox{j}}}_{\boldsymbol {\mbox{k}}}}d_{j_{k}is}+\sum _{j_{u}=0}^{M_{u}}{\boldsymbol {\mbox{C}}}_{{\boldsymbol {\mbox{j}}}_{\boldsymbol {\mbox{u}}}}d_{j_{u}is}\\B_{is}={\frac {1}{\Delta t}}\sum _{j_{1}=0}^{M_{1}}{\boldsymbol {\mbox{M}}}_{{\boldsymbol {\mbox{j}}}_{\boldsymbol {1}}}d_{j_{1}is}\end{array}}}$
( 36)

De tal manera que (35) se puede escribir de la forma (37):

 ${\displaystyle {\begin{array}{c}A_{00}{\boldsymbol {\phi }}_{0}^{n+1}+A_{10}{\boldsymbol {\phi }}_{1}^{n+1}+A_{20}{\boldsymbol {\phi }}_{2}^{n+1}+\cdots +A_{P0}{\boldsymbol {\phi }}_{P}^{n+1}=B_{00}{\boldsymbol {\phi }}_{0}^{n+1}+B_{10}{\boldsymbol {\phi }}_{1}^{n}+B_{20}{\boldsymbol {\phi }}_{2}^{n}+\cdots +B_{P0}{\boldsymbol {\phi }}_{P}^{n}-{\boldsymbol {\mbox{F}}}_{\boldsymbol {0}}\\A_{01}{\boldsymbol {\phi }}_{0}^{n+1}+A_{11}{\boldsymbol {\phi }}_{1}^{n+1}+A_{21}{\boldsymbol {\phi }}_{2}^{n+1}+\cdots +A_{P1}{\boldsymbol {\phi }}_{P}^{n+1}=B_{01}{\boldsymbol {\phi }}_{0}^{n+1}+B_{11}{\boldsymbol {\phi }}_{1}^{n}+B_{21}{\boldsymbol {\phi }}_{2}^{n}+\cdots +B_{P1}{\boldsymbol {\phi }}_{P}^{n}\\\vdots \\A_{0s}{\boldsymbol {\phi }}_{0}^{n+1}+A_{1s}{\boldsymbol {\phi }}_{1}^{n+1}+A_{2s}{\boldsymbol {\phi }}_{2}^{n+1}+\cdots +A_{Ps}{\boldsymbol {\phi }}_{P}^{n+1}=B_{0s}{\boldsymbol {\phi }}_{0}^{n+1}+B_{1s}{\boldsymbol {\phi }}_{1}^{n}+B_{2s}{\boldsymbol {\phi }}_{2}^{n}+\cdots +B_{Ps}{\boldsymbol {\phi }}_{P}^{n}\\\vdots \\A_{0P}{\boldsymbol {\phi }}_{0}^{n+1}+A_{1P}{\boldsymbol {\phi }}_{1}^{n+1}+A_{2P}{\boldsymbol {\phi }}_{2}^{n+1}+\cdots +A_{PP}{\boldsymbol {\phi }}_{P}^{n+1}=B_{0P}{\boldsymbol {\phi }}_{0}^{n+1}+B_{1P}{\boldsymbol {\phi }}_{1}^{n}+B_{2P}{\boldsymbol {\phi }}_{2}^{n}+\cdots +B_{PP}{\boldsymbol {\phi }}_{P}^{n}\end{array}}}$
( 37)

Y en términos matriciales (38):

 ${\displaystyle {\underset {\underline {\underline {\boldsymbol {\mbox{A}}}}}{\underbrace {\left({\begin{array}{ccc}A_{00}&\cdots &A_{0P}\\\vdots &\ddots &\vdots \\A_{0P}&\cdots &A_{PP}\end{array}}\right)} }}{\underset {\underline {\underline {\boldsymbol {\Phi }}}}{\underbrace {\left({\begin{array}{c}{\boldsymbol {\phi }}_{0}\\\vdots \\{\boldsymbol {\phi }}_{P}\end{array}}\right)} }}^{n+1}=\left({\begin{array}{ccc}B_{00}&\cdots &B_{0P}\\\vdots &\ddots &\vdots \\B_{0P}&\cdots &B_{PP}\end{array}}\right){\left({\begin{array}{c}{\boldsymbol {\phi }}_{0}\\\vdots \\{\boldsymbol {\phi }}_{P}\end{array}}\right)}^{n}={\underset {\underline {\underline {\boldsymbol {\mbox{R}}}}}{\underbrace {\left({\begin{array}{c}{\boldsymbol {\mbox{R}}}_{0}\\\vdots \\{\boldsymbol {\mbox{R}}}_{P}\end{array}}\right)} }}}$
( 38)

### 2.7. Sobre los autovalores y las autofunciones

De forma general, en 2 dimensiones, la ecuación (15) se puede escribir como:

 ${\displaystyle \int _{\Omega }R_{EE}(x_{1},y_{1};x_{2},y_{2})k_{j}(x_{2},y_{2})dx_{2}dy_{2}=}$${\displaystyle {\lambda }_{j}k_{j}(x_{1},y_{1})}$
( 39)

donde ${\textstyle R_{EE}({\boldsymbol {\mbox{x}}},{\boldsymbol {\mbox{y}}})}$ es la función de covarianza. En este trabajo se supone (40) [57] :

 ${\displaystyle R_{EE}(x_{1},y_{1};x_{2},y_{2})=C_{0}exp\left(-{\frac {\left|x_{1}-x_{2}\right|}{b_{1}}}-\right.}$${\displaystyle \left.{\frac {\left|y_{1}-y_{2}\right|}{b_{2}}}\right)}$
( 40)

donde ${\textstyle C_{0}}$ es el valor de la varianza [52] y ${\textstyle b_{1}}$ y ${\textstyle b_{2}}$ son longitudes de correlación del campo aleatorio [30]  and [52] . La ecuación (40) es válida en un dominio cuadrado ${\textstyle \Omega =\left(-a,a\right)\times \left(-\right.}$${\displaystyle \left.a,a\right)}$ . Siguiendo el procedimiento dado en [30] se puede llevar a cabo una separación de variables, de tal manera que la integral (39) se puede calcular separando términos según cada coordenada, por lo que se obtienen los autovalores y las autofunciones que se muestran en (41):

 ${\displaystyle {\begin{array}{l}k_{j}(x_{1},y_{1})=k_{j_{1}}(x_{1})\ast k_{j_{2}}(x_{1})\\{\lambda }_{j}={\lambda }_{j_{1}}\ast {\lambda }_{j_{2}}\end{array}}}$
( 41)

Por tanto, reemplazando (40) en (39) bajo la condición (41) se obtiene (42):

 ${\displaystyle \left(\int _{x}{\sqrt {C_{0}}}exp\left(-{\frac {\left|x_{1}-x_{2}\right|}{b_{1}}}\right)k_{j_{1}}(x_{2})dx_{2}\right)\left(\int _{y}{\sqrt {C_{0}}}exp\left(-\right.\right.}$${\displaystyle \left.\left.{\frac {\left|y_{1}-y_{2}\right|}{b_{2}}}\right)k_{j_{2}}(y_{2})dy_{2}\right)=}$${\displaystyle \left({\lambda }_{j_{1}}k_{j_{1}}(x_{1})\right)\left({\lambda }_{j_{2}}k_{j_{2}}(y_{2})\right)}$
( 42)

donde, de acuerdo con Ghanem y Spanos [30] (42) se puede resolver de forma independiente, así:

 ${\displaystyle {\begin{array}{l}\int _{x}{\sqrt {C_{0}}}exp\left(-{\frac {\left|x_{1}-x_{2}\right|}{b_{1}}}\right)k_{j_{1}}(x_{2})dx_{2}={\lambda }_{j_{1}}k_{j_{1}}(x_{1})\\\int _{y}{\sqrt {C_{0}}}exp\left(-{\frac {\left|y_{1}-y_{2}\right|}{b_{2}}}\right)k_{j_{2}}(y_{2})dy_{2}={\lambda }_{j_{2}}k_{j_{2}}(y_{2})\end{array}}}$
( 43)

donde ${\textstyle -a\leq x_{2}\leq a}$ y ${\textstyle -a\leq y_{2}\leq a}$ . El valor absoluto se puede remover mediante su definición, como se muestra a continuación —se lleva a cabo la demostración sobre la primera ecuación de (43)—:

 ${\displaystyle {\lambda }_{j_{1}}k_{j_{1}}(x_{1})={\sqrt {C_{0}}}{\int }_{-a}^{x_{1}}exp\left(-\right.}$${\displaystyle \left.{\frac {\left(x_{1}-x_{2}\right)}{b_{1}}}\right)k_{j_{1}}(x_{2})dx_{2}+}$${\displaystyle {\sqrt {C_{0}}}{\int }_{x_{1}}^{a}exp\left({\frac {\left(x_{1}-x_{2}\right)}{b_{1}}}\right)k_{j_{1}}(x_{2})dx_{2}}$
( 44)

Diferenciando con respecto a ${\textstyle x_{1}}$ , se tiene (45):

 ${\displaystyle {\lambda }_{j_{1}}{\frac {dk_{j_{1}}(x_{1})}{dx_{1}}}=-{\sqrt {C_{0}}}{\frac {1}{b_{1}}}exp\left(-\right.}$${\displaystyle \left.{\frac {\left(x_{1}\right)}{b_{1}}}\right){\int }_{-a}^{x_{1}}exp\left({\frac {x_{2}}{b_{1}}}\right)k_{j_{1}}(x_{2})dx_{2}+}$${\displaystyle {\sqrt {C_{0}}}{\frac {1}{b_{1}}}exp\left({\frac {\left(x_{1}\right)}{b_{1}}}\right){\int }_{x_{1}}^{a}exp\left(-\right.}$${\displaystyle \left.{\frac {\left(x_{2}\right)}{b_{1}}}\right)k_{j_{1}}(x_{2})dx_{2}}$
( 45)

Diferenciando, nuevamente, con respecto a ${\textstyle x_{1}}$ , se tiene (46):

 ${\displaystyle {\lambda }_{j_{1}}{\frac {d^{2}k_{j_{1}}(x_{1})}{d{x_{1}}^{2}}}=\left(h^{2}{\lambda }_{j_{1}}-\right.}$${\displaystyle \left.2{\sqrt {C_{0}}}h\right)k_{j_{1}}(x_{1})}$
( 46)

Definiendo ${\textstyle {\alpha }^{2}=\left(2{\sqrt {C_{0}}}h-\right.}$${\displaystyle \left.h^{2}{\lambda }_{j_{1}}\right)/{\lambda }_{j_{1}}}$ , la ecuación (46) se convierte en (47):

 ${\displaystyle {\frac {d^{2}k_{j_{1}}(x_{1})}{d{x_{1}}^{2}}}+{\alpha }^{2}k_{j_{1}}(x_{1})=}$${\displaystyle 0}$
( 47)

donde ${\textstyle -a\leq x_{1}\leq a}$ . La solución de (47) es de la forma (48):

 ${\displaystyle k_{j}(x_{1})=c_{1}e^{i\alpha x_{1}}+c_{2}e^{-i\alpha x_{1}},{\alpha }^{2}\geq 0}$
( 48)

donde ${\textstyle c_{1}}$ y ${\textstyle c_{2}}$ son constantes. Aplicando las condiciones de contorno se obtiene la solución dada por (49) [56] :

 ${\displaystyle {\lambda }_{j_{1}}={\frac {2{\sqrt {C_{0}}}h}{h^{2}+{{\alpha }_{j_{1}}}^{2}}}}$
( 49)

Y se obtienen para las autofunciones (50) [56] :

 ${\displaystyle k_{j_{1}}(x_{1})={\begin{array}{ll}{\frac {cos\left({\frac {\left(j_{1}-1\right)\pi }{2a}}x_{1}\right)}{\sqrt {a+\left({\frac {sin(2{\alpha }_{j_{1}}a)}{2{\alpha }_{j_{1}}}}\right)}}}&j_{1}{\mbox{ par}}\\{\frac {sin\left({\frac {\left(j_{1}-1\right)\pi }{2a}}x_{1}\right)}{\sqrt {a-\left({\frac {sin(2{\alpha }_{j_{1}}a)}{2{\alpha }_{j_{1}}}}\right)}}}&j_{1}{\mbox{ impar}}\end{array}}}$
( 50)

donde se tiene el arreglo ${\textstyle {\alpha }_{1}<{\alpha }_{2}<...<{\alpha }_{n}}$ , que permite encontrar los autovalores ordenados mediante ${\textstyle {\lambda }_{1}>{\lambda }_{2}>...>{\lambda }_{n}}$ .

### 2.8. Aspectos computacionales y condiciones de contorno e iniciales

Con el fin de evaluar la precisión del método de Petrov-Galerkin en conjunto con los elementos finitos estocásticos se analizaron 2 problemas, con advección dominante. Los problemas fueron implementados en lenguaje FORTRAN y solucionados para diferentes mallas, tal como se muestra a continuación.

Para los ejemplos que se han desarrollado en este artículo se utilizaron dominios cuadrados de lado unitario (figura 2 ). Además, se han utilizado mallas de 2.601 (2.500), 676 (625) y 169 (144) nodos (elementos). Los elementos utilizados son cuadriláteros con funciones de forma bilineales.

 Figura 2. Malla utilizada en el análisis por elementos finitos estocásticos espectrales.

Las condiciones iniciales son impuestas a la variable determinista según el problema. Para las variables estocásticas las condiciones iniciales son nulas sin importar las condiciones iniciales de la variable φ0[30] . En los ejemplos aquí desarrollados se utilizan unas condiciones iniciales según el tipo de ejemplo (φ0 ). De igual forma, las condiciones de Neumann (naturales) y Dirichlet son impuestas en los nodos de la frontera sobre la ecuación determinista y, para las variables estocásticas, se colocan condiciones Neumann y Dirichlet en los mismos nodos, con valor nulo, respectivamente. De otro lado, en la función de ponderación W* se utiliza la velocidad determinista para el método SUPG.

## 3. Experimentación numérica

### 3.1. Caso 1: movimiento convectivo de una función gaussiana

Este ejemplo, desarrollado entre otros autores por Wang et al. [28] , está definido por la ecuación diferencial (51), para un dominio ${\textstyle \Omega =\left(-0,5,0,5\right)\times \left(-\right.}$${\displaystyle \left.0,5,0,5\right)}$ .

 ${\displaystyle {\frac {\partial \varphi }{\partial t}}-k\Delta \varphi +}$${\displaystyle {\boldsymbol {\mbox{u}}}\cdot \nabla \varphi =0}$
( 51)

con ${\textstyle k=0,0001}$ y un campo de velocidad rotacional ${\textstyle {\boldsymbol {u}}=\left[-4y,4x\right]}$ . Las condiciones de borde definidas para las 4 fronteras del problema son condiciones de Dirichlet homogéneas, en tanto que las condiciones iniciales están definidas por la expresión (53), con ${\textstyle \sigma =0,0477}$ . La solución analítica del problema se plantea en la expresión (52).

 ${\displaystyle \varphi (x,y,t)={\frac {2{\sigma }^{2}}{2{\sigma }^{2}+4kt}}e^{-\lambda }}$
( 52)

Por tanto, las condiciones iniciales están dadas por (53):

 ${\displaystyle \varphi (x,y,0)=e^{-{\frac {{\left(x+0,25\right)}^{2}+y^{2}}{2{\sigma }^{2}}}}}$
( 53)

en esta última expresión:

 ${\displaystyle {\begin{array}{l}\lambda ={\frac {{\left({\overline {x}}+0,25\right)}^{2}+{\overline {y}}^{2}}{2{\sigma }^{2}+4kt}}\\{\overline {x}}=xCos(4t)+ySin(4t)\\{\overline {y}}=-xSin(4t)+yCos(4t)\end{array}}}$

Este problema se caracteriza por tener dentro del dominio espacial tanto zonas de convección dominante, ubicadas cerca de los bordes del problema, como zonas de difusión dominante, ubicadas en la región cercana al centro del dominio. Estas zonas, definidas espacialmente, no cambian en el tiempo, dado que el coeficiente de difusión y el campo de velocidades no son función de esta variable. Por lo anterior, este caso resulta especialmente interesante para evaluar la capacidad de las técnicas numéricas para aproximar problemas con condiciones convectivas variables sobre todo el dominio.

Los resultados deterministas que se obtuvieron en este ejemplo se muestran en la figura 3 (a-e). Se empleó una discretización temporal con pasos de tiempo ${\textstyle \Delta t=0,0005}$ y un tiempo final ${\textstyle t_{f}=0,5}$ . Se muestran los resultados para la malla de 2.500 (figs. 3 a y 3d), 625 (figs. 3 b y 3e) y 144 elementos (figs. 3 c y 3f), respectivamente. Además, la figura 3 muestra la solución Petrov-Galerkin (figs. 3 a-3c) y Bunbov-Galerkin (figs. 3 d-3f). Cada columna representa los siguientes instantes de tiempo (de izquierda a derecha): 0,005, 0,125, 0,250, 0,375 y 0,500 unidades de tiempo (adimensional).

 Figura 3. Resultados para el problema determinista con el uso de Petrov-Galerkin para a. 2.500 elementos; b. 625 elementos; y c. 144 elementos, y cuando se utiliza el método Bubnov Galerkin para d. 2.500 elementos; e. 625 elementos; y f. 144 elementos.

En la figura 3 se puede observar que las soluciones obtenidas por los 2 métodos son similares, cualitativamente, durante el desplazamiento de la función gaussiana . Sin embargo, se presentan inestabilidades que crecen de manera abrupta al aumentar el tamaño de la malla. Cuando se utiliza Bubnov-Galerkin ( figs. 3 d-3f), se presentan inestabilidades que se representan por líneas de contorno. Además, las oscilaciones alcanzan valores comparables con la altura de la función gaussiana . También se presenta una dispersión considerable de la variable, tanto que la concentración de la función gaussiana baja, incluso a la mitad de su valor, comparado con el método SUPG ( figs. 3 a-3c).

Para el caso estocástico se utilizaron los siguientes datos: ${\textstyle C_{0}=1.e-4}$ , ${\textstyle \left|a\right|=0,5}$ , ${\textstyle b_{1}=b_{2}=0,25}$ (longitud de correlación en x e y , respectivamente). El número de términos en la expansión de Karhunen-Loève es M  =  5 (en todos los parámetros), y en la expansión del polinomio de caos P  =  5   . Con los valores mencionados se encuentra que ${\textstyle {\alpha }_{1}=2,1537,{\alpha }_{2}=}$${\displaystyle 4,5779,{\alpha }_{3}=7,2868,{\alpha }_{4}=10,1732{\mbox{.}}}$ En la figura 4 se muestran los resultados para un solo instante de tiempo (t = 0,500): el tiempo final de la simulación. Ordenado mediante filas se encuentran cada una de las variables estocásticas ϕ0 , ϕ1 , ϕ2ϕ4 , respectivamente. Además, se ordena por columnas las diferentes mallas y los 2 tipos de solución: i a y i b para dominios con 144 elementos, ii a y ii b para 625 elementos y iii a y iii b para 2.500 elementos. Las columnas a muestran los resultados Bubnov Galerkin y las columnas b los resultados Petrov Galerkin.

 Figura 4. Patrones para la solución del problema estocástico de la ecuación (51) con los parámetros del primer ejemplo. Filas A, B, C…E corresponden a las variables estocásticas ${\textstyle {\varphi }_{0}}$ , ${\textstyle {\varphi }_{1}}$ , ${\textstyle {\varphi }_{2}}$ … ${\textstyle {\varphi }_{4}}$ , respectivamente. Las columnas i a, ii a, iii a corresponden a la solución con Bubnov-Galerkin y las columnas i b, ii b y iii b para Petrov-Galerkin. Los patrones corresponden al tiempo final de la simulación t = 0,500 (adimensional).

En este artículo se ha elegido un número de términos de la expansión M = 5, de forma similar al trabajo desarrollado por Ghanem y Spanos [30] . Según se establece en Ghanem y Spanos [30] y en Choi et al. [56] los primeros 5 términos contienen la mayor cantidad de información de la variación de la solución del problema. Adicionalmente, el coste computacional de aumentar el número de términos no justifica la calidad de la información obtenida del campo estocástico [30] .

Para el caso de las columnas Bubnov Galerkin (columnas del tipo a   ) las filas (${\textstyle {\varphi }_{0}}$ ) muestran que a medida que disminuye el tamaño del elemento de 144, 625 hasta 2.500 elementos (i, ii y iii , respectivamente) la respuesta se vuelve menos oscilante y dispersiva, como se comentó en el caso determinista. Para el caso de 144 elementos (columna i a) se muestran picos de concentración de la variable ${\textstyle {\varphi }_{0}}$ cerca de la frontera inferior izquierda y superior derecha. En la columna ii a se muestra que los picos de concentración y las inestabilidades disminuyen, pero se mantiene un gran número de sitios de inestabilidad (mostrado como líneas de contorno). En la columna iii a se muestra que los resultados son similares a los hallados mediante SUPG. Las figuras iii a y iii b son comparables debido a que el número de Peclet [14] es bajo y cercano a la unidad, por lo que el efecto de la función de ponderación perturbadora —ecuación (10)— es casi despreciable.

La siguiente fila de la figura 4 , para ${\textstyle {\varphi }_{1}}$ , muestra el primer coeficiente de la expansión del polinomio de caos [30] . Nuevamente las columas i a, ii a y iii a muestran los resultados para la solución Bubnov Galerkin. En dichas figuras se muestra una alta dispersión e inestabilidad de la variable estocástica, especialmente para el caso i a, que va disminuyendo para ii a y iii a. En la columna i a se observa que la variación del componente estadístico es alta y casi uniforme en todo el dominio, lo que indica que existe una alta incertidumbre en la variable debido a la solución oscilante que introduce el método Bubnov Galerkin. Una situación similar se observa en la columna ii a, donde se nota una alta variación en el dominio, aunque en la posición final de la campana gaussiana se presenta un pico negativo. En la figura iii a se observa un resultado similar al de la figura ii a. Se muestra un valor altamente localizado en el máximo de la campana gaussiana . Alrededor de la campana (en las columnas ii a y iii a) se presenta un gradiente de incertidumbre que se debe a la oscilación que se introduce por utilizar el método de Bubnov. En las figuras de la columna b (i b, ii b y iii b) para SUPG se muestran resultados similares a los de las columnas a . Sin embargo, en i b se observa la ubicación de la incertidumbre en el pico de la campana gaussiana , la cual aparece tenuemente. Esto indica que el método de Petrov-Galerkin logra disminuir la inestabilidad en la solución que causa variabilidad estadística. En las figuras ii b y iii b se muestran patrones similares a los del método de Bubnov Galerkin.

Las filas ${\textstyle {\varphi }_{2}}$ y ${\textstyle {\varphi }_{3}}$ muestran resultados similares que para ${\textstyle {\varphi }_{1}}$ . Nuevamente, en el caso de Bubnov-Galerkin (columnas i a, ii a y iii a) se nota la alta variabilidad en la respuesta del sistema. Además, en i a se muestra que el método Bubnov no es capaz de captar la zona ubicada en el centro de la campana gaussiana. Para los patrones hallados mediante el método SUPG se ubica un mínimo en el centro de la campana de Gauss. Alrededor de este mínimo se hallan altos valores de concentración de la variable que se deben al gradiente que se genera en la interfaz de la campana y el dominio restante, y que causa una incertidumbre en la solución.

Por último, en la fila ${\textstyle {\varphi }_{4}}$ se observan los patrones obtenidos para la cuarta componente del polinomio de caos. En las columnas i a y ii a se muestran puntos de alta concentración en la región inferior izquierda, donde se presenta la mayor convección (ver ecuación ${\textstyle {\overset {\rightarrow }{u}}=\left[-\right.}$${\displaystyle \left.4y,4x\right]}$ ). En i b, ii b, iii a y iii b se reconoce un valor alto en el centro de la campana gaussiana . Nuevamente, en las figuras iii a y iii b, donde la malla es más fina que en las columnas restantes, se observa una marcación importante del máximo de incertidumbre en el centro de la campana gaussiana .

### 3.2. Caso 2: movimiento convectivo de un frente de onda

Este ejemplo, desarrollado en Berzins [29] , modela la propagación de un frente de onda dentro de un dominio cuadrado ${\textstyle \Omega =\left(0,1\right)\times \left(0,1\right)}$ . El problema está definido por medio de la ecuación (51), para una constante ${\textstyle k=0,01}$ y un campo de velocidad definido por el vector ${\textstyle {\boldsymbol {u}}=\left[\omega (x,t),\omega (y,t)\right]}$ , en donde:

 ${\displaystyle \omega (r,t)={\frac {0,1A+0,5B+C}{A+B+C}},{\mbox{ }}{\mbox{con}}{\mbox{ }}r=}$${\displaystyle x,y}$
( 54)
 ${\displaystyle A(n,t)=e^{\frac {-0,05\left(n-0,5+4,95(t+0,3)\right)}{k}}}$
( 55)
 ${\displaystyle B(n,t)=e^{\frac {-0,25\left(n-0,5+0,75(t+0,3)\right)}{k}}}$
( 56)
 ${\displaystyle C(n,t)=e^{\frac {-0,5\left(n-0,375\right)}{k}}}$
( 57)

Las condiciones de borde definidas para este ejemplo son condiciones de Neumman homogéneas, en tanto que las condiciones iniciales están dadas por la ecuación (58).

 ${\displaystyle \varphi (x,y,0)=\omega (x,0).\omega (y,0)}$
( 58)

El ejemplo resulta interesante porque el campo de velocidades asociado al término convectivo es variable en el tiempo y en el espacio, formando una zona predominantemente convectiva que crece, a la vez que se reduce la zona de difusión dominante, tal como se muestra en la figura 5 . Esta condición exige una mayor capacidad de estabilización de la técnica numérica empleada, tanto en la dimensión espacial como en la temporal.

 Figura 5. Distribución de las zonas de alta advección y difusión para el caso 2.

Los resultados deterministas que se obtuvieron en este ejemplo se muestran en la figura 6 , donde se empleó una discretización temporal con pasos de tiempo ${\textstyle \Delta t=0,0005}$ y un tiempo final ${\textstyle t_{f}=0,5}$ .

 Figura 6. Resultados para el problema determinista del segundo caso con el uso de Petrov-Galerkin para a. 2.500 elementos; b. 625 elementos; y c. 144 elementos, y cuando se utiliza el método Bubnov Galerkin para d. 2.500 elementos; e. 625 elementos; y f. 144 elementos.

Se observa (fig. 6 ) que la malla más fina, formada por 2.500 elementos, no tiene diferencias significativas en las aproximaciones alcanzadas por las 2 técnicas. En especial se debe observar la solución con Bubnov Galerkin para 144 elementos, donde se nota una oscilación cercana al límite de la zona con convección fuerte (ver la fig. 5 para la referencia). Ambas técnicas, para 625 y 2.500 elementos, logran aproximaciones estables y exactas, aun en la zona de convección dominante.

Para el caso estocástico, al igual que en el anterior ejemplo, se utilizaron los siguientes datos: ${\textstyle C_{0}=1.e-4}$ , ${\textstyle \left|a\right|=0.5}$ , ${\textstyle b_{1}=b_{2}=0,25}$ (longitud de correlación en x e y , respectivamente). El número de términos en la expansión de Karhunen-Loève es M  =  5 , y en la expansión del polinomio de caos P  =  5   . Con los valores mencionados se encuentra que ${\textstyle {\alpha }_{1}=2,1537,{\alpha }_{2}=}$${\displaystyle 4,5779,{\alpha }_{3}=7,2868,{\alpha }_{4}=10,1732}$ . En la figura 7 se muestran los resultados para un solo instante de tiempo (t = 0,500): el tiempo final de la simulación. Los datos son ordenados de la misma forma que para el anterior caso, esto es: mediante filas se encuentran cada una de las variables estocásticas ${\textstyle {\varphi }_{0}}$ , ${\textstyle {\varphi }_{1}}$ , ${\textstyle {\varphi }_{2}}$${\textstyle {\varphi }_{4}}$ , respectivamente. Además, se ordenan por columnas las diferentes mallas y los 2 tipos de solución: i a y i b para dominios con 144 elementos, ii a y ii b para 625 elementos y iii a y iii b para 2.500 elementos. Las columnas a muestran los resultados Bubnov Galerkin y las columnas b los resultados Petrov Galerkin.

 Figura 7. Patrones para la solución del problema estocástico de la ecuación (51) con los parámetros del segundo ejemplo. Filas a, b, c…e corresponde a las variables estocásticas ${\textstyle {\varphi }_{0}}$ , ${\textstyle {\varphi }_{1}}$ , ${\textstyle {\varphi }_{2}}$ … ${\textstyle {\varphi }_{4}}$ , respectivamente. Las columnas i a, ii a y iii a corresponden a la solución con Bubnov-Galerkin y las columnas i b, ii b y iii b para Petrov-Galerkin. Los patrones corresponden al tiempo final de la simulación t = 0,500 (adimensional).

La fila a de la figura 7 muestra los resultados para ${\textstyle {\varphi }_{0}}$ . En la primera columna (i a) se observa el frente de onda con una oscilación apreciable en la zona de alta convección. Por detrás del frente de onda se notan 2 oscilaciones que producen máximos localizados que tienen el mismo perfil que el límite de la onda. En las figuras i b, ii a, ii b, iii a y iii b se observa un patrón continuo y suave en el dominio. La zona de alta convección no presenta oscilaciones.

La fila b   presenta el primer coeficiente del polinomio de caos: ${\textstyle {\varphi }_{1}}$ . La columna i a (cuya solución se obtuvo mediante Bubnov-Galerkin) presenta zonas de alta incertidumbre en la zona convectiva que se manifiesta con 2 oscilaciones. Para la siguiente columna i b se ha utilizado el método de Petrov Galerkin, lo que permite hallar patrones más suaves y sitios de menor incertidumbre. Las figuras i b, ii a, ii b, iii a y iii b muestran que cuando el refinamiento de la malla es mayor, y con el uso del método de estabilización SUPG, la zona de mayor gradiente de la variable estocástica (zonas de incertidumbre) se ubica en la interfaz del frente de onda. Obsérvese, además, que en la punta de la onda se presentan oscilaciones que dan cuenta de la alta variabilidad de la respuesta.

En la fila c   , para ${\textstyle {\varphi }_{2}}$ , columna i a, se observa, nuevamente, una gran dispersión en la solución. El patrón muestra una zona extendida de alta variabilidad. La columna i b muestra la formación de altas concentraciones de ${\textstyle {\varphi }_{2}}$ cerca y alrededor de la interfaz de la onda. Si se compara con el caso para ${\textstyle {\varphi }_{1}}$ se observa una oscilación en dirección longitudinal de la interfaz. A medida que se disminuye el tamaño del elemento y se utiliza el método de Petrov-Galerkin, el patrón de respuesta de ${\textstyle {\varphi }_{2}}$ se vuelve cada vez más concentrado en la interfaz de la onda. Adicionalmente se observa un patrón con altas concentraciones en la zona de convección (fig. 5 ).

Las filas d y e de la figura 7 muestran resultados para ${\textstyle {\varphi }_{3}}$ y ${\textstyle {\varphi }_{4}}$ , respectivamente. Una vez más se observa una gran variabilidad en la interfaz de la onda. Adicionalmente, se presentan 2 y 3 oscilaciones longitudinales en d y e , respectivamente. El patrón es de nuevo menos disperso en el caso de mallas con menor tamaño de elemento y con el uso de Petrov-Galerkin.

## 4. Discusión y conclusiones

Se analizaron 2 problemas de difusión-convección empleando 2 técnicas de solución numérica: SUPG y el método de Bubnov-Galerkin. En cada caso, adicionalmente, se tuvo en cuenta que los parámetros de la ecuación tienen un comportamiento probabilista, por lo que se utiliza el método de elementos finitos estocásticos espectrales. Se confirmó que las aproximaciones que se alcanzan con el método SUPG logran soluciones libres de oscilaciones falsas. Además, el método SUPG exhibe una difusión numérica elevada que crece rápidamente con el aumento en el tamaño de la malla.

De otro lado, las características de los patrones de la solución numérica son heredadas por el método de elementos finitos estocásticos. En dominios con discretización de malla gruesa se tienen patrones altamente dispersivos que se trasladan a cada uno de los coeficientes de los polinomios de caos. Igualmente, con el método de Bubnov-Galerkin se obtienen coeficientes estocásticos (${\textstyle {\varphi }_{i}}$ ) que exhiben una alta dispersión y variabilidad de la respuesta en todo el dominio. Por el contrario, conforme se disminuye el tamaño del elemento y se usa el método de Petrov-Galerkin, se hallan patrones con mayor concentración cerca de las zonas donde se presentan los mayores gradientes de difusión y convección. En general, aunque el método de Petrov-Galerkin y las mallas refinadas eliminen las oscilaciones falsas que se deben al término convectivo, los coeficientes estocásticos del polinomio de caos pueden captar aquellas zonas de oscilación leve.

La incertidumbre de las propiedades del material y de la velocidad convectiva puede ser incluida con funciones de covarianza que, en general, deben ser diferentes. Además, se debe tener en cuenta que la variabilidad de las propiedades puede tener diferentes escalas de incertidumbre. Es así como se puede incluir una mayor cantidad de términos en aquellos parámetros donde se requiere captar con mayor precisión la variabilidad en el espacio. En este artículo se utilizó la misma función de covarianza, funciones de correlación y valores de variación para todos los parámetros estocásticos. La utilización de un mayor número de términos en la expansión de Karhunen Loève indicaría que se presenta incertidumbre a escala micro y macroscópica. Por el contrario, una pequeña cantidad de términos indicaría que la incertidumbre se presenta a escala microscópica [30] . En estudios futuros se llevarán a cabo simulaciones que permitan modelar fenómenos descritos por ecuaciones de difusión convección, con parámetros que se encuentran en diferentes escalas.

En este trabajo se hicieron pruebas previas (no reportadas) donde se aumenta el valor del coeficiente de variación C0 y se aumenta la cantidad de términos en la expansión de Karhunen Loève y los polinomios de caos. Al incluir una gran cantidad de términos o utilizar un coeficiente de variación con alto valor, diverge la solución del problema. Este análisis permite, desde el punto de vista numérico, encontrar el límite en donde los parámetros tienen su valor máximo de variación y longitud de correlación que dan cuenta de la variabilidad. Unos valores altos de incertidumbre pueden dirigir la solución hacia patrones «caóticos», desde el punto de vista numérico.

En este trabajo se utilizó la función de ponderación perturbada de SUPG para estabilizar la solución del problema de difusión convección. Para calcular la parte de perturbación —ver ecuaciones (5) y (10)—, se propone el uso del coeficiente cero (determinista) de la velocidad del medio. De nuevo, si no se utilizan variables de campo y coeficientes estocásticos, la formulación se reduce al problema determinista.

En conclusión, el método de los elementos finitos estocásticos espectrales con aplicación a ecuación de difusión-convección permite incluir los parámetros de flujo (velocidad) y difusión dependientes del espacio y de la aleatoriedad. La dinámica de los patrones de los coeficientes del polinomio de caos está bajo la influencia del tipo de método de solución. El uso del método de Bubnov Galerkin deja al descubierto la obtención de patrones con alta dispersión de incertidumbre en todo el dominio, mientras que el método de Petrov-Galerkin permite obtener zonas localizadas. En general, la localización de la incertidumbre se observa en aquellos lugares de alto gradiente en la respuesta determinista. Estas zonas se pueden captar con mayor precisión mediante la refinación de la malla o mediante una técnica de estabilización como SUPG.

## Bibliografía

1. [1] I. Babuska, F. Ihlenburg, E. Paik, S. Sauter; A generalized finite element method for solving the Helmholtz equation in two dimensions with minimal pollution; Comput. Methods Appl. Mech. Engrg., 128 (1995), pp. 325–359
2. [2] D. Garzon-Alvarado, J. García-Aznar, M. Doblare; Appearance and location of secondary ossification centres may be explained by a reaction-diffusion mechanism; Comput. Biol. Med., 39 (6) (2009), pp. 554–561
3. [3] S. Ferreira, M. Martins, M. Vilela; Reaction-diffusion model for the growth of avascular tumor; Phys. Rev., 65 (2) (2002)
4. [4] M. Chaplain, A. Ganesh, I. Graham; Spatio-temporal pattern formation on spherical surfaces: Numerical simulation and application to solid tumor growth; J Math Biol, 42 (2001), pp. 387–423
5. [5] A. Madzvamuse; A numerical approach to the study of spatial pattern formation in the ligaments of arcoid bivalves; B. Math. Biol., 64 (2002), pp. 501–530
6. [6] S. Kondo, R. Asai; A reaction-diffusion wave on the skin of the marine anglefish, Pomacanthus; Nature, 376 (1995), pp. 765–768
7. [7] F. Crauste, M. Lhassan, A. Kacha; A delay reaction-diffusion model of the dynamics of botulinum in fish; Math. Biosci., 216 (2008), pp. 17–29
8. [8] F. Rossi, S. Ristori, M. Rustici, N. Marchettini, E. Tiezzi; Dynamics of pattern formation in biomimetic systems; J. Theor. Biol., 255 (2008), pp. 404–412
9. [9] B. Rothschild, J. Ault; Population-dynamic instability as a cause of patch structure; Ecol. Model., 93 (1996), pp. 237–239
10. [10] T. Nozakura, S. Ikeuchi; Formation of dissipative structures in galaxies; Astrophys J, 279 (1984), pp. 40–52
11. [11] R. Smith; Optimal and near-optimal advection-diffusion finite-difference schemes iii. Black-Scholes equation; R. Soc. Lond. Proc. Ser. A Math. Phys. Eng. Sci., 456 (2000), pp. 1019–1028
12. [12] O. Richter; Modelling dispersal of populations and genetic information by finite element methods; Environ. Modell. Softw., 23 (2) (2008), pp. 206–214
13. [13] L. Ferragut, M. Asensio, S. Monedero; A numerical method for solving convection–reaction–diffusion multivalued equations in fire spread modelling; Advances in Engineering Software, 38 (2007), pp. 366–371
14. [14] O. Zienkiewicz, R. Taylor; Problemas de convección dominante: aproximaciones de elementos finitos a la ecuación de difusión-advección; Finite Element Method, 3, Ed. Butterworth-Heinemann College, Oxford (2000), pp. 5–150O. Zienkiewicz, R. Löhner, K. Morgan, S. Nakazawa Finite elements in fluid mechanics- a decade of progress; S. Comp. Phys., 5 (1984), pp. 1–26
15. [15] F. Brezzi, D. Marini, A. Russo; Applications of the pseudo residual-free bubbles to the stabilization of convection-diffusion problems; Comput. Methods Appl. Mech. Engrg., 166 (1998), pp. 51–63E. Oñate Derivation of stabilized equations for numerical solution of advective diffusive transport and fluid flow problems; Comput. Methods Appl. Mech. Engrg., 151 (1998), pp. 233–265
16. [16] J. Chrispell, V. Ervin, E. Jenkins; A fractional step θ-method for convection–diffusion problems; J. Math. Anal. Appl., 333 (2007), pp. 204–218
17. [17] T. Hughes, L. Franca, G. Hulbert, Z. Johan, F. Sakhib; The Galerkin least square method for advective diffusion equations; AMD 94-ASME, 94 (1988), pp. 1–20
18. [18] O. Zienkiewicz, R. Gallagher, P. Hood; Newtonian and non-newtonian viscous impompressible flow; Temperature induced flows and finite elements solutions The mathematics of finite elements and applications, Ed. Academic Press, London (1975)
19. [19] I. Christie, D. Griffiths, O. Zienkiewicz; Finite element methods for second order differential equations with significant first derivatives; Internat. J. Numer. Methods Engrg., 10 (1976), pp. 1389–1396
20. [20] O. Zienkiewicz, J. Heinrich, P. Huyakorn, A. Mitchel; An upwind finite element scheme for two dimensional convective transport equations; Internat. J. Numer. Methods Engrg., 11 (1977), pp. 131–144
21. [21] T.J.R. Hughes, editor, Finite element methods for convection dominated flows, AMD 3 (1979) 19-35. ASME, New York.
22. [22] T. Hughes, L. Franca, G. Hulbert; A new finite element formulation for computational fluid dynamics: VIII. The Galerkin/least-squares method for advective diffusive equationsJT Comput. Methods Appl. Mech. Engrg; 73 (1989), pp. 173–189
23. [23] C. Baiocchi, F. Brezzi, L. Franca; Virtual bubbles and Galerkin/least-square type methodsJT Comput. Methods Appl. Mech. Engrg; 105 (1993), pp. 125–141
24. [24] I. Harari, T. Hughes; Stabilized finite element methods for steady advection-diffusion with production; Comput. Methods Appl. Mech. Engrg., 115 (1994), pp. 165–191
25. [25] A. Russo; Bubble stabilization of the finite element method for the linearized incompressible Navier-Stokes equation; Comput. Methods Appl. Mech. Engrg., 132 (1996), pp. 335–343
26. [26] B. Cockburn, G. Karniadakis, C. Shu; Discontinuos Galerkin methods, theory, computational and application; Ed. Springer, Berlin (2000)
27. [27] R. Araya, E. Behrens, R. Rodríguez; An adaptive stabilized finite element scheme for the advection–reaction–diffusion equation; Appl. Numer. Math., 54 (2005), pp. 491–503
28. [28] H. Wang, H. Dahle, R. Ewing, M. Espedal, R. Sharpley, S. Man; An ELLAM scheme for advection–diffusion equations in two dimensions; SIAM J Sci Comput, 20–6 (1979), pp. 2160–2194
29. [29] M. Berzins; Temporal error control for convection-dominated equations in two space dimensions; SIAM J Sci Comput, 16–3 (1995), pp. 558–580
30. [30] R. Ghanem, P. Spanos; Stochastic finite elements. A spectral approach; Dover Publications., Minneola, New York, USA (1991)
31. [31] E.W. Boyce, B.E. Goodwin; Random transverse vibration of elastic beams; SIAM J, 12–3 (1964), pp. 613–629
32. [32] J.D. Collins, W.T. Thompson; The eigenvalue problem for structural systems with uncertain parameters; AIAA J, 7–4 (1969), pp. 642–648
33. [33] G.C. Hart, J.D. Collins; The treatment of randomness in finite element modelling. SAE Shock Vibrations Symposium; Los Angeles, CA (1970), pp. 2509–2519
34. [34] T.K. Hasselman, G.C. Hart; Modal analysis of random structural systems; J. Eng. Mech., 98–EM3 (1972), pp. 561–579
35. [35] W.K. Liu, G. Besterfield, A. Mani; Probabilistic finite element methods in nonlinear structural dynamics; Comput. Methods Appl. Mech. Engrg., 57 (1986), pp. 61–81
36. [36] W.K. Liu, G. Besterfield, T. Belytschko; Transient probabilistic systems; Comput. Methods Appl. Mech. Engrg., 67 (1988), pp. 27–54
37. [37] M. Kaminski; Generalized perturbation-based stochastic finite element method in elastostatics; Comput. Struct., 85 (2007), pp. 586–594
38. [38] S. Lepage; Stochastic finite element method for the modeling of thermoelastic damping in micro-resonators. Ph. D. Thesis; Cranfield University (2006)
39. [39] L. Yao, P. He, S. Song; A perturbation stochastic finite-element method for groundwater flow models based on an undetermined-coefficients approach; Hydrogeol. J., 18–7 (2010), pp. 1603–1609
40. [40] A. Chaudhuri, M. Sekhar. Stochastic finite element method for probabilistic analysis of flow and transport in a three-dimensional heterogeneous porous formation. Water Resources Research. 41 (2005) 1–15.
41. [41] R.G. Ghanem; Stochastic finite elements with multiple random non-Gaussian properties; J Eng Mech ASCE, 125–1 (1999), pp. 26–40
42. [42] D. Xiu, G. Karniadakis; The Wiener-Askey polynomial chaos for stochastic differential equations; J Sci Comput, 24–2 (2002), pp. 619–644
43. [43] O. Le Maitre, O.M. Knio, H.N. Najm, R. Ghanem; A stochastic projection method for fluid flow; J Comp Phys, 173 (2001), pp. 481–511
44. [44] B. Van den Nieuwenhof; Stochastic finite elements for elastodynamics: Random fields and shape uncertainty modeling using direct and modal perturbation-based approaches. Ph. D. Thesis; Université Catholique de Louvain, Belgium (2003)
45. [45] D. Clouteau, R. Lafargue. An iterative solver for stochastic soil-structure interaction. In: P. Spanos, G. Deodatis, ed. Computational Stochastic Mechanics. Millpress, Rotterdam, 2003. p. 119-124.
46. [46] N.Z. Chen, C. Guedes Soares; Spectral stochastic finite element analysis for laminated composite planes; Comput. Methods Appl. Mech. Engrg., 197 (51-52) (2008), pp. 4830–4839
47. [47] K. Maute, G. Weickum, M. Eldred; A reduced-order stochastic finite element approach for design optimization under uncertainty; Struct. Saf., 31 (6) (2009), pp. 450–459
48. [48] B. Moller, W. Graf, M. Beer, J. Sickert; Fuzzy stochastic finite element method. Computational fluid and solid mechanics 2003; Proceedings second MIT conference on computational fluid and solid mechanics, June 17-20 (2003), pp. 2074–2077
49. [49] B. Nicolai, J. Egea, N. Scheerlinck, J. Banga, A. Datta; Fuzzy finite element analysis of heat conduction problems with uncertain parameters; J. Food Eng., 103 (1) (2011), pp. 38–46
50. [50] A. Nouy, A. Clement, F. Schoefs, N. Moes; An extended stochastic finite element method for solving stochastic partial differential equations on random domains; Comput. Methods Appl. Mech. Engrg., 197 (51–52) (2008), pp. 4663–4682
51. [51] Z. Lei, C. Qiu; Neumann dynamic stochastic finite element method of vibration for structures with stochastic parameters to random excitation; Comput. Struct., 77 (6) (2000), pp. 651–657
52. [52] F. Xu; A multiscale stochastic finite element method on elliptic problems involving uncertainties; Comput. Methods Appl. Mech. Engrg., 196 (25–28) (2007), pp. 2723–2736
53. [53] P. Malliavin; Stochastic analysis; Springer Verlag, New York (1997)
54. [54] M. Loève; Probability theory; (4th ed)Springer, New York (1977)
55. [55] R. Cameron, W. Martin; The orthogonal development of nonlinear functional in series of Fourier-Hermite functionals; Ann Math, 48 (1947), pp. 385–392
56. [56] S.K. Choi, R.V. Grandhi, R.A. Canfield; Reliability-based Structural design; Springer Verlag, London, UK (2007)
57. [57] B. Sudret, M. Berveiller; Application of a stochastic finite element procedure to reliability analysis; Proc. 11th IFIP WG7.5 Conference, Banff, Canada (2003)

### Document information

Published on 01/12/13
Accepted on 11/02/13
Submitted on 11/03/11

Volume 29, Issue 4, 2013
DOI: 10.1016/j.rimni.2013.07.001
Licence: Other

### Document Score

0

Views 11
Recommendations 0