Mauricio Villa Amaya, Alejandro Sánchez Yali
Introducción
Considere la transferencia de calor en una placa rectangular
representada por el dominio $\Omega = [a, b]\times [c, d] \subset \mathbb{R}^2$, cuyas
fronteras se mantienen a temperatura constante. Asumiendo que la temperatura
$T(x, y, t)$ en un punto $(x, y)\in \Omega$ para un tiempo $t$ satisface la
siguiente ecuación:
\begin{eqnarray*} \partial_{t}T(x, y, t)-\alpha\nabla^{2}T(x, y, t)&=&0\hspace{0.5cm} \forall (x, y)\in \mathring{\Omega} \hspace{0.5cm} \forall t >0 \hspace{0.5cm} \alpha \neq 0,\\ T(x, y, t) &=& h(x, y) \hspace{0.5cm} \forall (x, y)\in \partial\Omega \hspace{0.5cm} \forall t >0,\\ T(x, y, 0)&=& T_{_0}(x, y) \hspace{0.5cm}\forall (x, y)\in \mathring{\Omega} \end{eqnarray*}
\begin{eqnarray*} \partial_{t}T(x, y, t)-\alpha\nabla^{2}T(x, y, t)&=&0\hspace{0.5cm} \forall (x, y)\in \mathring{\Omega} \hspace{0.5cm} \forall t >0 \hspace{0.5cm} \alpha \neq 0,\\ T(x, y, t) &=& h(x, y) \hspace{0.5cm} \forall (x, y)\in \partial\Omega \hspace{0.5cm} \forall t >0,\\ T(x, y, 0)&=& T_{_0}(x, y) \hspace{0.5cm}\forall (x, y)\in \mathring{\Omega} \end{eqnarray*}
Donde el laplaciano $\nabla^2$ actúa únicamente sobre las variables
$x$, $y$, la función $h(x, y)$ representa la temperatura en
la frontera que no varÃa con el tiempo y la función $T_{_0}(x, y)$ representa
la distribución inicial de temperatura sobre la placa.
Para valores pequeños de $t$ se puede modelar la variación de la temperatura en la placa,
pero dado a que el sistema tiende al equilibrio, la distribución de temperatura se mantendrá constante
para valores muy grandes de $t$, esto puede ser expresado matemáticamente con el lÃmite:
\begin{equation} T(x, y)=\lim_{t\to \infty}T(x, y, t). \end{equation}
\begin{equation} T(x, y)=\lim_{t\to \infty}T(x, y, t). \end{equation}
Por lo tanto es razonable esperar que la función de temperatura $T$ dependa únicamente de la posición,
entonces en el estado estacionario se satisface:
\begin{equation}\label{eq:problemadirichlet} \left\{\begin{array}{c} \nabla^2 T(x, y)=0\hspace{0.5cm}\forall (x, y) \in \mathring{\Omega},\\ T(x, y) = h(x, y)\hspace{0.5cm} \forall (x, y) \in \partial \Omega.\\ \end{array}\right. \end{equation}
\begin{equation}\label{eq:problemadirichlet} \left\{\begin{array}{c} \nabla^2 T(x, y)=0\hspace{0.5cm}\forall (x, y) \in \mathring{\Omega},\\ T(x, y) = h(x, y)\hspace{0.5cm} \forall (x, y) \in \partial \Omega.\\ \end{array}\right. \end{equation}
La ecuación de $\nabla^2 T(x, y)=0$ y en conjunto el problema con valores en la
frontera \eqref{eq:problemadirichlet} se conoce como el problema de Dirichlet.
El objetivo a continuación será encontrar una solución aproximada por
el método de Euler al problema \eqref{eq:problemadirichlet}.
Método de Euler
Sean $I_{m} = \{0, 1, \dots, m\}$ y $I_{n} = \{0, 1, \dots, n\}$.
Entonces la discretización del conjunto $\Omega = [a, b] \times
[c, d]$ es dada por la malla o partición $P(\Omega)=\{(x_{_j},
y_{_k}):(j,k)\in I_{_m}\times I_{_n}\}$, donde $x_{_j}= a+j\Delta x$ con
$\Delta x = (b-a)/m$ y $y_{_k}= c+k\Delta y$ con $\Delta y = (d-c)/n
$. Los números $m$ y $n$ se interpretan como el número de subintervalos en los
que se dividen $[a, b]$ y $[c, d]$ respectivamente.
Suponiendo que para el problema \eqref{eq:problemadirichlet} existe
una única solución $T(x, y)$ se emplea la serie de Taylor para aproximar las
derivadas parciales:
\begin{equation}\label{eq:segundaDerivada} \dfrac{\partial T^{2}}{\partial x^{2}} \approx \left(\dfrac{a_{_j}-2id +s_{_j}}{\Delta x^{2}}\right)(T_{_{jk}}),\hspace{0.5cm} \dfrac{\partial T^{2}}{\partial y^{2}} \approx \left(\dfrac{a_{_k}-2id +s_{_k}}{\Delta y^{2}}\right)(T_{_{jk}}) \end{equation}
donde $a_{_j}$, $s_{_j}$, $id$, $a_{_k}$, $s_{_k}$ son operadores definidos por:
\begin{equation}\label{eq:operadores} a_{_j}(T_{_{jk}})=T_{_{j-1k}},\hspace{0.5cm} s_{_j}(T_{_{jk}})=T_{_{j+1k}},\hspace{0.5cm} id(T_{_{jk}})=T_{_{jk}},\hspace{0.5cm} a_{_k}(T_{_{jk1}})=T_{_{jk-1}},\hspace{0.5cm} s_{_k}(T_{_{jk}})=T_{_{jk+1}}, \end{equation}
Hay que tener presente que $T_{_{jk}}$ representa el valor aproximado de la
temperatura para el problema \eqref{eq:problemadirichlet} en la partición
$P(\Omega)$ en el punto $(x_{_j}, y_{_k})$.
Al remplazar las expresiones \eqref{eq:segundaDerivada} en
\eqref{eq:problemadirichlet}, el problema de Dirichlet discretizado es:
\begin{equation}\label{eq:2.3}
\left\{\begin{array}{cc}
\left(\dfrac{a_{_j}-2id +s_{_j}}{\Delta x^{2}}\right)(T_{_{jk}})
+\left(\dfrac{a_{_k}-2id +s_{_k}}{\Delta y^{2}}\right)(T_{_{jk}})=0 &\hspace{0.5cm}
\forall (x_{_j},y_{_{k}}) \in P(\mathring{\Omega})\subset \mathbb{R}^2,\\
T_{_{jk}} = h_{_{jk}}&\hspace{0.5cm} \forall (x_{_j},y_{_k}) \in P(\partial \Omega).\\
\end{array}\right.
\end{equation}
Donde:
\[ P(\mathring{\Omega})=P(\Omega)\cap\mathring{\Omega} \]
\[ P(\partial \Omega) = P(\Omega)\cap \partial \Omega \]
Observe que $P(\mathring{\Omega})$ representa los puntos internos de la
malla y $P(\partial \Omega)$ los puntos de la frontera, asà que
$P(\Omega)=P(\mathring{\Omega})\cup P(\partial \Omega)$.
Al aplicar los operadores \eqref{eq:operadores} la ecuación
\eqref{eq:segundaDerivada} es equivalente a:
\begin{equation}\label{eq:sistema matricial indicial} (2\Delta y^2+2\Delta x^2)T_{_{jk}}-(\Delta y^2T_{_{j-1k}}+\Delta y^2T_{_{j+1k}}+ \Delta x^2T_{_{jk-1}}+\Delta x^2T_{_{jk+1}})=0, \hspace{0.5cm} \forall (x_{_j},y_{_{k}}) \in P(\mathring{\Omega}) \end{equation}
\begin{equation}\label{eq:sistema matricial indicial} (2\Delta y^2+2\Delta x^2)T_{_{jk}}-(\Delta y^2T_{_{j-1k}}+\Delta y^2T_{_{j+1k}}+ \Delta x^2T_{_{jk-1}}+\Delta x^2T_{_{jk+1}})=0, \hspace{0.5cm} \forall (x_{_j},y_{_{k}}) \in P(\mathring{\Omega}) \end{equation}
Al despejar $T_{_{jk}}$:
\begin{equation}\label{eq:moditer1} T_{_{jk}}=\dfrac{\Delta y^2(T_{_{j-1k}}+T_{_{j+1k}}) + \Delta x^2(T_{_{jk-1}} + T_{_{jk+1}})} {(2\Delta y^2+2\Delta x^2)}, \hspace{0.5cm} \forall (x_{_j},y_{_{k}}) \in P(\mathring{\Omega}) \end{equation}
\begin{equation}\label{eq:moditer1} T_{_{jk}}=\dfrac{\Delta y^2(T_{_{j-1k}}+T_{_{j+1k}}) + \Delta x^2(T_{_{jk-1}} + T_{_{jk+1}})} {(2\Delta y^2+2\Delta x^2)}, \hspace{0.5cm} \forall (x_{_j},y_{_{k}}) \in P(\mathring{\Omega}) \end{equation}
o equivalentemente:
\begin{equation}\label{eq:moditer2} T_{_{jk}}=\underbrace{\frac{1}{2}\left(\dfrac{\Delta y^{2}(a_{_j}+s_{j}) + \Delta x^{2}(a_{_k}+s_{k})}{\Delta y^{2}+\Delta x^{2}}\right)}_{\mathcal{D}}(T_{_{jk}}) = \mathcal{D}(T_{_{jk}}). \end{equation}
\begin{equation}\label{eq:moditer2} T_{_{jk}}=\underbrace{\frac{1}{2}\left(\dfrac{\Delta y^{2}(a_{_j}+s_{j}) + \Delta x^{2}(a_{_k}+s_{k})}{\Delta y^{2}+\Delta x^{2}}\right)}_{\mathcal{D}}(T_{_{jk}}) = \mathcal{D}(T_{_{jk}}). \end{equation}
El modelo \eqref{eq:moditer2} sugiere el siguiente algoritmo para calcular la
temperatura en la placa $\Omega$.
- Inicialización
- Iteración \[it_{_0}(T_{_{jk}}) = T_{_{jk}},\hspace{0.5cm} \forall (x_{_j}, y_{_k})\in P(\Omega)\] \[it_{_{r+1}}(T_{_{jk}}) = \mathcal{D}(it_{_r}(T_{_{jk}}))\]
- Terminación \[|it_{_{r+1}}(T_{_{jk}})-it_{_{r}}(T_{_{jk}})|<\epsilon.\]
\[T_{_{jk}}=T_{_0}(x_{_j}, y_{_k}),\hspace{0.5cm} \forall (x_{_j}, y_{_k}\in P(\mathring{\Omega}) \] \[T_{_{jk}}=h(x_{_j}, y_{_k}),\hspace{0.5cm} \forall (x_{_j}, y_{_k}\in P(\partial \Omega)\]
Ejemplo
El modelo iterativo anterior se puede aplicar para cualquier placa
rectangular. Por ejemplo, se puede considerar una placa de dimensiones $8\,cm$
en el $x$ y $10\,cm$ en el eje $y$. En este caso el dominio $\Omega$ es el
arreglo rectangular $[0,8]\times [0,10]$. Cuando se aplica una partición de
orden $m = 40$ y $n = 20$, los pasos para cada eje son $\Delta x=1/4$ y $\Delta
y=1/2$. Considere para todos los nodos de $\Omega$ una temperatura inicial de
$25\,$°C y las siguientes condiciones de frontera: Sobre la frontera superior
e izquierda una temperatura variable linealmente de $100\,$°C a $25\,$°C.
Al implementar la ecuación \eqref{eq:moditer2} en Python
con las condiciones anteriores se obtienen los siguientes resultados como se
muestra en la Figura 1.
Conclusiones
Con el modelo iterativo \eqref{eq:moditer2}, el cálculo de la
temperatura para cada punto $(x_{_j}, y_{_k})$ depende de los valores de los
puntos adyacentes por filas y por columnas, como puede apreciarse en la figura 2.
Las iteraciones solo se realizan sobre los puntos internos de la malla
mientras que la temperatura en los puntos de la frontera se mantiene constante.
Es posible actualizar la frontera con la temperatura de los puntos adyacentes
internos variando un poco el método dado que no siempre se cuentan con los puntos
adyacentes requeridos (ver figura 2). También se puede considerar
una ecuación en donde las condiciones de frontera varÃan en cada iteración debido
a interacciones del sistema con el entorno. Estos casos serán tratados en
próximos post.
Los errores en las aproximaciones de las derivadas en diferencias finitas
se deben al truncamiento en la serie de Taylor y depende de la naturaleza de la
ecuación en diferenciales parciales que se requiera resolver.
Referencias
- TJ Barth, M Griebel, DE Keyes, RM Nieminen, D Roose, and T Schlick. Finite difference computing with pdes.
- Randall J LeVeque. Finite difference methods for ordinary and partial differential equa- tions: steady-state and time-dependent problems. SIAM, 2007.
- M Necati Ozisik. Boundary value problems of heat conduction. Courier Corporation, 2013.
