Ecuaciones de Navier-Stokes en flujos incompresibles bidimensionales

Se van a discutir en esta sección algunos aspectos de la aplicación del MEF a un problema descrito por las ecuaciones estacionarias de Navier-Stokes para flujos incompresibles de fluidos con viscosidad constante. Como ya se comentó anteriormente, en este tipo de flujo las ecuaciones de conservación de la masa y de la cantidad de movimiento están desacopladas de la ecuación de conservación de la energía, por lo que en lo que sigue se considerarán sólo las dos primeras. La ecuación de conservación de la energía podría resolverse a continuación, una vez conocida la distribución de velocidad en el flujo. Las ecuaciones a resolver son, por tanto, las ecuaciones (5.38) a (5.39):

$$\vec{\nabla}\cdot\vec{v^{*}}=0$$ (6.68)

$$\vec{v^{*}}\cdot\vec{\nabla}\vec{v^{*}}+ \vec{\nabla}p^{*} = \nu \nabla^{2}\vec{v^{*}}$$ (6.69)

($p$ se supondrá que es una presión reducida).

Utilizando una función de ponderación $s$ asociada a la ecuación (6.68) y otra función vectorial $\vec{W}$ asociada a la ecuación (6.69), la formulación de residuos ponderados puede escribirse de la forma

$$\int_{\Omega}\vec{W}\cdot (\vec{v}\cdot\vec{\nabla}\vec{v}) \Omega + \int_{\Omega}\vec{W}\cdot\vec{\nabla}p \Omega - \int_{\Omega}s\vec{\nabla}\cdot\vec{v} \Omega = \int_{\Omega}\nu\vec{W}\cdot\nabla^{2}\vec{v} \Omega.$$ (6.70)

Integrando por partes los términos segundo y cuarto de esta ecuación, se obtiene

$$\begin{split}\int_{\Omega}\vec{W}\cdot (\vec{v}\cdot\vec{\nab... ...\int_{\Gamma}p\vec{W}\cdot\vec{n} \mbox{d}\Gamma, \end{split}$$ (6.71)

siendo $\vec{n}$ el vector unitario normal al contorno.

Los términos del segundo miembro de la ecuación (6.71) para un problema bidimensional pueden ponerse de la forma

$$\int_{\Gamma}[\nu W_{x}(\vec{\nabla}u\cdot\vec{n}) + \nu W_{y}(\vec{\nabla}v\cdot\vec{n})-p(W_{x}n_{x}+W_{y}n_{y})] \Gamma.$$ (6.72)

Las condiciones de contorno para un problema correspondiente a un flujo interno en un conducto de forma arbitraria, por ejemplo, para el que se utiliza un dominio de cálculo limitado por las paredes del conducto y por las secciones de entrada y salida, serán las siguientes:

• En las paredes se impone la condición de no deslizamiento:

  $$u=0,\; v=0.$$ (6.73)

  
  

• En la sección de entrada se impone una determinada distribución de velocidad. Suponiendo que el flujo entra sólo con componente $u$, las condiciones son:

  $$u=u_{0},\; v=0.$$ (6.74)

  
  
  Tanto estas condiciones como las anteriores son obviamente condiciones de contorno esenciales.
• En la sección de salida (suponiendo que ésta es perpendicular al eje $x$ y que, por tanto, $n_{x}=1$ y $n_{y}=0$), los términos del segundo miembro de la ecuación (6.71) toman la forma

  $$\int_{\Gamma}\left[-pW_{x}+\nu\left(W_{x} \frac{\partial u}{\partial x}+W_{y} \frac{\partial v}{\partial x}\right)\right] \Gamma.$$ (6.75)

  
  
  Si se comparan estos términos con los correspondientes de la ecuación (6.37), se observa que como condiciones de contorno naturales pueden imponerse valores de las magnitudes $-p+\nu \partial u/\partial x$ y $\nu \partial v/\partial x$. Si se los gradientes de $u$ y $v$ según $x$ pueden suponerse despreciables en la sección de salida, lo que suele ser aceptable, y se tiene en cuenta que la presión tan sólo aparece en las ecuaciones en el término $\vec{\nabla}p$, de forma que la solución para la distribución de presiones satisface las ecuaciones si se le suma una constante, se deduce que pueden imponerse valores nulos para las magnitudes mencionadas. De esta forma, los términos del segundo miembro de la ecuación (6.71) pueden despreciarse.

Dado que la presión no interviene en la ecuación (6.71) dentro de una derivada, a diferencia de lo que ocurre con las componentes de la velocidad, cuyas derivadas primeras sí aparecen, parece razonable elegir diferentes funciones de forma para la presión y las componentes de la velocidad. Introduciendo las funciones de forma en la ecuación (6.71) se obtiene el sistema de ecuaciones

$$\mathbf{\underline{K}}(\vec{U})\vec{U}+\mathbf{\underline{C}}\vec{P}=\vec{f}$$ (6.76)

$$\mathbf{\underline{C}}^{T}\vec{U}=\vec{g},$$ (6.77)

donde $\vec{U}$ es el vector de componentes de la velocidad en los nodos, $\vec{P}$ el vector de presiones en los nodos, $\vec{f}$ y $\vec{g}$ vectores que incorporan las condiciones de contorno, $\mathbf{\underline{K}}(\vec{U})$ una matriz simétrica definida positiva, y $\mathbf{\underline{C}}$ una matriz asimétrica indefinida. Debido al carácter indefinido de la matriz $\mathbf{\underline{C}}$, es posible, por ejemplo, que para $u=v=0$, $\vec{f}=0$ y $\vec{g}=0$, la ecuación $\mathbf{\underline{C}}\vec{P}=0$ permita soluciones no triviales para $\vec{P}$. Para evitar problemas de este tipo se recurre a utilizar una interpolación mixta, empleando esquemas de interpolación diferentes para las componentes de la velocidad y de la presión. Al ser inferiores las condiciones de regularidad para la presión que para la velocidad, la presión se aproxima mediante polinomios de interpolación de un grado más bajo que los utilizados para la velocidad.

En cuanto al tipo de elementos empleados para este problema, suelen preferirse elementos triangulares con siete grados de libertad para las componentes de la velocidad y tres grados de libertad para la presión, o bien cuadriláteros con nueve y tres grados de libertad para componentes de la velocidad y presión respectivamente. Puede consultarse el texto de Cuvelier et al. (1986) indicado en la bibliografía donde se discute la elección de los elementos.

Por último, cabe indicar que la resolución de los sistemas (6.76) y (6.77) no es sencilla al contener la matriz correspondiente ceros en la diagonal. Para facilitar la resolución se suele hacer uso del método de penalización, ya mencionado en otras asignaturas. Una revisión sobre el uso del método de los elementos finitos en flujos incompresibles viscosos puede encontrase en la publicación de Gresho (1989) indicada en la bibliografía.