Solution on square domains of reaction-convection-diffusion equations using spectral stochastic finite element and streamline upwind Galerkin Petrov

==Resumen==

Se desarrolla la solución numérica de 2 problemas de difusión-convección (DC), mediante el empleo del método de elementos finitos de Petrov-Galerkin en contracorriente (SUPG). Los parámetros que definen el comportamiento de las ecuaciones se modelan como campos estocásticos, y son los siguientes: la velocidad convectiva, la difusión y la capacidad calorífica como valores de tipo aleatorio. Por tanto, se combina el método SUPG para ecuaciones de DC, con convección dominante, con el método de los elementos finitos estocásticos espectrales. Los parámetros de cada ecuación se han descrito mediante la expansión de Karhunen-Loève, mientras que la incógnita se representa mediante la expansión de los polinomios de caos. Los objetivos del artículo son: en primer lugar, estudiar la influencia de los campos estocásticos en la solución de problemas de DC con SUPG, y en segundo lugar obtener los patrones de cada coeficiente de la expansión en polinomios de caos. Los resultados muestran la versatilidad del método para solucionar diferentes problemas físicos gracias a la generalidad en la descripción estadística de la solución, y la riqueza en la representación de las zonas donde se halla la mayor variabilidad en la respuesta. Los patrones muestran 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.

==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 [[#bib0005|[1]]] , [[#bib0010|[2]]] , [[#bib0015|[3]]] , [[#bib0020|[4]]] , [[#bib0025|[5]]] , [[#bib0030|[6]]] , [[#bib0035|[7]]] , [[#bib0040|[8]]] , [[#bib0045|[9]]] , [[#bib0050|[10]]] , [[#bib0055|[11]]] , [[#bib0060|[12]]]  and [[#bib0065|[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).

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\frac{\partial \varphi }{\partial t}-\nabla \cdot (\boldsymbol{\mbox{k}}\nabla \varphi )+</math><math>\boldsymbol{\mbox{u}}\cdot \nabla \varphi -s\varphi =</math><math>f(\boldsymbol{\mbox{x}})\mbox{ }en\mbox{ }\Omega </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 1)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\varphi (t=0)={\varphi }_o</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 2)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\varphi (\boldsymbol{\mbox{x}})=g(\boldsymbol{\mbox{x}})\mbox{ }\mbox{en}\mbox{ }{\Gamma }_{\varphi }</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 3)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\overset{\rightarrow}{\nabla }(\varphi )\bullet \boldsymbol{\mbox{n}}=</math><math>h(\boldsymbol{\mbox{x}})\mbox{ }\mbox{en}\mbox{ }{\Gamma }_{\nabla }</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 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), <math display="inline">f(\boldsymbol{\mbox{x}})</math>  es la función de generación, <math display="inline">g(\boldsymbol{\mbox{x}})</math>  es la función que define el valor del campo escalar <math display="inline">\varphi </math>  sobre la frontera <math display="inline">{\Gamma }_{\varphi }</math>  y <math display="inline">h(\boldsymbol{\mbox{x}})</math>  es la función que define el valor del flujo sobre la frontera <math display="inline">{\Gamma }_{\nabla }</math> , cuya normal es '''n''' .      

El comportamiento y los patrones de la solución de los ejemplos señalados [[#bib0005|[1]]] , [[#bib0010|[2]]] , [[#bib0015|[3]]] , [[#bib0020|[4]]] , [[#bib0025|[5]]] , [[#bib0030|[6]]] , [[#bib0035|[7]]] , [[#bib0040|[8]]] , [[#bib0045|[9]]] , [[#bib0050|[10]]] , [[#bib0055|[11]]] , [[#bib0060|[12]]]  and [[#bib0065|[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 <math display="inline">s</math> . En casos más generales, en ecuaciones de reacción-difusión-convección, donde <math display="inline">s</math>  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 [[#bib0010|[2]]]  and [[#bib0025|[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 [[#bib0070|[14]]] , el método de cálculo finito [[#bib0075|[15]]] , el método de paso fraccional <math display="inline">\theta </math>[[#bib0080|[16]]] , el método de mínimos cuadrados de Galerkin [[#bib0085|[17]]]  y el método Petrov-Galerkin de contracorriente [[#bib0090|[18]]] , [[#bib0095|[19]]] , [[#bib0100|[20]]] , [[#bib0105|[21]]] , [[#bib0110|[22]]] , [[#bib0115|[23]]] , [[#bib0120|[24]]] , [[#bib0125|[25]]] , [[#bib0130|[26]]] , [[#bib0135|[27]]] , [[#bib0140|[28]]]  and [[#bib0145|[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 [[#bib0090|[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 [[#bib0150|[30]]] .      

El método de la perturbación es el utilizado con más frecuencia debido a su simplicidad matemática [[#bib0150|[30]]] . Este método consiste en expandir todas las cantidades aleatorias alrededor de su media utilizando las series de Taylor. El principal inconveniente del método surge de aquellos términos que causan inestabilidad que, en general, son los términos de alto orden. Además, también de manera general, los términos que siguen al primer y segundo componente de la serie de Taylor no son de gran utilidad práctica, por lo cual se limita la respuesta. Por esta razón, el método es utilizado para aquellos problemas en los cuales la respuesta de los parámetros y las variables tienen una «pequeña» amplitud aleatoria alrededor de la media [[#bib0150|[30]]] . Las aplicaciones del método de la perturbación han sido frecuentes, y se han utilizado en ingeniería estructural y en modelamiento matemático. Por ejemplo, Boyce y Goodwin [[#bib0155|[31]]]  utilizan el método para hallar la solución de autovalores y autovectores para vigas y cuerdas sometidas a vibración aleatoria. Collins y Thompson [[#bib0160|[32]]]  utilizan la metodología de perturbaciones para calcular los autovalores en problemas estructurales que tienen campos aleatorios que son descritos por su matriz de covarianza. Hart y Collins [[#bib0165|[33]]]  y Hasselman y Hart [[#bib0170|[34]]]  hacen un desarrollo compatible del método de las perturbaciones y del método de los elementos finitos. Este importante avance permite introducir las características aleatorias de las variables en el método de los elementos finitos (MEF); es así como Liu et al. [[#bib0175|[35]]]  introducen de manera general el método de las perturbaciones en problemas lineales y no lineales para mecánica de sólidos. En Liu et al. [[#bib0180|[36]]]  se amplía el método para problemas transientes de mecánica de sólidos. Recientemente, Kaminski [[#bib0185|[37]]]  ha hecho un desarrollo general del método de la perturbación para problemas de orden superior, con lo que se supera el problema de los momentos estadísticos que generan inestabilidad. Igualmente Lepage [[#bib0190|[38]]]  utiliza el método conjugado de perturbación y elementos finitos, ''Perturbation Stochastic Finite Element Method'' , denominado ahora PSFEM, para modelar el amortiguamiento termoelástico para el diseño de microrresonadores. Además de las aplicaciones señaladas se puede encontrar trabajos con la técnica de perturbación en el análisis de flujos de canales abiertos  [[#bib0195|[39]]] , en medios porosos [[#bib0200|[40]]]  y en diseño estructural [[#bib0175|[35]]] , entre otros.      

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 [[#bib0150|[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 [[#bib0205|[41]]]  and [[#bib0210|[42]]]  y campos lognormales [[#bib0205|[41]]] . Desarrollos recientes han mostrado la versatilidad del método para diversas aplicaciones, entre las que se incluyen: mecánica de fluidos [[#bib0215|[43]]] , análisis modal [[#bib0220|[44]]] , interacción fluido estructura [[#bib0225|[45]]] , análisis de placas laminadas [[#bib0230|[46]]]  y optimización [[#bib0235|[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 [[#bib0150|[30]]] , [[#bib0155|[31]]] , [[#bib0160|[32]]] , [[#bib0165|[33]]] , [[#bib0170|[34]]] , [[#bib0175|[35]]] , [[#bib0180|[36]]] , [[#bib0185|[37]]] , [[#bib0190|[38]]] , [[#bib0195|[39]]] , [[#bib0200|[40]]] , [[#bib0205|[41]]] , [[#bib0210|[42]]] , [[#bib0215|[43]]] , [[#bib0220|[44]]] , [[#bib0225|[45]]] , [[#bib0230|[46]]]  and [[#bib0235|[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 [[#bib0150|[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 [[#bib0240|[48]]]  y transferencia de calor [[#bib0245|[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 [[#bib0250|[50]]] , el método de los elementos finitos estocásticos dinámicos de Neumann [[#bib0255|[51]]] , o el método de elementos finitos estocásticos multiescala [[#bib0260|[52]]] .      

Todos los métodos aquí mencionados se han comparado típicamente con el método de Monte Carlo [[#bib0150|[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 [[#bib0150|[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 [[#bib0265|[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 [[#bib0105|[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 (<math display="inline">W_i</math> ), tal como se muestra en (5) [[#bib0030|[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)  [[#bib0035|[7]]] :

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>W_i^\ast =N_i+\frac{\alpha h}{2}\frac{dN_i}{dx}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 5)
|}

donde:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\alpha =min\left(Coth\left|P_e\right|-\frac{1}{\left|P_e\right|}\mbox{ },\mbox{ }1-\right. </math><math>\left. \frac{1}{\left|P_e\right|}\right)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 6)
|}

En esta última expresión <math display="inline">P_e</math>  es el número adimensional de Peclet definido como <math display="inline">P_e=</math><math>\frac{\left\|\boldsymbol{u}\right\|h}{2k}</math> . Nótese que la función de peso solo debe ser perturbada si <math display="inline">P_e>1</math> . 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 <math display="inline">\left[K\right]</math>  de un elemento ''e'' , calculados de acuerdo con el método SUPG unidimensional, están definidos por medio de la expresión (7)  [[#bib0070|[14]]] :

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>K_{lm}^e=\int_e\frac{dW_l^\ast }{dx}k_{im}\frac{dN_m}{dx}dx+</math><math>\int_eW_l^\ast u\frac{dN_m}{dx}dx</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 7)
|}

En esta última expresión, ''N''<sub>''m''</sub>  son las funciones de forma empleadas para aproximar la solución [[#bib0070|[14]]] , tal como se observa en (8):

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\varphi }^e=\sum_mN_m{\varphi }_m</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 8)
|}

De forma análoga, los términos del vector <math display="inline">\left[F\right]</math>  para un elemento <math display="inline">e</math>  se definen de acuerdo con (9) [[#bib0070|[14]]] :

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>F_l^e=\int_eW_l^\ast f(x)dx</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 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) [[#bib0105|[21]]] :

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>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. </math><math>\left. \frac{u_y}{\left\|\boldsymbol{\mbox{u}}\right\|}\left(\frac{\partial N_i}{\partial y}\right)\right]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 10)
|}

donde <math display="inline">u_x</math>  y <math display="inline">u_y</math>  son las componentes de la velocidad en las direcciones globales <math display="inline">x-</math><math>y</math> , en tanto que α es el coeficiente de perturbación calculado de acuerdo con (6), y <math display="inline">h</math>  es la longitud característica del elemento, como se muestra en la [[#fig0005|figura 1]] .

<span id='fig0005'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_García_865974548-1-s2.0-S0213131513000357-gr1.jpg|center|320px|Línea de corriente al interior de un elemento bidimensional.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 1.

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

</span>
|}

===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):

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\int_{\Omega }{W_l}^\ast \frac{\partial \phi }{\partial t}d\Omega +</math><math>\int_{\Omega }\overset{\rightarrow}{\nabla }{W_l}^\ast \cdot \boldsymbol{\mbox{k}}\overset{\rightarrow}{\nabla }\phi d\Omega +</math><math>\int_{\Omega }{W_l}^\ast \boldsymbol{\mbox{u}}\cdot \overset{\rightarrow}{\nabla }\phi d\Omega =</math><math>\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 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 11)
|}

donde <math display="inline">l=1,...,nnod</math> , <math display="inline">\overset{\rightarrow}{n}</math>  es el vector normal al borde de flujo o borde de Neumman y <math display="inline">nnod</math>  es el número de nodos del dominio discretizado. Empleando ahora aproximaciones por tramos del tipo <math display="inline">{\varphi }^e=</math><math>\sum_{m=1}^{nnod}N_m{\varphi }_m</math> , así como la función de peso modificada (10), se puede conducir de la ecuación global a una expresión a nivel elemental [[#bib0070|[14]]] , como la mostrada en (12):

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{\mbox{M}}\frac{\Delta \boldsymbol{\phi }}{\Delta t}+</math><math>\boldsymbol{\mbox{K}}\boldsymbol{\phi }+\boldsymbol{\mbox{C}}\boldsymbol{\phi }=</math><math>\boldsymbol{\mbox{F}}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 12)
|}

donde:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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)}^Tfd{\Omega }^e+c.c.\\
\Delta \boldsymbol{\phi }={\boldsymbol{\phi }}^{t+\Delta t}-{\boldsymbol{\phi }}^t\\
\boldsymbol{\phi }={\boldsymbol{\phi }}^{t+\Delta t}
\end{array}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | 
|}

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

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | 
|}

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:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\left[\int_{{\Omega }^e}{\left(\boldsymbol{\mbox{W}}^\boldsymbol{\ast }\right)}^T1(\boldsymbol{\mbox{x}},\omega )\boldsymbol{\mbox{N}}d\Omega \right]\frac{\Delta \boldsymbol{\phi }}{\Delta t}+</math><math>\left[\int_{{\Omega }^e}{\left(\nabla \boldsymbol{\mbox{W}}^\boldsymbol{\ast }\right)}^Tk(\boldsymbol{\mbox{x}},\omega )\nabla \boldsymbol{\mbox{N}}d\Omega \right]\boldsymbol{\phi }+</math><math>\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 }=</math><math>\boldsymbol{\mbox{F}}(x,\omega )</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 13)
|}

donde los parámetros del modelo (1, <math display="inline">k</math> , <math display="inline">\overset{\rightarrow}{u}</math>  y '''F''' ) dependen del espacio y poseen naturaleza aleatoria, la cual está indicada por el argumento <math display="inline">\omega </math>[[#bib0205|[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 [[#bib0270|[54]]]  de un proceso estocástico, por ejemplo <math display="inline">k(\boldsymbol{\mbox{x}},\omega )</math> , se basa en la expansión espectral de su función de covarianza <math display="inline">R_{EE}(\boldsymbol{\mbox{x}}_\boldsymbol{1},\boldsymbol{\mbox{x}}_\boldsymbol{2})</math> , donde <math display="inline">\boldsymbol{\mbox{x}}_\boldsymbol{1}</math>  y <math display="inline">\boldsymbol{\mbox{x}}_\boldsymbol{2}</math>  denotan las coordenadas espaciales y <math display="inline">\omega </math>  denota la aleatoriedad. Por definición, la función de covarianza es simétrica y definida positiva [[#bib0150|[30]]]  and [[#bib0270|[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 <math display="inline">k(\boldsymbol{\mbox{x}},\omega )</math>  toman la forma (14):

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>k(\boldsymbol{\mbox{x}},\omega )=\overline{k}(\boldsymbol{\mbox{x}})+</math><math>\sum_{j=1}^{\infty }\sqrt{{\lambda }_j}k_j(\boldsymbol{\mbox{x}}){\xi }_j(\omega )</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 14)
|}

donde <math display="inline">\overline{k}(\boldsymbol{\mbox{x}})</math>  es el valor promedio del proceso estocástico, en este caso, del coeficiente de difusión, los términos <math display="inline">\left\{{\xi }_j(\omega )\right\}</math>  forman un conjunto ortogonal de variables aleatorias, <math display="inline">\left\{{\lambda }_j\right\}</math>  es el conjunto de autovalores y <math display="inline">\left\{k_j(\boldsymbol{\mbox{x}})\right\}</math>  es el conjunto de autofunciones de la función de covarianza, los cuales pueden ser evaluados como la solución de la siguiente integral:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\int_{\Omega }R_{EE}(\boldsymbol{\mbox{x}},\boldsymbol{\mbox{y}})k_j(\boldsymbol{\mbox{y}})d\boldsymbol{\mbox{y}}=</math><math>{\lambda }_jk_j(\boldsymbol{\mbox{x}})</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 15)
|}

donde <math display="inline">\Omega </math>  denota el dominio espacial sobre el cual se presenta el proceso <math display="inline">k(\boldsymbol{\mbox{x}},\omega )</math> . 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 [[#bib0230|[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 <math display="inline">\left\{{\xi }_j(\omega )\right\}</math>  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  [[#bib0150|[30]]] .      

En el problema aquí desarrollado las variables estocásticas (1, <math display="inline">k</math> , <math display="inline">\overset{\rightarrow}{\boldsymbol{\mbox{u}}}</math>  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 [[#bib0205|[41]]] :

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>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 )</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 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 [[#bib0150|[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):

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{\mbox{k}}(\boldsymbol{\mbox{x}},\omega )=</math><math>\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 )</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 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):

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 18)
|}

Las funciones <math display="inline">1_{j_1}(\boldsymbol{\mbox{x}})</math> , <math display="inline">k_{j_k}(\boldsymbol{\mbox{x}})</math> , <math display="inline">v_{xj_x}</math>  y <math display="inline">v_{yj_y}(\boldsymbol{\mbox{x}})</math>  corresponden a las autofunciones de la función de covarianza de la variable de campo. Igualmente, los valores <math display="inline">{\lambda }_{j_1}^1</math> , <math display="inline">{\lambda }_{j_k}^k</math> , <math display="inline">{\lambda }_{j_x}^{v_x}</math>  y <math display="inline">{\lambda }_{j_y}^{v_y}</math>  son los autovalores. De esta forma, cada conjunto de autofunciones y autovalores conforman el denominado autopar [[#bib0150|[30]]] .      

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

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\left\{\left[\int_{{\Omega }^e}{\left(\boldsymbol{\mbox{W}}^\boldsymbol{\ast }\right)}^T1\boldsymbol{\mbox{N}}d\Omega \right]+\right. </math><math>\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}+</math><math>\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. </math><math>\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 }+</math><math>\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. </math><math>\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 }=</math><math>\boldsymbol{\mbox{F}}(x,\omega )</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 19)
|}

Al definir <math display="inline">{\xi }_0=1</math> , <math display="inline">{\lambda }_0^1=</math><math>1</math> , <math display="inline">{\lambda }_0^k=1</math> , <math display="inline">{\lambda }_0^{v_x}=</math><math>1</math> , <math display="inline">{\lambda }_0^{v_y}=</math><math>1</math> , <math display="inline">1_0(\boldsymbol{\mbox{x}})=</math><math>1</math> , <math display="inline">k_0(\boldsymbol{\mbox{x}})=</math><math>\overline{k}(\boldsymbol{\mbox{x}})</math> , <math display="inline">v_{x_0}(\boldsymbol{\mbox{x}})=</math><math>\overline{v_x}(\boldsymbol{\mbox{x}})</math>  y <math display="inline">v_{y_0}(\boldsymbol{\mbox{x}})=</math><math>\overline{v_y}(\boldsymbol{\mbox{x}})</math>  y definir vectorialmente la velocidad (<math display="inline">\boldsymbol{\mbox{u}}(\boldsymbol{\mbox{x}})=</math><math>{\left[v_x(\boldsymbol{\mbox{x}}),v_y(\boldsymbol{\mbox{x}})\right]}^T</math> ) se llega a (20):

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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}+</math><math>\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 }+</math><math>\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 }=</math><math>\boldsymbol{\mbox{F}}(x,\omega )</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 20)
|}

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

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\left\{\sum_{j_1=0}^{\infty }\boldsymbol{\mbox{M}}_{\boldsymbol{\mbox{j}}_\boldsymbol{1}}{\xi }_{k_1}\right\}\frac{\Delta \boldsymbol{\phi }}{\Delta t}+</math><math>\left\{\sum_{j_k=0}^{\infty }\boldsymbol{\mbox{K}}_{\boldsymbol{\mbox{j}}_\boldsymbol{\mbox{k}}}{\xi }_{j_k}\right\}\boldsymbol{\phi }+</math><math>\left\{\sum_{j_u=0}^{\infty }\boldsymbol{\mbox{C}}_{\boldsymbol{\mbox{j}}_\boldsymbol{\mbox{u}}}{\xi }_{j_u}\right\}\boldsymbol{\phi }=</math><math>\boldsymbol{\mbox{F}}(x,\omega )</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 21)
|}

donde:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | 
|}

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

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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}+</math><math>\left\{\sum_{j_k=0}^{\infty }\boldsymbol{\mbox{K}}_{\boldsymbol{\mbox{j}}_\boldsymbol{\mbox{k}}}{\xi }_{j_k}\right\}{\boldsymbol{\phi }}^{t+\Delta t}+</math><math>\left\{\sum_{j_u=0}^{\infty }\boldsymbol{\mbox{C}}_{\boldsymbol{j}_\boldsymbol{u}}{\xi }_{j_u}\right\}{\boldsymbol{\phi }}^{t+\Delta t}=</math><math>\boldsymbol{\mbox{F}}(x,\omega )</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 22)
|}

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

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\left[\frac{1}{\Delta t}\left\{\sum_{j_1=0}^{\infty }\boldsymbol{\mbox{M}}_{\boldsymbol{\mbox{j}}_\boldsymbol{1}}{\xi }_{k_1}\right\}+\right. </math><math>\left. \left\{\sum_{j_k=0}^{\infty }\boldsymbol{\mbox{K}}_{\boldsymbol{\mbox{j}}_\boldsymbol{\mbox{k}}}{\xi }_{j_k}\right\}+\right. </math><math>\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}=</math><math>\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</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 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 [[#bib0220|[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 <math display="inline">\left\{{\xi }_j(\omega )\right\}</math> . Se ha demostrado [[#bib0275|[55]]]  que la dependencia funcional de las variables aleatorias puede ser expresada en términos de polinomios de caos, descritos como (24):

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\phi (\omega )=a_0{\Gamma }_0+\sum_{i=1}^{\infty }a_i{\Gamma }_1({\xi }_i(\omega ))+</math><math>\sum_{j=1}^{\infty }\sum_{i=1}^{\infty }a_{ij}{\Gamma }_2({\xi }_i(\omega ),{\xi }_j(\omega ))+</math><math>...\mbox{.}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 24)
|}

en donde <math display="inline">{\Gamma }_n({\xi }_{i_1}(\omega ),{\xi }_{i_2}(\omega ),...,{\xi }_{i_n}(\omega ))</math>  son los polinomios de caos de orden ''n  ''  en las variables <math display="inline">\left\{{\xi }_j(\omega )\right\}</math> . Estos polinomios son la generalización de los polinomios multidimensionales de Hermite [[#bib0150|[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) [[#bib0150|[30]]] . Si se desarrollan los polinomios de caos, se puede obtener un conjunto de polinomios bien ordenados indicialmente descritos por <math display="inline">{\Psi }_i(\omega )</math>[[#bib0150|[30]]]  and [[#bib0205|[41]]] , de tal forma que (24) se puede reescribir como (25):

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\varphi (\omega )=\sum_{i=1}^{\infty }{\varphi }_i{\Psi }_i(\omega )</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 25)
|}

donde la variable <math display="inline">\varphi (\omega )</math>  representa un valor en «un» punto nodal. Estos polinomios (<math display="inline">{\Psi }_i(\omega )</math> ) son mutuamente ortogonales en el sentido del producto interno <math display="inline">\left\langle {\Psi }_i(\omega ),{\Psi }_j(\omega )\right\rangle</math> , el cual está dado por (26):

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\overline{{\Psi }_i(\omega ){\Psi }_j(\omega )}=\left\langle {\Psi }_i(\omega ),{\Psi }_j(\omega )\right\rangle</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 26)
|}

donde puede apreciarse que el producto interno está definido como el promedio (esperanza estadística) del producto de los 2 polinomios [[#bib0205|[41]]] . Este producto interno será cero en el caso en que <math display="inline">i\not =</math><math>j</math> . La caracterización completa del proceso aleatorio tiene lugar cuando se encuentran los coeficientes deterministas, dados por <math display="inline">{\varphi }_i</math>[[#bib0150|[30]]] .      

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

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\left[\frac{1}{\Delta t}\left\{\sum_{j_1=0}^{\infty }\boldsymbol{\mbox{M}}_{\boldsymbol{\mbox{j}}_\boldsymbol{1}}{\xi }_{k_1}\right\}+\right. </math><math>\left. \left\{\sum_{j_k=0}^{\infty }\boldsymbol{\mbox{K}}_{\boldsymbol{\mbox{j}}_\boldsymbol{\mbox{k}}}{\xi }_{j_k}\right\}+\right. </math><math>\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)=</math><math>\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)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 27)
|}

donde el vector <math display="inline">{\boldsymbol{\phi }}_\boldsymbol{\mbox{i}}=</math><math>[{\varphi }_{i1},{\varphi }_{i2},....,{\varphi }_{innodel}]</math>  representa la i-ésima incógnita estocástica.      

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

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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. </math><math>\left. \left\{\sum_{j_k=0}^{\infty }\boldsymbol{\mbox{K}}_{\boldsymbol{\mbox{j}}_\boldsymbol{\mbox{k}}}{\xi }_{j_k}\right\}+\right. </math><math>\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)-</math><math>\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)=</math><math>0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 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):

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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. </math><math>\left. \left\{\sum_{j_k=0}^{M_k}\boldsymbol{\mbox{K}}_{\boldsymbol{\mbox{j}}_\boldsymbol{\mbox{k}}}{\xi }_{j_k}\right\}+\right. </math><math>\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)-</math><math>\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</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 29)
|}

donde <math display="inline">M_1</math> , <math display="inline">M_K</math>  y <math display="inline">M_u</math>  determinan el truncamiento de la expansión de Karhunen-Loève para cada uno de los campos aleatorios. De igual forma, <math display="inline">P</math>  determina el truncamiento de la expansión de los polinomios de caos [[#bib0150|[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 [[#bib0150|[30]]] :

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\left\langle \boldsymbol{\mbox{r}},{\Psi }_l(\omega )\right\rangle=</math><math>0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 30)
|}

Para <math display="inline">l=0,1,2,...,P</math> . Mediante manipulación algebraica se obtiene la ecuación (31):

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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. </math><math>\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. </math><math>\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. </math><math>\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. </math><math>\left. \left\langle \boldsymbol{\mbox{F}}(x,\omega ){\Psi }_l(\omega )\right\rangle=\right. </math><math>\left. 0\right]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 31)
|}

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

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\left\langle {\xi }_{k_1}{\Psi }_i(\omega ){\Psi }_l(\omega )\right\rangle=</math><math>\overline{{\xi }_{k_1}{\Psi }_i(\omega ){\Psi }_l(\omega )}=</math><math>d_{j_kil}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 32)
|}

donde el término <math display="inline">d_{j_kil}</math>  está dado por (33) [[#bib0280|[56]]]  and [[#bib0285|[57]]] :

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>d_{j_kil}=\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}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 33)
|}

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

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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_1il}{\boldsymbol{\phi }}_i^{t+\Delta t}\right\}+\right. </math><math>\left. \sum_{i=1}^P\left\{\sum_{j_k=0}^{M_k}\boldsymbol{\mbox{K}}_{\boldsymbol{\mbox{j}}_\boldsymbol{\mbox{k}}}d_{j_kil}{\boldsymbol{\phi }}_i^{t+\Delta t}\right\}+\right. </math><math>\left. \sum_{i=1}^P\left\{\sum_{j_u=0}^{M_u}\boldsymbol{\mbox{C}}_{\boldsymbol{\mbox{j}}_\boldsymbol{\mbox{u}}}d_{j_uil}{\boldsymbol{\phi }}_i^{t+\Delta t}\right\}-\right. </math><math>\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_1il}{\boldsymbol{\phi }}_i^t\right\}-\right. </math><math>\left. \left\langle \boldsymbol{\mbox{F}}(x,\omega ){\Psi }_l(\omega )\right\rangle=\right. </math><math>\left. 0\right]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 34)
|}

donde <math display="inline">l=0,1,2,...,P</math> .      

Expandiendo (34) y considerando que el término <math display="inline">\boldsymbol{F}(x,\omega )</math>  es determinista (esto es: <math display="inline">\left\langle \boldsymbol{\mbox{F}}(x,\omega ){\Psi }_l(\omega )\right\rangle=</math><math>\boldsymbol{\mbox{F}}_\boldsymbol{0}</math> ), la ecuación se convierte en el sistema (35):

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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_1i0}+\sum_{j_k=0}^{M_k}\boldsymbol{\mbox{K}}_{j_k}d_{j_ki0}+\sum_{j_u=0}^{M_u}\boldsymbol{\mbox{C}}_{j_u}d_{j_ui0}\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_1i0}{\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_1i1}+\sum_{j_k=0}^{M_k}\boldsymbol{\mbox{K}}_{j_k}d_{j_ki1}+\sum_{j_u=0}^{M_u}\boldsymbol{\mbox{C}}_{j_u}d_{j_ui1}\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_1i1}{\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_kis}+\sum_{j_u=0}^{M_u}\boldsymbol{\mbox{C}}_{j_u}d_{j_uis}\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_1is}{\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_kiP}+\sum_{j_u=0}^{M_u}\boldsymbol{\mbox{C}}_{j_u}d_{j_uiP}\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_1iP}{\boldsymbol{\phi }}_i^t\right\}
\end{array}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 35)
|}

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

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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_kis}+\sum_{j_u=0}^{M_u}\boldsymbol{\mbox{C}}_{\boldsymbol{\mbox{j}}_\boldsymbol{\mbox{u}}}d_{j_uis}\\
B_{is}=\frac{1}{\Delta t}\sum_{j_1=0}^{M_1}\boldsymbol{\mbox{M}}_{\boldsymbol{\mbox{j}}_\boldsymbol{1}}d_{j_1is}
\end{array}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 36)
|}

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

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 37)
|}

Y en términos matriciales (38):

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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)}}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 38)
|}

===2.7. Sobre los autovalores y las autofunciones===

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

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\int_{\Omega }R_{EE}(x_1,y_1;x_2,y_2)k_j(x_2,y_2)dx_2dy_2=</math><math>{\lambda }_jk_j(x_1,y_1)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 39)
|}

donde <math display="inline">R_{EE}(\boldsymbol{\mbox{x}},\boldsymbol{\mbox{y}})</math>  es la función de covarianza. En este trabajo se supone (40) [[#bib0285|[57]]] :

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>R_{EE}(x_1,y_1;x_2,y_2)=C_0exp\left(-\frac{\left|x_1-x_2\right|}{b_1}-\right. </math><math>\left. \frac{\left|y_1-y_2\right|}{b_2}\right)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 40)
|}

donde <math display="inline">C_0</math>  es el valor de la varianza [[#bib0260|[52]]]  y <math display="inline">b_1</math>  y <math display="inline">b_2</math>  son longitudes de correlación del campo aleatorio [[#bib0150|[30]]]  and [[#bib0260|[52]]] . La ecuación (40) es válida en un dominio cuadrado <math display="inline">\Omega =</math><math>\left(-a,a\right)\times \left(-a,a\right)</math> . Siguiendo el procedimiento dado en [[#bib0150|[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):

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 41)
|}

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

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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. </math><math>\left. \left. \frac{\left|y_1-y_2\right|}{b_2}\right)k_{j_2}(y_2)dy_2\right)=</math><math>\left({\lambda }_{j_1}k_{j_1}(x_1)\right)\left({\lambda }_{j_2}k_{j_2}(y_2)\right)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 42)
|}

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

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 43)
|}

donde <math display="inline">-a\leq x_2\leq a</math>  y <math display="inline">-</math><math>a\leq y_2\leq a</math> . 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)—:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\lambda }_{j_1}k_{j_1}(x_1)=\sqrt{C_0}{\int }_{-a}^{x_1}exp\left(-\right. </math><math>\left. \frac{\left(x_1-x_2\right)}{b_1}\right)k_{j_1}(x_2)dx_2+</math><math>\sqrt{C_0}{\int }_{x_1}^aexp\left(\frac{\left(x_1-x_2\right)}{b_1}\right)k_{j_1}(x_2)dx_2</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 44)
|}

Diferenciando con respecto a <math display="inline">x_1</math> , se tiene (45):

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\lambda }_{j_1}\frac{dk_{j_1}(x_1)}{dx_1}=-\sqrt{C_0}\frac{1}{b_1}exp\left(-\right. </math><math>\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+</math><math>\sqrt{C_0}\frac{1}{b_1}exp\left(\frac{\left(x_1\right)}{b_1}\right){\int }_{x_1}^aexp\left(-\right. </math><math>\left. \frac{\left(x_2\right)}{b_1}\right)k_{j_1}(x_2)dx_2</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 45)
|}

Diferenciando, nuevamente, con respecto a <math display="inline">x_1</math> , se tiene (46):

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\lambda }_{j_1}\frac{d^2k_{j_1}(x_1)}{d{x_1}^2}=\left(h^2{\lambda }_{j_1}-\right. </math><math>\left. 2\sqrt{C_0}h\right)k_{j_1}(x_1)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 46)
|}

Definiendo <math display="inline">{\alpha }^2=\left(2\sqrt{C_0}h-\right. </math><math>\left. h^2{\lambda }_{j_1}\right)/{\lambda }_{j_1}</math> , la ecuación (46) se convierte en (47):

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\frac{d^2k_{j_1}(x_1)}{d{x_1}^2}+{\alpha }^2k_{j_1}(x_1)=</math><math>0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 47)
|}

donde <math display="inline">-a\leq x_1\leq a</math> . La solución de (47) es de la forma (48):

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>k_j(x_1)=c_1e^{i\alpha x_1}+c_2e^{-i\alpha x_1},{\alpha }^2\geq 0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 48)
|}

donde <math display="inline">c_1</math>  y <math display="inline">c_2</math>  son constantes. Aplicando las condiciones de contorno se obtiene la solución dada por (49) [[#bib0280|[56]]] :

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\lambda }_{j_1}=\frac{2\sqrt{C_0}h}{h^2+{{\alpha }_{j_1}}^2}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 49)
|}

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

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>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}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 50)
|}

donde se tiene el arreglo <math display="inline">{\alpha }_1<{\alpha }_2<...<{\alpha }_n</math> , que permite encontrar los autovalores ordenados mediante <math display="inline">{\lambda }_1>{\lambda }_2>...>{\lambda }_n</math> .      

===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 ([[#fig0010|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.

<span id='fig0010'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_García_865974548-1-s2.0-S0213131513000357-gr2.jpg|center|372px|Malla utilizada en el análisis por elementos finitos estocásticos espectrales.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 2.

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

</span>
|}

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 '''φ'''<sub>'''0'''</sub>[[#bib0150|[30]]] . En los ejemplos aquí desarrollados se utilizan unas condiciones iniciales según el tipo de ejemplo ('''φ'''<sub>'''0'''</sub> ). 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''<sup>''*''</sup>  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. [[#bib0140|[28]]] , está definido por la ecuación diferencial (51), para un dominio <math display="inline">\Omega =</math><math>\left(-0,5,0,5\right)\times \left(-0,5,0,5\right)</math> .

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\frac{\partial \varphi }{\partial t}-k\Delta \varphi +</math><math>\boldsymbol{\mbox{u}}\cdot \nabla \varphi =0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 51)
|}

con <math display="inline">k=0,0001</math>  y un campo de velocidad rotacional <math display="inline">\boldsymbol{u}=</math><math>\left[-4y,4x\right]</math> . 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 <math display="inline">\sigma =</math><math>0,0477</math> . La solución analítica del problema se plantea en la expresión (52).

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\varphi (x,y,t)=\frac{2{\sigma }^2}{2{\sigma }^2+4kt}e^{-\lambda }</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 52)
|}

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

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\varphi (x,y,0)=e^{-\frac{{\left(x+0,25\right)}^2+y^2}{2{\sigma }^2}}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 53)
|}

en esta última expresión:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | 
|}

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 [[#fig0015|figura 3]]  (a-e). Se empleó una discretización temporal con pasos de tiempo <math display="inline">\Delta t=</math><math>0,0005</math>  y un tiempo final <math display="inline">t_f=</math><math>0,5</math> . Se muestran los resultados para la malla de 2.500 ([[#fig0015|figs. 3]] a y 3d), 625 ([[#fig0015|figs. 3]] b y 3e) y 144 elementos ([[#fig0015|figs. 3]] c y 3f), respectivamente. Además, la [[#fig0015|figura 3]]  muestra la solución Petrov-Galerkin ([[#fig0015|figs. 3]] a-3c) y Bunbov-Galerkin ([[#fig0015|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).

<span id='fig0015'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_García_865974548-1-s2.0-S0213131513000357-gr3.jpg|center|320px|Resultados para el problema determinista con el uso de Petrov-Galerkin para a. ...]]


|-
| <span style="text-align: center; font-size: 75%;">

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.

</span>
|}

En la [[#fig0015|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 ( [[#fig0015|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 ( [[#fig0015|figs. 3]] a-3c).      

Para el caso estocástico se utilizaron los siguientes datos: <math display="inline">C_0=</math><math>1.e-4</math> , <math display="inline">\left|a\right|=</math><math>0,5</math> , <math display="inline">b_1=b_2=0,25</math>  (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 <math display="inline">{\alpha }_1=</math><math>2,1537,{\alpha }_2=4,5779,{\alpha }_3=7,2868,{\alpha }_4=</math><math>10,1732\mbox{.}</math>  En la [[#fig0020|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 ''ϕ''<sub>0</sub> , ''ϕ''<sub>1</sub> , ''ϕ''<sub>2</sub> …''ϕ''<sub>4</sub> , 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.

<span id='fig0020'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_García_865974548-1-s2.0-S0213131513000357-gr4.jpg|center|367px|Patrones para la solución del problema estocástico de la ecuación (51) con los ...]]


|-
| <span style="text-align: center; font-size: 75%;">

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 <math display="inline">{\varphi }_0</math> , <math display="inline">{\varphi }_1</math> , <math display="inline">{\varphi }_2</math>  … <math display="inline">{\varphi }_4</math> , 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).                  

</span>
|}

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 [[#bib0150|[30]]] . Según se establece en Ghanem y Spanos [[#bib0150|[30]]]  y en Choi et al. [[#bib0280|[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 [[#bib0150|[30]]] .      

Para el caso de las columnas Bubnov Galerkin (columnas del tipo ''a  '' ) las filas (<math display="inline">{\varphi }_0</math> ) 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 <math display="inline">{\varphi }_0</math>  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 [[#bib0070|[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 [[#fig0020|figura 4]] , para <math display="inline">{\varphi }_1</math> , muestra el primer coeficiente de la expansión del polinomio de caos [[#bib0150|[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 <math display="inline">{\varphi }_2</math>  y <math display="inline">{\varphi }_3</math>  muestran resultados similares que para <math display="inline">{\varphi }_1</math> . 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 <math display="inline">{\varphi }_4</math>  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 <math display="inline">\overset{\rightarrow}{u}=</math><math>\left[-4y,4x\right]</math> ). 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 [[#bib0145|[29]]] , modela la propagación de un frente de onda dentro de un dominio cuadrado <math display="inline">\Omega =</math><math>\left(0,1\right)\times \left(0,1\right)</math> . El problema está definido por medio de la ecuación (51), para una constante <math display="inline">k=</math><math>0,01</math>  y un campo de velocidad definido por el vector <math display="inline">\boldsymbol{u}=</math><math>\left[\omega (x,t),\omega (y,t)\right]</math> , en donde:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\omega (r,t)=\frac{0,1A+0,5B+C}{A+B+C},\mbox{ }\mbox{con}\mbox{ }r=</math><math>x,y</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 54)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>A(n,t)=e^{\frac{-0,05\left(n-0,5+4,95(t+0,3)\right)}{k}}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 55)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>B(n,t)=e^{\frac{-0,25\left(n-0,5+0,75(t+0,3)\right)}{k}}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 56)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>C(n,t)=e^{\frac{-0,5\left(n-0,375\right)}{k}}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 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).

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\varphi (x,y,0)=\omega (x,0).\omega (y,0)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 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 [[#fig0025|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.

<span id='fig0025'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_García_865974548-1-s2.0-S0213131513000357-gr5.jpg|center|103px|Distribución de las zonas de alta advección y difusión para el caso 2.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 5.

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

</span>
|}

Los resultados deterministas que se obtuvieron en este ejemplo se muestran en la [[#fig0030|figura 6]] , donde se empleó una discretización temporal con pasos de tiempo <math display="inline">\Delta t=</math><math>0,0005</math>  y un tiempo final <math display="inline">t_f=</math><math>0,5</math> .

<span id='fig0030'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_García_865974548-1-s2.0-S0213131513000357-gr6.jpg|center|320px|Resultados para el problema determinista del segundo caso con el uso de ...]]


|-
| <span style="text-align: center; font-size: 75%;">

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.

</span>
|}

Se observa ([[#fig0030|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 [[#fig0025|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: <math display="inline">C_0=</math><math>1.e-4</math> , <math display="inline">\left|a\right|=</math><math>0.5</math> , <math display="inline">b_1=b_2=0,25</math>  (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 <math display="inline">{\alpha }_1=</math><math>2,1537,{\alpha }_2=4,5779,{\alpha }_3=7,2868,{\alpha }_4=</math><math>10,1732</math> . En la [[#fig0035|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 <math display="inline">{\varphi }_0</math> , <math display="inline">{\varphi }_1</math> , <math display="inline">{\varphi }_2</math>  … <math display="inline">{\varphi }_4</math> , 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.

<span id='fig0035'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_García_865974548-1-s2.0-S0213131513000357-gr7.jpg|center|360px|Patrones para la solución del problema estocástico de la ecuación (51) con los ...]]


|-
| <span style="text-align: center; font-size: 75%;">

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 <math display="inline">{\varphi }_0</math> , <math display="inline">{\varphi }_1</math> , <math display="inline">{\varphi }_2</math>  … <math display="inline">{\varphi }_4</math> , 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).                  

</span>
|}

La fila ''a''  de la  [[#fig0035|figura 7]]  muestra los resultados para <math display="inline">{\varphi }_0</math> . 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: <math display="inline">{\varphi }_1</math> . 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 <math display="inline">{\varphi }_2</math> , 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 <math display="inline">{\varphi }_2</math>  cerca y alrededor de la interfaz de la onda. Si se compara con el caso para <math display="inline">{\varphi }_1</math>  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 <math display="inline">{\varphi }_2</math>  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 ([[#fig0025|fig. 5]] ).      

Las filas ''d''  y ''e''  de la  [[#fig0035|figura 7]]  muestran resultados para <math display="inline">{\varphi }_3</math>  y <math display="inline">{\varphi }_4</math> , 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 (<math display="inline">{\varphi }_i</math> ) 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 [[#bib0150|[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 C<sub>0</sub>  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==

<ol style='list-style-type: none;margin-left: 0px;'><li><span id='bib0005'></span>
[[#bib0005|[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</li>
<li><span id='bib0010'></span>
[[#bib0010|[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</li>
<li><span id='bib0015'></span>
[[#bib0015|[3]]] S. Ferreira, M. Martins, M. Vilela; Reaction-diffusion model for the growth of avascular tumor; Phys. Rev., 65 (2) (2002)</li>
<li><span id='bib0020'></span>
[[#bib0020|[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</li>
<li><span id='bib0025'></span>
[[#bib0025|[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</li>
<li><span id='bib0030'></span>
[[#bib0030|[6]]] S. Kondo, R. Asai; A reaction-diffusion wave on the skin of the marine anglefish, Pomacanthus; Nature, 376 (1995), pp. 765–768</li>
<li><span id='bib0035'></span>
[[#bib0035|[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</li>
<li><span id='bib0040'></span>
[[#bib0040|[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</li>
<li><span id='bib0045'></span>
[[#bib0045|[9]]] B. Rothschild, J. Ault; Population-dynamic instability as a cause of patch structure; Ecol. Model., 93 (1996), pp. 237–239</li>
<li><span id='bib0050'></span>
[[#bib0050|[10]]] T. Nozakura, S. Ikeuchi; Formation of dissipative structures in galaxies; Astrophys J, 279 (1984), pp. 40–52</li>
<li><span id='bib0055'></span>
[[#bib0055|[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</li>
<li><span id='bib0060'></span>
[[#bib0060|[12]]] O. Richter; Modelling dispersal of populations and genetic information by finite element methods; Environ. Modell. Softw., 23 (2) (2008), pp. 206–214</li>
<li><span id='bib0065'></span>
[[#bib0065|[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</li>
<li><span id='bib0070'></span>
[[#bib0070|[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</li>
<li><span id='bib0075'></span>
[[#bib0075|[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</li>
<li><span id='bib0080'></span>
[[#bib0080|[16]]] J. Chrispell, V. Ervin, E. Jenkins; A fractional step θ-method for convection–diffusion problems; J. Math. Anal. Appl., 333 (2007), pp. 204–218</li>
<li><span id='bib0085'></span>
[[#bib0085|[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</li>
<li><span id='bib0090'></span>
[[#bib0090|[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)</li>
<li><span id='bib0095'></span>
[[#bib0095|[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</li>
<li><span id='bib0100'></span>
[[#bib0100|[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</li>
<li><span id='bib0105'></span>
[[#bib0105|[21]]] T.J.R. Hughes, editor, Finite element methods for convection dominated flows, AMD 3 (1979) 19-35. ASME, New York.</li>
<li><span id='bib0110'></span>
[[#bib0110|[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</li>
<li><span id='bib0115'></span>
[[#bib0115|[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</li>
<li><span id='bib0120'></span>
[[#bib0120|[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</li>
<li><span id='bib0125'></span>
[[#bib0125|[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</li>
<li><span id='bib0130'></span>
[[#bib0130|[26]]] B. Cockburn, G. Karniadakis, C. Shu; Discontinuos Galerkin methods, theory, computational and application; Ed. Springer, Berlin (2000)</li>
<li><span id='bib0135'></span>
[[#bib0135|[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</li>
<li><span id='bib0140'></span>
[[#bib0140|[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</li>
<li><span id='bib0145'></span>
[[#bib0145|[29]]] M. Berzins; Temporal error control for convection-dominated equations in two space dimensions; SIAM J Sci Comput, 16–3 (1995), pp. 558–580</li>
<li><span id='bib0150'></span>
[[#bib0150|[30]]] R. Ghanem, P. Spanos; Stochastic finite elements. A spectral approach; Dover Publications., Minneola, New York, USA (1991)</li>
<li><span id='bib0155'></span>
[[#bib0155|[31]]] E.W. Boyce, B.E. Goodwin; Random transverse vibration of elastic beams; SIAM J, 12–3 (1964), pp. 613–629</li>
<li><span id='bib0160'></span>
[[#bib0160|[32]]] J.D. Collins, W.T. Thompson; The eigenvalue problem for structural systems with uncertain parameters; AIAA J, 7–4 (1969), pp. 642–648</li>
<li><span id='bib0165'></span>
[[#bib0165|[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</li>
<li><span id='bib0170'></span>
[[#bib0170|[34]]] T.K. Hasselman, G.C. Hart; Modal analysis of random structural systems; J. Eng. Mech., 98–EM3 (1972), pp. 561–579</li>
<li><span id='bib0175'></span>
[[#bib0175|[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</li>
<li><span id='bib0180'></span>
[[#bib0180|[36]]] W.K. Liu, G. Besterfield, T. Belytschko; Transient probabilistic systems; Comput. Methods Appl. Mech. Engrg., 67 (1988), pp. 27–54</li>
<li><span id='bib0185'></span>
[[#bib0185|[37]]] M. Kaminski; Generalized perturbation-based stochastic finite element method in elastostatics; Comput. Struct., 85 (2007), pp. 586–594</li>
<li><span id='bib0190'></span>
[[#bib0190|[38]]] S. Lepage; Stochastic finite element method for the modeling of thermoelastic damping in micro-resonators. Ph. D. Thesis; Cranfield University (2006)</li>
<li><span id='bib0195'></span>
[[#bib0195|[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</li>
<li><span id='bib0200'></span>
[[#bib0200|[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.</li>
<li><span id='bib0205'></span>
[[#bib0205|[41]]] R.G. Ghanem; Stochastic finite elements with multiple random non-Gaussian properties; J Eng Mech ASCE, 125–1 (1999), pp. 26–40</li>
<li><span id='bib0210'></span>
[[#bib0210|[42]]] D. Xiu, G. Karniadakis; The Wiener-Askey polynomial chaos for stochastic differential equations; J Sci Comput, 24–2 (2002), pp. 619–644</li>
<li><span id='bib0215'></span>
[[#bib0215|[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</li>
<li><span id='bib0220'></span>
[[#bib0220|[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)</li>
<li><span id='bib0225'></span>
[[#bib0225|[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.</li>
<li><span id='bib0230'></span>
[[#bib0230|[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</li>
<li><span id='bib0235'></span>
[[#bib0235|[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</li>
<li><span id='bib0240'></span>
[[#bib0240|[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</li>
<li><span id='bib0245'></span>
[[#bib0245|[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</li>
<li><span id='bib0250'></span>
[[#bib0250|[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</li>
<li><span id='bib0255'></span>
[[#bib0255|[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</li>
<li><span id='bib0260'></span>
[[#bib0260|[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</li>
<li><span id='bib0265'></span>
[[#bib0265|[53]]] P. Malliavin; Stochastic analysis; Springer Verlag, New York (1997)</li>
<li><span id='bib0270'></span>
[[#bib0270|[54]]] M. Loève; Probability theory;  (4th ed)Springer, New York (1977)</li>
<li><span id='bib0275'></span>
[[#bib0275|[55]]] R. Cameron, W. Martin; The orthogonal development of nonlinear functional in series of Fourier-Hermite functionals; Ann Math, 48 (1947), pp. 385–392</li>
<li><span id='bib0280'></span>
[[#bib0280|[56]]] S.K. Choi, R.V. Grandhi, R.A. Canfield; Reliability-based Structural design; Springer Verlag, London, UK (2007)</li>
<li><span id='bib0285'></span>
[[#bib0285|[57]]] B. Sudret, M. Berveiller; Application of a stochastic finite element procedure to reliability analysis; Proc. 11th IFIP WG7.5 Conference, Banff, Canada (2003)</li>
</ol>