Método de volumen finito para flujo no estacionario

En mecánica de fluidos , un flujo no estacionario es el transporte de un fluido cuyas propiedades dependen del tiempo . Este tipo de fenómeno pertenece al paradigma de la mecánica de medios continuos , lo que significa que las ecuaciones que lo rigen involucran diferencias infinitesimales , es decir, diferenciales de funciones en las variables espaciales y temporales continuas.

El método de volumen finito es un enfoque de dinámica de fluidos computacional mediante el cual el sistema puede resolverse numéricamente a través de la discretización de la geometría simulada , lo que permite que el cálculo paso a paso de la solución también tenga en cuenta las derivadas temporales de las propiedades del fluido . ^{[ 1 ]}

La ecuación de conservación para el transporte de un escalar en flujo no estacionario tiene la forma general ^{[ 2 ].}

${\displaystyle {\frac {\partial \rho \phi }{\partial t}}+\operatorname {div} \left(\rho \phi \upsilon \right)=\operatorname {div} \left(\Gamma \operatorname {grad} \phi \right)+S_{\phi },}$

dónde${\displaystyle \rho }$es la densidad y${\displaystyle \phi }$es la forma conservativa de todo flujo de fluidos;${\displaystyle \Gamma }$es el coeficiente de difusión y${\displaystyle S}$es el término fuente. La tasa neta de flujo de${\displaystyle \phi }$fuera del elemento fluido ( convección ) es${\displaystyle \operatorname {div} \left(\rho \phi \upsilon \right),}$y${\displaystyle \operatorname {div} \left(\Gamma \operatorname {grad} \phi \right)}$es la tasa de aumento de${\displaystyle \phi }$debido a la difusión ;${\displaystyle S_{\phi }}$es la tasa de aumento de${\displaystyle \phi }$debido a fuentes y${\displaystyle {\frac {\partial \rho \phi }{\partial t}}}$es la tasa de aumento de${\displaystyle \phi }$del elemento fluido (transitorio).

El primer término de la ecuación refleja la inestabilidad del flujo y está ausente en el caso de flujos estacionarios. La integración de volumen finito de la ecuación gobernante se realiza sobre un volumen de control y también sobre un intervalo de tiempo finito ∆t.

${\displaystyle \int \limits _{cv}\!\!\!\int _{t}^{t+\Delta t}\left({\frac {\partial \rho \phi }{\partial t}}\,\mathrm {d} t\right)\,\mathrm {d} V+\int _{t}^{t+\Delta t}\!\!\!\int \limits _{A}\left(n\cdot {\rho \phi u}\,\mathrm {d} A\right)\,\mathrm {d} t=\int _{t}^{t+\Delta t}\!\!\!\int \limits _{A}\left(n\cdot \left(\Gamma \operatorname {grad} \phi \right)\,\mathrm {d} A\right)\,\mathrm {d} t+\int _{t}^{t+\Delta t}\!\!\!\int \limits _{cv}S_{\phi }\,\mathrm {d} V\,\mathrm {d} t}$

La integración del volumen de control de la parte estacionaria de la ecuación es similar a la integración de la ecuación gobernante del estado estacionario . Debemos centrarnos en la integración del componente transitorio de la ecuación. Para comprender mejor la técnica de integración, nos referimos a la ecuación de conducción de calor transitoria unidimensional . ^{[ 3 ]}

${\displaystyle \rho c{\frac {\partial T}{\partial t}}={\frac {\partial }{\partial x}}{\frac {k\partial T}{\partial x}}+S}$

${\displaystyle \int _{t}^{t+\Delta t}\!\!\!\int \limits _{cv}\rho c{\frac {\partial T}{\partial t}}\,\mathrm {d} V\,\mathrm {d} t=\int _{t}^{t+\Delta t}\!\!\!\int \limits _{cv}{\frac {\partial }{\partial x}}{\frac {k\partial T}{\partial x}}\,\mathrm {d} V\,\mathrm {d} t+\int _{t}^{t+\Delta t}\!\!\!\int \limits _{cv}S\,\mathrm {d} V\,\mathrm {d} t}$

${\displaystyle \int _{e}^{w}\!\!\!\int _{t}^{t+\Delta t}\left(\rho c{\frac {\partial T}{\partial t}}\,\mathrm {d} t\right)\,\mathrm {d} V=\int _{t}^{t+\Delta t}\left[\left(kA{\frac {\partial T}{\partial x}}\right)_{e}-\left(kA{\frac {\partial T}{\partial x}}\right)_{w}\right]\,\mathrm {d} t+\int _{t}^{t+\Delta t}{\bar {S}}\Delta V\,\mathrm {d} t}$

Ahora, manteniendo el supuesto de que la temperatura en el nodo prevalece en todo el volumen de control, el lado izquierdo de la ecuación se puede escribir como ^{[ 4 ].}

${\displaystyle \int \limits _{cv}\!\!\!\int _{t}^{t+\Delta t}\left(\rho c{\frac {\partial T}{\partial t}}\,\mathrm {d} t\right)\,\mathrm {d} V=\rho c\left(T_{P}-{T_{P}}^{O}\right)\Delta V}$

Al utilizar un esquema de diferencias finitas hacia atrás de primer orden , podemos escribir el lado derecho de la ecuación como:

${\displaystyle \rho c\left(T_{P}-{T_{P}}^{0}\right)\Delta V=\int _{t}^{t+\Delta t}\left[\left(K_{e}A{\frac {T_{E}-T_{P}}{\delta x_{PE}}}\right)-\left(K_{w}A{\frac {T_{P}-T_{W}}{\delta x_{WP}}}\right)\right]\,\mathrm {d} t+\int _{t}^{t+\Delta t}{\bar {S}}\Delta V\,\mathrm {d} t}$

Ahora, para evaluar el lado derecho de la ecuación, utilizamos un parámetro de ponderación.${\displaystyle \theta }$entre 0 y 1, y escribimos la integral de${\displaystyle T_{P}}$

${\displaystyle I_{T}=\int _{t}^{t+\Delta t}T_{P}\,\mathrm {d} t=\left[\theta T_{P}+\left(1-\theta \right){T_{P}}^{0}\right]\Delta t}$

Ahora bien, la forma exacta de la ecuación discretizada final depende del valor de${\displaystyle \Theta }$. Como la varianza de${\displaystyle \Theta }$es${\displaystyle 0<\Theta <1}$, el esquema que se utilizará para calcular${\displaystyle T_{P}}$depende del valor de la${\displaystyle \Theta .}$De este modo

${\displaystyle \rho c{\frac {\left(T_{P}-{T_{P}}^{0}\right)}{\Delta t}}\Delta x=\theta \left[\left(K_{e}{\frac {T_{E}-T_{P}}{\delta x_{PE}}}\right)-\left(K_{w}{\frac {T_{P}-T_{W}}{\delta x_{WP}}}\right)\right]+(1-\theta )\left[\left(K_{e}{\frac {T_{E}-T_{P}}{\delta x_{PE}}}\right)-\left(K_{w}{\frac {T_{P}-T_{W}}{\delta x_{WP}}}\right)\right]+{\bar {S}}\Delta x}$

• Esquema explícito : en el esquema explícito, el término fuente se linealiza como${\displaystyle b=S_{u}+{S_{P}}{T_{P}}^{0}}$Sustituimos.${\displaystyle \theta =0}$para obtener la discretización explícita, es decir: ^{[ 5 ]}

${\displaystyle a_{P}T_{P}=a_{w}{T_{w}}^{0}+a_{e}{T_{e}}^{0}+\left[{a_{P}}^{0}-\left(a_{w}+a_{e}-S_{P}\right)\right]{T_{P}}^{0}+S_{u}}$

dónde${\displaystyle a_{P}={a_{P}}^{0}}$. Un aspecto a destacar es que el lado derecho contiene valores en el paso de tiempo anterior y, por lo tanto, el lado izquierdo se puede calcular mediante la coincidencia hacia adelante en el tiempo. El esquema se basa en la diferenciación hacia atrás y su error de truncamiento de la serie de Taylor es de primer orden con respecto al tiempo. Todos los coeficientes deben ser positivos. Para k constante y espaciado de cuadrícula uniforme,${\displaystyle \delta x_{PE}=\delta x_{WP}=\Delta x}$Esta condición puede escribirse como

${\displaystyle \rho c{\frac {\Delta x}{\Delta t}}>{\frac {2K}{\Delta x}}}$

Esta desigualdad impone una condición estricta sobre el paso de tiempo máximo que se puede utilizar y representa una limitación seria en el esquema. Resulta muy costoso mejorar la precisión espacial porque el paso de tiempo máximo posible debe reducirse como el cuadrado de${\displaystyle \Delta x}$^{[ 6 ]}

• Esquema de Crank-Nicolson  : el método de Crank-Nicolson resulta de establecer${\displaystyle \theta ={\frac {1}{2}}}$La ecuación discretizada de conducción de calor transitoria se convierte en:

${\displaystyle a_{P}T_{P}=a_{E}\left[{\frac {T_{E}+{T_{E}}^{0}}{2}}\right]+a_{W}\left[{\frac {T_{W}+{T_{W}}^{0}}{2}}\right]+\left[{a_{P}}^{0}-{\frac {a_{E}}{2}}-{\frac {a_{W}}{2}}\right]{T_{P}}^{0}+b,}$

dónde${\displaystyle a_{P}={\frac {a_{W}+a_{E}}{2}}+{a_{P}}^{0}-{\frac {S_{P}}{2}}.}$

Dado que en la ecuación hay más de un valor desconocido de T en el nuevo nivel de tiempo, el método es implícito y se deben resolver ecuaciones simultáneas para todos los puntos nodales en cada paso de tiempo. Aunque los esquemas con${\displaystyle {\frac {1}{2}}<\theta <1}$incluyendo el esquema de Crank-Nicolson, son incondicionalmente estables para todos los valores del paso de tiempo, es más importante asegurar que todos los coeficientes sean positivos para obtener resultados físicamente realistas y acotados. Este es el caso si el coeficiente de${\displaystyle {T_{P}}^{0}}$satisface la siguiente condición

${\displaystyle {a_{P}}^{0}=\left[{\frac {a_{E}+a_{W}}{2}}\right],}$

lo que conduce a

${\displaystyle \Delta t<\rho c{\frac {\Delta x^{2}}{K}},}$

Por lo tanto, el método de Crank-Nicolson se basa en diferencias centrales y, por consiguiente, tiene una precisión de segundo orden en el tiempo. La precisión general de un cálculo también depende de la práctica de diferencias espaciales, por lo que el método de Crank-Nicolson se suele utilizar junto con diferencias centrales espaciales.

3. Esquema totalmente implícito : cuando el valor de Ѳ se establece en 1, obtenemos el esquema totalmente implícito. La ecuación discretizada es: ^{[ 7 ]}

${\displaystyle a_{P}T_{P}=a_{W}T_{W}+a_{E}T_{E}+{a_{P}}^{0}{T_{P}}^{0}+S_{u}}$

${\displaystyle a_{P}={a_{P}}^{0}+a_{W}+a_{E}-S_{P}}$

Ambos lados de la ecuación contienen temperaturas en el nuevo paso de tiempo, y se debe resolver un sistema de ecuaciones algebraicas en cada nivel de tiempo. El procedimiento de avance temporal comienza con un campo inicial dado de temperaturas. ${\displaystyle T^{0}}$El sistema de ecuaciones se resuelve después de seleccionar el paso de tiempo.${\displaystyle \Delta t}$A continuación, la solución.${\displaystyle T}$se le asigna a${\displaystyle T^{0}}$El procedimiento se repite para avanzar la solución con un paso de tiempo adicional. Se observa que todos los coeficientes son positivos, lo que hace que el esquema implícito sea incondicionalmente estable para cualquier tamaño de paso de tiempo. Dado que la precisión del esquema es solo de primer orden en el tiempo, se necesitan pasos de tiempo pequeños para garantizar la precisión de los resultados. El método implícito se recomienda para cálculos transitorios de propósito general debido a su robustez y estabilidad incondicional.