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.
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.
Muchos problemas de la física, la química, la economía, la biología, la bioingeniería, e incluso la ecología, entre otros campos, pueden ser modelados a través del balance de 3 fenómenos: la difusión, la convección y la reacción [1–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).
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), f(x) es la función de generación, g(x) es la función que define el valor del campo escalar φ sobre la frontera Γφ y h(x) es la función que define el valor del flujo sobre la frontera Γ∇, cuya normal es n.
El comportamiento y los patrones de la solución de los ejemplos señalados [1–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 s. En casos más generales, en ecuaciones de reacción-difusión-convección, donde s puede depender de más de una variable, del espacio y del tiempo, se puede obtener un conjunto de relaciones entre los parámetros que definen el comportamiento de la solución, por ejemplo, se pueden obtener patrones de Turing o de Hopf [2,5]. Por tanto, los parámetros definen el tipo de solución de la ecuacion y, además, definen el tipo de método numérico que se debe utilizar. En efecto, aunque la formulación convencional de elementos finitos, o método de Bubnov-Galerkin, resulta útil y adecuada para el tratamiento de muchos problemas de la ingeniería, especialmente en el campo de la mecánica de sólidos, presenta problemas de estabilidad cuando en la ecuación diferencial aparecen operadores no autoadjuntos, tal como el término convectivo en la ecuación (1). Este término, bajo la formulación convencional (Bubnov-Galerkin), tiene un efecto desestabilizador de la matriz de rigidez, ya que introduce asimetría en esta y produce oscilaciones falsas en la aproximación. La eliminación de estas oscilaciones, o estabilización de la solución, se logra regresando el carácter simétrico a la matriz de rigidez, lo cual se puede alcanzar mediante diversos métodos, entre los cuales cabe citar: el método de las líneas de características [14], el método de cálculo finito [15], el método de paso fraccional θ[16], el método de mínimos cuadrados de Galerkin [17] y el método Petrov-Galerkin de contracorriente [18–29].
De otro lado, los parámetros mencionados pueden ser deterministas o estocásticos. Cuando son deterministas el método de Petrov-Galerkin es adecuado para encontrar la solución aproximada [18]. En el caso de que los parámetros sean estocásticos (variables de campo aleatorio), no existen, hasta donde los autores conocen, reportes sobre la solución de ecuaciones de difusión-convección (con convección dominante) que permitan incluir la variabilidad de campo en los parámetros de la ecuación (1). Para solucionar el problema estocástico se han utilizado varios métodos, entre los que se cuentan el método de la perturbación y los métodos de elementos finitos estocásticos espectrales [30].
El método de la perturbación es el utilizado con más frecuencia debido a su simplicidad matemática [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 [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 [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 [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 [33] y Hasselman y Hart [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. [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. [36] se amplía el método para problemas transientes de mecánica de sólidos. Recientemente, Kaminski [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 [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 [39], en medios porosos [40] y en diseño estructural [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 [30]. El método fue originalmente desarrollado para problemas de elasticidad lineal en los que las propiedades del material presentaban campos aleatorios gaussianos. Luego, utilizando el mismo principio se extendió a problemas de transferencia de calor [41,42] y campos lognormales [41]. Desarrollos recientes han mostrado la versatilidad del método para diversas aplicaciones, entre las que se incluyen: mecánica de fluidos [43], análisis modal [44], interacción fluido estructura [45], análisis de placas laminadas [46] y optimización [47]. El método consiste en representar el campo aleatorio (que, en general, son las propiedades del material) en términos de la descomposición de Karhunen-Loeve. Posteriormente, la variable que se debe encontrar, que también tiene características estocásticas, se representa mediante los polinomios de caos. Ambas representaciones se discretizan y se implementan en el método de los elementos finitos, para lo cual se utiliza la proyección de Galerkin [30–47]. Este procedimiento se puede programar en un paquete estándar de elementos finitos, lo cual representa una ventaja económica en términos computacionales con respecto a otras técnicas como los métodos de Monte Carlo [30] y el método de la perturbación.
Pueden encontrarse algunas técnicas alternativas a las ya explicadas en la literatura. Por ejemplo, el método de los elementos finitos estocásticos difusos se ha utilizado para modelar problemas en mecánica de fluidos y sólidos [48] y transferencia de calor [49], entre otros. Este método introduce en las variables estocásticas el uso de conjuntos difusos que describen su comportamiento. De igual manera, existen otras tantas alternativas para la descripción estocástica, como el método de los elementos finitos extendidos estocásticos [50], el método de los elementos finitos estocásticos dinámicos de Neumann [51], o el método de elementos finitos estocásticos multiescala [52].
Todos los métodos aquí mencionados se han comparado típicamente con el método de Monte Carlo [30]. Este método requiere un exhaustivo cálculo con una gran población de valores aleatorios que simulan el espectro de eventos de un parámetro o una variable. La principal desventaja, motivo por el que se han desarrollado otros métodos, radica en el alto coste computacional, por tanto, los métodos alternos, principalmente el método de los elementos finitos estocásticos espectrales, han sido una buena alternativa.
Los métodos estocásticos han sido aplicados al diseño y al análisis en ingeniería con el objeto de representar fenómenos físicos de forma más precisa y fiable [30]. Esta metodología, en conjunto con las nuevas tecnologías computacionales, provee información de los fenómenos físicos en cuanto a los límites de desempeño real, las variaciones y sus implicaciones en la respuesta de cada fenómeno. En diseño en ingeniería, por ejemplo, el método de los elementos finitos estocásticos permite obtener modelos robustos sobre el comportamiento estructural y el desempeño de máquinas. El proceso de análisis estocástico provee ventajas sobre los métodos de diseño y análisis porque se pueden obtener valores medios, varianzas e intervalos de confianza que aportan mayor información sobre el comportamiento de una estructura, máquina o fenómeno físico. Adicionalmente, el análisis estocástico puede ayudar a desarrollar una guía inicial para identificar aquellas zonas en las cuales se debe hacer especial énfasis en el estudio de un fenómeno físico [53].
Desde esta perspectiva, el presente artículo desarrolla numéricamente 2 problemas de difusión-convección (DC), empleando el método de Petrov-Galerkin en contracorriente (SUPG). Los parámetros de la ecuación se modelan como variables de campo aleatorio. Así pues, se combina el método SUPG y el método de los elementos finitos estocásticos espectrales. Los objetivos del artículo son: por un lado, estudiar la influencia de los campos estocásticos en la solución de problemas de DC con SUPG y, por otro, obtener los patrones de cada coeficiente del polinomio de caos. Aunque la propuesta teórica se ha desarrollado incluyendo el término reactivo y de fuente, los ejemplos aquí solucionados no consideran estos términos. Sin embargo, los resultados muestran que la incógnita se ve afectada por los campos estocásticos y que el patrón obtenido en los términos de la expansión con el polinomio de caos muestra la incertidumbre en la incógnita, que depende de la dinámica de la difusión, la velocidad convectiva y el tipo de solución utilizado.
2Materiales y métodos2.1El método de estabilización de Petrov-Galerkin en contracorrienteEl SUPG se basa en la estabilización de la matriz de rigidez mediante el empleo de funciones de ponderación modificadas, de modo que se otorgue un mayor peso a la información de los nodos ubicados aguas arriba. Este método de estabilización fue planteado por primera vez por Hughes y Brooks [21].
En problemas unidimensionales la estabilización con el método de Petov-Galerkin se puede alcanzar añadiendo un término perturbador a la función de peso estándar (Wi), tal como se muestra en (5) [6], donde h es el tamaño característico del elemento y α es un parámetro de perturbación positivo calculado por medio de la expresión (6) [7]:
donde:
En esta última expresión Pe es el número adimensional de Peclet definido como Pe=uh2k. Nótese que la función de peso solo debe ser perturbada si Pe>1. En otras palabras, para valores inferiores a la unidad del número de Peclet, la formulación SUPG lleva implícito el método de Bubnov-Galerkin.
Para una ecuación de difusión-convección, los términos de la matriz de rigidez [K] de un elemento e, calculados de acuerdo con el método SUPG unidimensional, están definidos por medio de la expresión (7) [14]:
En esta última expresión, Nm son las funciones de forma empleadas para aproximar la solución [14], tal como se observa en (8):
De forma análoga, los términos del vector [F] para un elemento e se definen de acuerdo con (9) [14]:
Adicionalmente, este método resulta especialmente útil para problemas en 2 o 3 dimensiones. Por tanto, de forma análoga, se puede perturbar la función de ponderación para estabilizar la convección para el caso bidimensional, como se observa en (10) [21]:
donde ux y uy son las componentes de la velocidad en las direcciones globales x−y, en tanto que α es el coeficiente de perturbación calculado de acuerdo con (6), y h es la longitud característica del elemento, como se muestra en la figura 1.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):
donde l=1,...,nnod, n→ es el vector normal al borde de flujo o borde de Neumman y nnod es el número de nodos del dominio discretizado. Empleando ahora aproximaciones por tramos del tipo φe=∑m=1nnodNmφm, así como la función de peso modificada (10), se puede conducir de la ecuación global a una expresión a nivel elemental [14], como la mostrada en (12):donde:donde se ha utilizado una integración temporal del tipo Backward-Euler [15]. Además, a nivel elemental se tienen los siguientes arreglos: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 espectralesEn el caso de que los parámetros de la ecuación se puedan representar como campos aleatorios, la ecuación (11) se convierte en:
donde los parámetros del modelo (1, k, u→ y F) dependen del espacio y poseen naturaleza aleatoria, la cual está indicada por el argumento ω[41]. A continuación, se muestra la metodología que se debe seguir para resolver el problema con elementos finitos estocásticos.2.4. Expansión de Karhunen-Loève
La expansión de Karhunen-Loève [54] de un proceso estocástico, por ejemplo k(x,ω), se basa en la expansión espectral de su función de covarianza REE(x1,x2), donde x1 y x2 denotan las coordenadas espaciales y ω denota la aleatoriedad. Por definición, la función de covarianza es simétrica y definida positiva [30,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 k(x,ω) toman la forma (14):
donde k¯(x) es el valor promedio del proceso estocástico, en este caso, del coeficiente de difusión, los términos {ξj(ω)} forman un conjunto ortogonal de variables aleatorias, λj es el conjunto de autovalores y kj(x) es el conjunto de autofunciones de la función de covarianza, los cuales pueden ser evaluados como la solución de la siguiente integral:donde Ω denota el dominio espacial sobre el cual se presenta el proceso k(x,ω). El punto más importante en la formulación radica en que las fluctuaciones aleatorias espaciales se pueden descomponer en un conjunto de funciones deterministas en el espacio que están multiplicadas por los coeficientes aleatorios que, a su vez, son independientes de dichas funciones [46]. En este punto se deben mencionar 2 características importantes. La primera se refiere al tipo de proceso: si el proceso es gaussiano, las variables aleatorias {ξj(ω)} forman un vector gaussiano ortonormal. En segundo lugar, si el proceso tiene un alto ruido se deben incluir un gran número de términos en la sumatoria. Por el contrario, se puede incluir únicamente un término para describir el proceso. En general, la mayoría de los fenómenos físicos tienen materiales con propiedades que varían suavemente en las escalas de interés en ingeniería y, por tanto, se requerirán pocos términos en la expansión de Karhunen Loève [30].
En el problema aquí desarrollado las variables estocásticas (1, k, u→ y F) dadas en (13) se pueden expresar según (16), (17) y (18). En el caso de (16) se hace una expansión del coeficiente de la matriz de masa M —ver ecuación (13)— de forma similar al desarrollo hecho en Ghanem [41]:
El coeficiente 1(x,t) tiene en cuenta la variabilidad del término temporal de la ecuación. Este coeficiente es homólogo al valor de k(x,t) que se encuentra multiplicando al laplaciano de la concentración de la variable. Un acercamiento similar se ha hecho en Ghanem y Spanos [30], donde el parámetro asociado a la derivada temporal se divide en toda la ecuación donde se pierde de vista su significado físico que, en ecuaciones diferenciales parabólicas del tipo de transferencia de calor, tiene en cuenta el calor específico de la sustancia en estudio. Por tanto, el coeficiente 1(x,t) permite estudiar la variación de este coeficiente, aunque su valor promedio se encuentre en el término k(x,t).
Para el caso del coeficiente de difusión se tiene (17):
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):
Las funciones 1j1(x), kjk(x), vxjx y vyjy(x) corresponden a las autofunciones de la función de covarianza de la variable de campo. Igualmente, los valores λj11, λjkk, λjxvx y λjyvy son los autovalores. De esta forma, cada conjunto de autofunciones y autovalores conforman el denominado autopar [30].
Reemplazando (16), (17) y (18) en (13) se tiene (19):
Al definir ξ0=1, λ01=1, λ0k=1, λ0vx=1, λ0vy=1, 10(x)=1, k0(x)=k¯(x), vx0(x)=vx¯(x) y vy0(x)=vy¯(x) y definir vectorialmente la velocidad (u(x)=vx(x),vy(x)T) se llega a (20):
La ecuación (20) se puede escribir como (21):
donde:
Utilizando el método de Backward Euler para la integración temporal, la ecuación (21) se convierte en (22):
Además, mediante manipulación algebraica se llega a (23).
2.5. Expansión de los polinomios de caos
La función de covarianza de la solución, en general, no es conocida de antemano [44]. Por tanto, no se puede llevar a cabo una expansión de Karhunen Loève para describir la incógnita que hay que encontrar. En este caso la solución se puede representar por una serie de términos en una combinación no lineal que es función de las variables aleatorias {ξj(ω)}. Se ha demostrado [55] que la dependencia funcional de las variables aleatorias puede ser expresada en términos de polinomios de caos, descritos como (24):
en donde Γn(ξi1(ω),ξi2(ω),...,ξin(ω)) son los polinomios de caos de orden n en las variables {ξj(ω)}. Estos polinomios son la generalización de los polinomios multidimensionales de Hermite [30]. La validez de la convergencia de esta expansión es independiente del comportamiento del material. En general, cualquier variable aleatoria, con comportamiento estocástico desconocido, puede ser expresada como un polinomio en variables gaussianas según (24) [30]. Si se desarrollan los polinomios de caos, se puede obtener un conjunto de polinomios bien ordenados indicialmente descritos por Ψi(ω)[30,41], de tal forma que (24) se puede reescribir como (25):donde la variable φ(ω) representa un valor en «un» punto nodal. Estos polinomios (Ψi(ω)) son mutuamente ortogonales en el sentido del producto interno Ψi(ω),Ψj(ω), el cual está dado por (26):donde puede apreciarse que el producto interno está definido como el promedio (esperanza estadística) del producto de los 2 polinomios [41]. Este producto interno será cero en el caso en que i≠j. La caracterización completa del proceso aleatorio tiene lugar cuando se encuentran los coeficientes deterministas, dados por φi[30].
Por tanto, reemplazando (25) en (23) se llega a (27):
donde el vector ϕi=[φi1,φi2,....,φinnodel] representa la i-ésima incógnita estocástica.
Con un poco de álgebra se puede definir el error r dado por (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):
donde M1, MK y Mu determinan el truncamiento de la expansión de Karhunen-Loève para cada uno de los campos aleatorios. De igual forma, P determina el truncamiento de la expansión de los polinomios de caos [30]. Nótese que debido al truncamiento el error no es exactamente cero, y se convierte en un valor apreciable. Para reducir el error se hace una proyección sobre los polinomios de caos, de tal forma que el producto interno queda definido por [30]:
Para l=0,1,2,...,P. Mediante manipulación algebraica se obtiene la ecuación (31):
De Ghanem y Spanos [30] y de Choi et al. [56] se pueden obtener los productos internos indicados en (31). Por ejemplo, en el segundo término de (31) se obtiene (32) [57]:
donde el término djkil está dado por (33) [56,57]:
De igual manera se obtienen los demás coeficientes. Por tanto, el sistema toma la forma (34):
donde l=0,1,2,...,P.
Expandiendo (34) y considerando que el término F(x,ω) es determinista (esto es: F(x,ω)Ψl(ω)=F0), la ecuación se convierte en el sistema (35):
Por tanto, para la s-ésima fila se puede definir (36):
De tal manera que (35) se puede escribir de la forma (37):
Y en términos matriciales (38):
2.7. Sobre los autovalores y las autofunciones
De forma general, en 2 dimensiones, la ecuación (15) se puede escribir como:
donde REE(x,y) es la función de covarianza. En este trabajo se supone (40) [57]:donde C0 es el valor de la varianza [52] y b1 y b2 son longitudes de correlación del campo aleatorio [30,52]. La ecuación (40) es válida en un dominio cuadrado Ω=−a,a×−a,a. Siguiendo el procedimiento dado en [30] se puede llevar a cabo una separación de variables, de tal manera que la integral (39) se puede calcular separando términos según cada coordenada, por lo que se obtienen los autovalores y las autofunciones que se muestran en (41):
Por tanto, reemplazando (40) en (39) bajo la condición (41) se obtiene (42):
donde, de acuerdo con Ghanem y Spanos [30] (42) se puede resolver de forma independiente, así:donde −a≤x2≤a y −a≤y2≤a. El valor absoluto se puede remover mediante su definición, como se muestra a continuación —se lleva a cabo la demostración sobre la primera ecuación de (43)—:
Diferenciando con respecto a x1, se tiene (45):
Diferenciando, nuevamente, con respecto a x1, se tiene (46):
Definiendo α2=2C0h−h2λj1/λj1, la ecuación (46) se convierte en (47):
donde −a≤x1≤a. La solución de (47) es de la forma (48):donde c1 y c2 son constantes. Aplicando las condiciones de contorno se obtiene la solución dada por (49) [56]:
Y se obtienen para las autofunciones (50) [56]:
donde se tiene el arreglo α1<α2<...<αn, que permite encontrar los autovalores ordenados mediante λ1>λ2>...>λn.2.8. Aspectos computacionales y condiciones de contorno e iniciales
Con el fin de evaluar la precisión del método de Petrov-Galerkin en conjunto con los elementos finitos estocásticos se analizaron 2 problemas, con advección dominante. Los problemas fueron implementados en lenguaje FORTRAN y solucionados para diferentes mallas, tal como se muestra a continuación.
Para los ejemplos que se han desarrollado en este artículo se utilizaron dominios cuadrados de lado unitario (figura 2). Además, se han utilizado mallas de 2.601 (2.500), 676 (625) y 169 (144) nodos (elementos). Los elementos utilizados son cuadriláteros con funciones de forma bilineales.
Las condiciones iniciales son impuestas a la variable determinista según el problema. Para las variables estocásticas las condiciones iniciales son nulas sin importar las condiciones iniciales de la variable φ0[30]. En los ejemplos aquí desarrollados se utilizan unas condiciones iniciales según el tipo de ejemplo (φ0). De igual forma, las condiciones de Neumann (naturales) y Dirichlet son impuestas en los nodos de la frontera sobre la ecuación determinista y, para las variables estocásticas, se colocan condiciones Neumann y Dirichlet en los mismos nodos, con valor nulo, respectivamente. De otro lado, en la función de ponderación W* se utiliza la velocidad determinista para el método SUPG.
3Experimentación numérica3.1Caso 1: movimiento convectivo de una función gaussianaEste ejemplo, desarrollado entre otros autores por Wang et al. [28], está definido por la ecuación diferencial (51), para un dominio Ω=−0,5,0,5×−0,5,0,5.
con k=0,0001 y un campo de velocidad rotacional u=−4y,4x. 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 σ=0,0477. La solución analítica del problema se plantea en la expresión (52).
Por tanto, las condiciones iniciales están dadas por (53):
en esta última expresión:
Este problema se caracteriza por tener dentro del dominio espacial tanto zonas de convección dominante, ubicadas cerca de los bordes del problema, como zonas de difusión dominante, ubicadas en la región cercana al centro del dominio. Estas zonas, definidas espacialmente, no cambian en el tiempo, dado que el coeficiente de difusión y el campo de velocidades no son función de esta variable. Por lo anterior, este caso resulta especialmente interesante para evaluar la capacidad de las técnicas numéricas para aproximar problemas con condiciones convectivas variables sobre todo el dominio.
Los resultados deterministas que se obtuvieron en este ejemplo se muestran en la figura 3 (a-e). Se empleó una discretización temporal con pasos de tiempo Δt=0,0005 y un tiempo final tf=0,5. Se muestran los resultados para la malla de 2.500 (figs. 3a y 3d), 625 (figs. 3b y 3e) y 144 elementos (figs. 3c y 3f), respectivamente. Además, la figura 3 muestra la solución Petrov-Galerkin (figs. 3a-3c) y Bunbov-Galerkin (figs. 3d-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).
En la figura 3 se puede observar que las soluciones obtenidas por los 2 métodos son similares, cualitativamente, durante el desplazamiento de la función gaussiana. Sin embargo, se presentan inestabilidades que crecen de manera abrupta al aumentar el tamaño de la malla. Cuando se utiliza Bubnov-Galerkin (figs. 3d-3f), se presentan inestabilidades que se representan por líneas de contorno. Además, las oscilaciones alcanzan valores comparables con la altura de la función gaussiana. También se presenta una dispersión considerable de la variable, tanto que la concentración de la función gaussiana baja, incluso a la mitad de su valor, comparado con el método SUPG (figs. 3a-3c).
Para el caso estocástico se utilizaron los siguientes datos: C0=1.e−4, a=0,5, b1=b2=0,25 (longitud de correlación en x e y, respectivamente). El número de términos en la expansión de Karhunen-Loève es M=5 (en todos los parámetros), y en la expansión del polinomio de caos P=5. Con los valores mencionados se encuentra que α1=2,1537,α2=4,5779,α3=7,2868,α4=10,1732. En la figura 4 se muestran los resultados para un solo instante de tiempo (t=0,500): el tiempo final de la simulación. Ordenado mediante filas se encuentran cada una de las variables estocásticas ϕ0, ϕ1, ϕ2…ϕ4, respectivamente. Además, se ordena por columnas las diferentes mallas y los 2 tipos de solución: ia y ib para dominios con 144 elementos, iia y iib para 625 elementos y iiia y iiib para 2.500 elementos. Las columnas a muestran los resultados Bubnov Galerkin y las columnas b los resultados Petrov Galerkin.
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 φ0, φ1, φ2 … φ4, respectivamente. Las columnas ia, iia, iiia corresponden a la solución con Bubnov-Galerkin y las columnas ib, iib y iiib para Petrov-Galerkin. Los patrones corresponden al tiempo final de la simulación t=0,500 (adimensional).
En este artículo se ha elegido un número de términos de la expansión M=5, de forma similar al trabajo desarrollado por Ghanem y Spanos [30]. Según se establece en Ghanem y Spanos [30] y en Choi et al. [56] los primeros 5 términos contienen la mayor cantidad de información de la variación de la solución del problema. Adicionalmente, el coste computacional de aumentar el número de términos no justifica la calidad de la información obtenida del campo estocástico [30].
Para el caso de las columnas Bubnov Galerkin (columnas del tipo a) las filas (φ0) muestran que a medida que disminuye el tamaño del elemento de 144, 625 hasta 2.500 elementos (i, ii y iii, respectivamente) la respuesta se vuelve menos oscilante y dispersiva, como se comentó en el caso determinista. Para el caso de 144 elementos (columna ia) se muestran picos de concentración de la variable φ0 cerca de la frontera inferior izquierda y superior derecha. En la columna iia 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 iiia se muestra que los resultados son similares a los hallados mediante SUPG. Las figuras iiia y iiib son comparables debido a que el número de Peclet [14] es bajo y cercano a la unidad, por lo que el efecto de la función de ponderación perturbadora —ecuación (10)— es casi despreciable.
La siguiente fila de la figura 4, para φ1, muestra el primer coeficiente de la expansión del polinomio de caos [30]. Nuevamente las columas ia, iia y iiia 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 ia, que va disminuyendo para iia y iiia. En la columna ia 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 iia, 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 iiia se observa un resultado similar al de la figura iia. Se muestra un valor altamente localizado en el máximo de la campana gaussiana. Alrededor de la campana (en las columnas iia y iiia) 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 (ib, iib y iiib) para SUPG se muestran resultados similares a los de las columnas a. Sin embargo, en ib 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 iib y iiib se muestran patrones similares a los del método de Bubnov Galerkin.
Las filas φ2 y φ3 muestran resultados similares que para φ1. Nuevamente, en el caso de Bubnov-Galerkin (columnas ia, iia y iiia) se nota la alta variabilidad en la respuesta del sistema. Además, en ia 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 φ4 se observan los patrones obtenidos para la cuarta componente del polinomio de caos. En las columnas ia y iia se muestran puntos de alta concentración en la región inferior izquierda, donde se presenta la mayor convección (ver ecuación u→=−4y,4x). En ib, iib, iiia y iiib se reconoce un valor alto en el centro de la campana gaussiana. Nuevamente, en las figuras iiia y iiib, 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.2Caso 2: movimiento convectivo de un frente de ondaEste ejemplo, desarrollado en Berzins [29], modela la propagación de un frente de onda dentro de un dominio cuadrado Ω=0,1×0,1. El problema está definido por medio de la ecuación (51), para una constante k=0,01 y un campo de velocidad definido por el vector u=ω(x,t),ω(y,t), en donde:
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).
El ejemplo resulta interesante porque el campo de velocidades asociado al término convectivo es variable en el tiempo y en el espacio, formando una zona predominantemente convectiva que crece, a la vez que se reduce la zona de difusión dominante, tal como se muestra en la figura 5. Esta condición exige una mayor capacidad de estabilización de la técnica numérica empleada, tanto en la dimensión espacial como en la temporal.
Los resultados deterministas que se obtuvieron en este ejemplo se muestran en la figura 6, donde se empleó una discretización temporal con pasos de tiempo Δt=0,0005 y un tiempo final tf=0,5.
Se observa (fig. 6) que la malla más fina, formada por 2.500 elementos, no tiene diferencias significativas en las aproximaciones alcanzadas por las 2 técnicas. En especial se debe observar la solución con Bubnov Galerkin para 144 elementos, donde se nota una oscilación cercana al límite de la zona con convección fuerte (ver la fig. 5 para la referencia). Ambas técnicas, para 625 y 2.500 elementos, logran aproximaciones estables y exactas, aun en la zona de convección dominante.
Para el caso estocástico, al igual que en el anterior ejemplo, se utilizaron los siguientes datos: C0=1.e−4, a=0.5, b1=b2=0,25 (longitud de correlación en x e y, respectivamente). El número de términos en la expansión de Karhunen-Loève es M=5, y en la expansión del polinomio de caos P=5. Con los valores mencionados se encuentra que α1=2,1537,α2=4,5779,α3=7,2868,α4=10,1732. En la figura 7 se muestran los resultados para un solo instante de tiempo (t=0,500): el tiempo final de la simulación. Los datos son ordenados de la misma forma que para el anterior caso, esto es: mediante filas se encuentran cada una de las variables estocásticas φ0, φ1, φ2 … φ4, respectivamente. Además, se ordenan por columnas las diferentes mallas y los 2 tipos de solución: ia y ib para dominios con 144 elementos, iia y iib para 625 elementos y iiia y iiib para 2.500 elementos. Las columnas a muestran los resultados Bubnov Galerkin y las columnas b los resultados Petrov Galerkin.
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 φ0, φ1, φ2 … φ4, respectivamente. Las columnas ia, iia y iiia corresponden a la solución con Bubnov-Galerkin y las columnas ib, iib y iiib para Petrov-Galerkin. Los patrones corresponden al tiempo final de la simulación t=0,500 (adimensional).
La fila a de la figura 7 muestra los resultados para φ0. En la primera columna (ia) 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 ib, iia, iib, iiia y iiib 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: φ1. La columna ia (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 ib se ha utilizado el método de Petrov Galerkin, lo que permite hallar patrones más suaves y sitios de menor incertidumbre. Las figuras ib, iia, iib, iiia y iiib 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 φ2, columna ia, se observa, nuevamente, una gran dispersión en la solución. El patrón muestra una zona extendida de alta variabilidad. La columna ib muestra la formación de altas concentraciones de φ2 cerca y alrededor de la interfaz de la onda. Si se compara con el caso para φ1 se observa una oscilación en dirección longitudinal de la interfaz. A medida que se disminuye el tamaño del elemento y se utiliza el método de Petrov-Galerkin, el patrón de respuesta de φ2 se vuelve cada vez más concentrado en la interfaz de la onda. Adicionalmente se observa un patrón con altas concentraciones en la zona de convección (fig. 5).
Las filas d y e de la figura 7 muestran resultados para φ3 y φ4, respectivamente. Una vez más se observa una gran variabilidad en la interfaz de la onda. Adicionalmente, se presentan 2 y 3 oscilaciones longitudinales en d y e, respectivamente. El patrón es de nuevo menos disperso en el caso de mallas con menor tamaño de elemento y con el uso de Petrov-Galerkin.
4Discusión y conclusionesSe 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 (φi) que exhiben una alta dispersión y variabilidad de la respuesta en todo el dominio. Por el contrario, conforme se disminuye el tamaño del elemento y se usa el método de Petrov-Galerkin, se hallan patrones con mayor concentración cerca de las zonas donde se presentan los mayores gradientes de difusión y convección. En general, aunque el método de Petrov-Galerkin y las mallas refinadas eliminen las oscilaciones falsas que se deben al término convectivo, los coeficientes estocásticos del polinomio de caos pueden captar aquellas zonas de oscilación leve.
La incertidumbre de las propiedades del material y de la velocidad convectiva puede ser incluida con funciones de covarianza que, en general, deben ser diferentes. Además, se debe tener en cuenta que la variabilidad de las propiedades puede tener diferentes escalas de incertidumbre. Es así como se puede incluir una mayor cantidad de términos en aquellos parámetros donde se requiere captar con mayor precisión la variabilidad en el espacio. En este artículo se utilizó la misma función de covarianza, funciones de correlación y valores de variación para todos los parámetros estocásticos. La utilización de un mayor número de términos en la expansión de Karhunen Loève indicaría que se presenta incertidumbre a escala micro y macroscópica. Por el contrario, una pequeña cantidad de términos indicaría que la incertidumbre se presenta a escala microscópica [30]. En estudios futuros se llevarán a cabo simulaciones que permitan modelar fenómenos descritos por ecuaciones de difusión convección, con parámetros que se encuentran en diferentes escalas.
En este trabajo se hicieron pruebas previas (no reportadas) donde se aumenta el valor del coeficiente de variación C0 y se aumenta la cantidad de términos en la expansión de Karhunen Loève y los polinomios de caos. Al incluir una gran cantidad de términos o utilizar un coeficiente de variación con alto valor, diverge la solución del problema. Este análisis permite, desde el punto de vista numérico, encontrar el límite en donde los parámetros tienen su valor máximo de variación y longitud de correlación que dan cuenta de la variabilidad. Unos valores altos de incertidumbre pueden dirigir la solución hacia patrones «caóticos», desde el punto de vista numérico.
En este trabajo se utilizó la función de ponderación perturbada de SUPG para estabilizar la solución del problema de difusión convección. Para calcular la parte de perturbación —ver ecuaciones (5) y (10)—, se propone el uso del coeficiente cero (determinista) de la velocidad del medio. De nuevo, si no se utilizan variables de campo y coeficientes estocásticos, la formulación se reduce al problema determinista.
En conclusión, el método de los elementos finitos estocásticos espectrales con aplicación a ecuación de difusión-convección permite incluir los parámetros de flujo (velocidad) y difusión dependientes del espacio y de la aleatoriedad. La dinámica de los patrones de los coeficientes del polinomio de caos está bajo la influencia del tipo de método de solución. El uso del método de Bubnov Galerkin deja al descubierto la obtención de patrones con alta dispersión de incertidumbre en todo el dominio, mientras que el método de Petrov-Galerkin permite obtener zonas localizadas. En general, la localización de la incertidumbre se observa en aquellos lugares de alto gradiente en la respuesta determinista. Estas zonas se pueden captar con mayor precisión mediante la refinación de la malla o mediante una técnica de estabilización como SUPG.