Ver el código de «Dibujar un sólido (Grupo 5)»

{{ TrabajoED | Ecuación del calor. Grupo GRwM | [[:Categoría:EDP|EDP]]|[[:Categoría:EDP23/24|2023-24]] | Guillermo Gómez Tejedor,  Marina Jiménez Barrantes y Rocío Tajuelo Díaz}}

= Introducción=
En el trabajo que se presenta a continuación, vamos a estudiar la ecuación del calor en una dimensión. Para ello, vamos a considerar una varilla metálica que se encuentra aislada por su superficie lateral, de modo que la conducción de calor se produzca en la dirección longitudinal. 
Partiendo de esto y añadiendo ciertas condiciones de frontera, que iremos modificando a lo largo del trabajo, vamos a calcular la solución de la ecuación del calor y la vamos a representar en los distintos casos.

Posteriormente y empleando la solución fundamental de la ecuación del calor, estudiaremos la solución para dimensión 2.

Cabe destacar que los cálculos y visualizaciones han sido programados con MatLab.

= Conceptos previos=
En esta sección, vamos a presentar algunos conceptos esenciales para comprender la obtención del sistema EDP que permite modelizar el problema.

'''Flujo de calor''': Se define flujo de transferencia de calor, o simplemente flujo de calor, a la cantidad de calor 
PONER BIE

'''Ley de Fourier''': Establece que el flujo de transferencia de calor por conducción es proporcional y de sentido contrario al gradiente de temperatura (<math> T </math>) en esa dirección. En nuestro caso, al trabajar en una dimensión, <math> \nabla T(\textbf{x},t)= T_{x} </math> y la ley de Fourier es
<center><math> q=-k\cdot T_{x} 

</math></center>
siendo <math> q </math> el flujo de transferencia de calor y <math> k </math> el coeficiente de conductividad térmica, que es un valor que indica la capacidad de un material para transferir calor a través de él cuando existe una diferencia de temperatura.

'''Energía calorífica''': La energía calorífica <math> e </math> es una forma de energía que se produce como resultado del movimiento de átomos en un objeto y es una medida de la cantidad total de energía cinética de estos. La expresión que relaciona la temperatura del objeto con la energía calorífica es:

<center><math> e= c \cdot T

</math></center>
siendo <math> c </math> el calor específico, que se define como la cantidad de calor que hay que suministrar a la unidad de masa de una sustancia para elevar su temperatura en una unidad.

'''Principio de conservación de la energía''':Establece que la cantidad total de energía térmica para cierto volumen de control permanece constante cuando se tiene en cuenta el flujo de calor entrante y saliente con respecto al tiempo.

Por otro lado, vamos a comprobar que el principio del máximo se verifica en este problema:

'''Principio del máximo''': sea <math> u \in C^{2,1} (Q_T) \cap C(\overline{Q_T})</math> tal que <math>u_t - \Delta u \leq 0 </math> en <math> Q_T = \Omega \times (0,T)</math>. Entonces <math>u</math> alcanza un máximo en la frontera parabólica, es decir, <math> \max \limits_{(x,t) \in \overline{Q_T}} u = \max \limits_{(x,t) \in \partial _P Q_T} u </math>.


REVISAAAR

CONECTARLO CON EL TRABAJO ANTERIOR DE SERIES DE FOURIER (INTENTAR PONER LINK O ALGO)

=Planteamiento del problema=
Para comenzar con el estudio de la ecuación del calor, primero debemos plantear el problema a resolver, que involucra esta ecuación junto a ciertas condiciones de frontera y condición inicial.
Como ya se ha mencionado, en nuestro estudio vamos a considerar una varilla metálica que se encuentra aislada por su superficie lateral. De esta manera, la conducción de calor se produce únicamente en la dirección longitudinal.

Además, vamos a considerar que la temperatura inicial de la varilla es 0 ºC. También vamos a fijar la temperatura  en el extremo izquierdo en 0ºC y en el derecho a 1ºC.

Teniendo en cuenta el principio de conservación la energía y la definición de la energía en función de la temperatura, así como las condición inicial y de frontera, se obtiene el siguiente sistema de EDP:
<center><math>\left \{ \begin{array}{ll} 

\frac{\partial T}{\partial t}-\frac{k}{c} \frac{\partial ^2 T}{\partial x^2}=0 & \quad 0 < x < 1, 0 < t < T, \\ 

T(x,0)=0, & \quad 0 < x < 1, \\ 

T(0,t)=0, & \quad 0 < t < T, \\

T(1,t)=1 & \quad 0 < t < T. 

\end{array} \right. 

</math></center>


Suponemos que tanto la conductividad térmica <math>k</math> como el calor específico <math>c</math> toman el valor constante <math>1</math>. Entonces, el sistema de EDP final queda de la siguiente forma:

<center><math>\left \{ \begin{array}{ll} 

\frac{\partial T}{\partial t}-\frac{\partial ^2 T}{\partial x^2}=0 & \quad 0 < x < 1, 0 < t < T, \\ 

T(x,0)=0, & \quad 0 < x < 1, \\ 

T(0,t)=0, & \quad 0 < t < T, \\

T(1,t)=1 & \quad 0 < t < T. 

\end{array} \right. 

</math></center>

=Resolución del sistema EDP=
Una vez planteado el sistema, procedemos a resolverlo. Para ello, con el objetivo de homogeneizar las condiciones de frontera, vamos a comenzar obteniendo la solución estacionaria.



Para calcular la ecuación del estado estacionario, vamos a tomar el tiempo <math> t </math> infinito. Esto implica que la variación de la temperatura con respecto al tiempo desaparezca, de modo que la ecuación del sistema HACER REFERENCIA es ahora <math> T_{xx}=0</math>.

Considerando además las condiciones <math> T(0)=0</math> y <math> T(1)=0</math>, provenientes de las condiciones frontera del problema original, se tiene que la solución de la ecuación del estado estacionario es:
<center><math> \begin{array}{ll}  T_{est}(x)=x, & \quad 0 < x < 1 \end{array}. </math></center>

En la siguiente gráfica, se representa esta solución:










Consideramos ahora el problema equivalente con condiciones de frontera homogéneas, donde se define <math> T_{hom}(x,t):= T_{est}(x)-T(x,t)</math>:


<center><math>\left \{ \begin{array}{ll} 

\frac{\partial T_{hom}}{\partial t}-\frac{\partial ^2 T_{hom}}{\partial x^2}=0 & \quad 0 < x < 1, 0 < t < T, \\ 

T_{hom}(0,t)=0, & \quad 0 < t < T, \\

T_{hom}(1,t)=0 & \quad 0 < t < T, \\ 

T_{hom}(x,0)=x, & \quad 0 < x < 1.

\end{array} \right. 

</math></center>


Obtenemos la solución de este problema aplicando separación de variables. Para ello, vamos a suponer que la solución es de la forma:

<center><math> T_{hom}(x,t)=u(x)v(t) </math></center>

Haciendo las cuentas pertinentes se tiene que

<center><math> T_{hom}(x,t)=\sum_{k=1}^{\infty} \frac{2(-1)^{k+1}}{k \pi}\sin{(k \pi x)} e^{-k^2\pi^2t} </math>.</center>

A continuación se representa en el intervalo de tiempo <math> t \in [0,1] </math>, tomando los 10 primeros términos de la serie:









Finalmente, la solución del problema original es <center><math> T(x,t)=x-\sum_{k=1}^{\infty} \frac{2(-1)^{k+1}}{k \pi}\sin{(k \pi x)} e^{-k^2\pi^2t} </math></center>

=Interpretación del flujo en los extremos=
Para finalizar el estudio de la solución del sistema de EDP PONER REFERECIA, vamos obtener e interpretar el flujo de calor saliente y entrante en ambos extremos a lo
largo del tiempo. 

Para ello, por la propia definición de flujo, debemos obtener la derivada de la solución con respecto a la variable x, que es:

<center><math> T_{x}(x,t)=\sum_{k=1}^{\infty} 2(-1)^{k+1}\cos{(k \pi x)} e^{-k^2\pi^2t} </math>.</center>
Además, para su interpretación, debemos tener en cuenta que el flujo es positivo cuando su sentido es el del eje longitudinal, es decir, de izquierda a derecha. De modo que será entrante o saliente dependiendo de su signo en los extremos. Calculamos así la derivada en los extremos:

<center><math> T_{x}(0,t)=\sum_{k=1}^{\infty} 2(-1)^{k+1} e^{-k^2\pi^2t} </math>.</center>
<center><math> T_{x}(1,t)=\sum_{k=1}^{\infty} -2 e^{-k^2\pi^2t} </math>.</center>

Si consideramos únicamente el primer término de la serie, obtenemos que en ambos extremos el flujo es saliente. Pues, por un lado se tiene que 
<center><math> T_{x}(0,t) = PONER >0 </math>.</center>

Lo que implica que el flujo sea negativo en el extremo izquierdo, y que por tanto vaya en contra del sentido del eje longitudinal, es saliente. 
Por otro lado, se tiene que 

<center><math> T_{x}(1,t) < 0 </math>.</center>
es decir, el flujo en el lado derecho va en la dirección del eje longitudinal, luego también es saliente.


REVISAR : Sin embargo, a pesar de que los primeros términos de la serie tienen una mayor influencia cuando <math> t </math> es cercano a cero, a medida que el tiempo aumenta, el resto de términos adquieren relevancia y por tanto no son despreciables. 
Por esta razón, hemos aumentado el número de términos considerados. De esta manera, obtenemos valores de flujo más precisos.

En la siguiente gráfica se ilustra el flujo para 10 y 1000 términos de la serie:







Como podemos observar, y pese a lo que habíamos deducido al considerar un único término de la serie, obtenemos que el flujo es entrante en el lado izquierdo y saliente en el derecho.

=Solución al considerar otro coeficiente de conductividad=
Si consideramos ahora el coeficiente de conductividad térmica igual a <math>\frac{1}{2}</math>, se tiene el siguiente problema:

<center><math> \left \{ \begin{array}{ll}
\frac{\partial T}{\partial t}-\frac{1}{2} \frac{\partial ^2 T}{\partial x^2}=0 & \quad 0 < x < 1, 0 < t < T, \\ 

T(x,0)=0, & \quad 0 < x < 1, \\ 

T(0,t)=0, & \quad 0 < t < T, \\

T(1,t)=1 & \quad 0 < t < T. 

\end{array} \right. 

</math></center>

Siguiendo los mismos pasos que el problema inicial, obtenemos que la solución es
<center><math> T_{x}(x,t)=\sum_{k=1}^{\infty} 2(-1)^{k+1}\cos{(k \pi x)} e^\frac{-k^2\pi^2t}{2} </math>.</center>


Además, con el objetivo de comparar esta solución con la obtenida en el anterior apartado, vamos a dibujar la diferencia entre la solución y el estado estacionario en <math>x= \frac{1}{2} </math> para <math> t \in [0, 1]</math>, es decir, <math> u(1/2, t) - u_{est}(1/2)</math>.




PONEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEER




Uno de los aspectos que se puede observar es que, al considerar un coeficiente de conductividad menor, en este caso <math>\frac{1}{2}</math>, aumenta el tiempo necesario para que se alcance el valor <math> -1</math> en el flujo. Esto se debe a que, como se mencionaba en la sección de conceptos previos, la constante <math> k</math> indica la capacidad del material para transferir calor a través de él cuando existe una diferencia de temperatura. Por tanto, si este valor es menor, el tiempo en alcanzar el valor del flujo para la solución estacionaria es mayor, concretamente el doble.

=Solución al considerar otra condición inicial=

==Replanteamos el problema==
En este nuevo caso, vamos a suponer que la temperatura en el extremo derecho es también <math>0</math>ºC, es decir, tenemos unas condiciones frontera homogéneas. Además, vamos a cambiar la condición inicial por una nueva función. Mantendremos tanto la conductividad térmica <math>k</math> como el calor específico <math>c</math> con el valor constante <math>1</math>, siendo así el sistema EDP:

<center><math>\left \{ \begin{array}{ll} 

\frac{\partial T}{\partial t}- \frac{\partial ^2 T}{\partial x^2}=0 & \quad 0 < x < 1, 0 < t < T, \\ 

T(x,0)=max\{0,1-4|x-\frac{1}{2}|\}, & \quad 0 < x < 1, \\ 

T(0,t)=T(1,t)=0 & \quad 0 < t < T

\end{array} \right. 

</math></center>

==Resolución del sistema EDP==
Al haber considerado condiciones de frontera nulas se tiene que la solución estacionaria del problema es <math>0</math>. Por tanto, no es necesario el proceso de homogeneización.

Aplicando el método de separación de variables obtenemos que la solución de la ecuación del calor es de la forma:


<center><math> u(x,t)=\sum_{k=1}^{\infty} c_{k}\sin{(k \pi x)} e^{-k^2 \pi^2t} </math>.</center>


En este caso, en lugar de obtener los coeficientes de Fourier de forma analítica, vamos a aproximar las integrales por la fórmula del trapecio.
Teniendo esto en cuenca, la gráfica de la solución para <math>t \in [0,1]</math> es:


Como podemos observar, la solución a medida que avanza el tiempo se aproxima al estado estacionario.
También es destacable que la solución es estrictamente positiva para <math> t > 0 </math>. Esto se interpreta como que la difusión ’transporta’ la energía con velocidad infinita.

==Interpretación del flujo en los extremos==
Para finalizar el estudio de la solución del sistema de EDP PONER REFERECIA, vamos obtener e interpretar el flujo de calor saliente y entrante en ambos extremos a lo
largo del tiempo. 

Para ello obtenemos la derivada de la solución con respecto a la variable x, que es:

<center><math> T_{x}(x,t)=\sum_{k=1}^{\infty} c_{k} \cos{(k \pi x)} e^{-k^2\pi^2t} </math>.</center>

Como ya pudimos apreciar en el anterior cálculo del flujo, no podemos restringir nuestro estudio a considerar el primer término de la serie. Por ello, también en esta ocasión vamos a considerar 10 y 1000 términos de la serie. Siendo el resultado el siguiente:

=Cambio de condiciones frontera=
Ahora suponemos que el extremo derecho de la barra está aislada térmicamente, en lugar de que la temperatura sea de <math>0 ~ ºC</math>. Entonces, tenemos el siguiente problema:

<center><math>\left \{ \begin{array}{ll} 

\frac{\partial T}{\partial t}- \frac{\partial ^2 T}{\partial x^2}=0 & \quad 0 < x < 1, 0 < t < T, \\ 

T(0,t)=0, & \quad 0 < t < T, \\

T_{x}(1,t)=0 & \quad 0 < t < T, \\ 

T(x,0)=max \{0,1-4|x-1/2|\}, & \quad 0 < x < 1.

\end{array} \right. 

</math></center>

que tiene condiciones de frontera nulas y por tanto, no es necesario el proceso de homogeneización.

A continuación, obtenemos la solución de este problema aplicando separación de variables. En este caso, obtenemos que la solución es de la forma:


<center><math> u(x,t)=\sum_{k=1}^{\infty} c_{k}\sin{((\frac{\pi}{2}+k\pi)x)} e^{-(\frac{\pi}{2}+k\pi)^2t} </math>.</center>

Con esta solución, no podemos emplear la base trigonométrica usual <math>\{ \frac{1}{2}, sen(k \pi x), cos(k \pi x) \}_{k=1}^{\infty} </math>. En su lugar, vamos a trabajar con la base <math>\{ \frac{1}{2}, sen((\frac{\pi}{2}+k\pi)x), cos((\frac{\pi}{2}+k\pi)x) \}_{k=1}^{\infty} </math>. Más concretamente, al extender de forma impar la función <math>max \{0,1-4|x-1/2|\}</math>, vamos a tener la base <math>\{sen((\frac{\pi}{2}+k\pi)x)\}_{k=1}^{\infty} </math> pues  forma un conjunto completo en <math>[0,1]</math> y es ortogonal. Comprobemos esto último:


<center><math> < sen((\frac{\pi}{2}+n\pi)x), sen((\frac{\pi}{2}+m\pi)x)>_{L^2(0,1)}= \int_{0}^{1}sen((\frac{\pi}{2}+n\pi)x) \cdot sen((\frac{\pi}{2}+m\pi)x) dx= \left[ \frac{sen(((n-m) \pi )x)}{\pi (2n-2m)}- \frac{sen((\pi(n+m+1)x)}{\pi(2n+2m+2)} \right]_{0}^{1}=0 </math>.</center>


Finalmente, para obtener la solución, obtenemos los coeficientes de Fourier. De modo que la solución es :


PONEEER


Sin embargo, para la representación en Matlab, vamos a emplear el método del trapecio.

=Principio del máximo=
Consideramos el problema del apartado anterior. Se verifica que <math> T \in C^{2,1} (Q_T) \cap C(\overline{Q_T})</math> y <math>T_t - \Delta T \leq 0 </math>, donde <math> Q_T=(0,1)\times(0,T)</math>. Entonces, por el principio del máximo, se cumple que  <math>T</math> alcanza un máximo en la frontera parabólica, es decir, <math> \max \limits_{(x,t) \in \overline{Q_T}} T = \max \limits_{(x,t) \in \partial _P Q_T} T </math>. Este resultado se puede apreciar en la siguiente gráfica:


PONER GRÁFICA DEL MÁXIMO


Como podemos observar, el máximo, cuyo valor es PONER, se alcanza en PONER, i.e, en la frontera.

=Soluciones de la ecuación del calor usando la solución fundamental=
En esta sección vamos a estudiar la solución de la ecuación del calor en dimensión <math> n</math> empleando la solución fundamental, que viene dada por la siguiente expresión:

<center><math> \Phi_D (x,t) = \frac{1}{(4 \pi Dt)^{n/2}}e^{-\frac{|x|^2}{4Dt}}</math></center>

donde <math>D=\frac{k}{c}</math>.

==Solución fundamental de la ecuación del calor en dimensión <math>1</math>==
Como ejemplo en dimensión <math>n=1</math>, vamos a considerar <math>x \in [-1,1]</math>, <math>t \in [10^{-2},1]</math> y <math>k=c=1</math>, es decir, <math>D=1</math>. Por tanto, la expresión de la solución fundamental queda de la siguiente manera:

<center><math> \Phi_1 (x,t) = \frac{1}{(4 \pi t)^{1/2}}e^{-\frac{|x|^2}{4t}}</math>.</center>

A continuación podemos ver la representación gráfica de dicha función:

PONER GRÁFICA

==Ecuación del calor en una dimensión en el semiespacio <math> x>0</math>.==

En esta sección, vamos a considerar el problema:

<center><math>\left \{ \begin{array}{ll} 

\frac{\partial T}{\partial t}- \frac{\partial ^2 T}{\partial x^2}=0 & \quad 0 < x , 0 < t  \\ 

T(0,t)=1 & \quad 0 < t  \\ 


T(x,t)=0 & \quad 0 < t , x \rightarrow \infty \\ 

T(x,0)=0 & \quad x>0  


\end{array} \right. 

</math></center>

La forma de resolver esta ecuación diferencial es muy diferente, y reside en suponer que <math> T(x,t)= u(\frac{x}{\sqrt{t}}) </math>. 
Recalculamos la ecuación diferencial para <math>u</math>:

<center><math>

0=\frac{\partial T}{\partial t}- \frac{\partial ^2 T}{\partial x^2}=u_t(\frac{x}{\sqrt{t}})-u_{xx}(\frac{x}{\sqrt{t}})= -\frac{x}{2t\sqrt{t}}u'(\frac{x}{\sqrt{t}}) -\frac{1}{t}u''(\frac{x}{\sqrt{t}}) \Longleftrightarrow 0= \frac{1}{2}\frac{x}{\sqrt{t}}u'(\frac{x}{\sqrt{t}}) + u''(\frac{x}{\sqrt{t}})

</math></center>
Llamamos <math> \zeta = \frac{x}{\sqrt{t}} </math> y obtenemos que <math>u</math> cumple la siguiente ecuación diferencial:
<center><math>
0=\frac{1}{2}\zeta u'(\zeta) +u''(\zeta) 
</math></center>
Resolviendo la ecuación diferencial, obtenemos que <math> u(\zeta)=\int^{\zeta}_0e^{-\frac{s^2}{4}}ds = \sqrt{\pi}\left(\frac{2}{\sqrt{\pi}}\int^{\frac{\zeta}{2}}_0e^{-z^2}dz \right) = \sqrt{\pi}~erf(\frac{\zeta}{2})</math> luego <math> T_1(x,t)=u(\frac{x}{\sqrt{t}})= \sqrt{\pi} ~ erf(\frac{x}{2\sqrt{t}}) </math> (siendo <math>erf(x)=\frac{2}{\sqrt{\pi}}\int^x_0e^{-z^2}dz</math>). Sin embargo, esta solución cumple:

<center><math>\left \{ \begin{array}{ll} 

\frac{\partial T_1}{\partial t}- \frac{\partial ^2 T_1}{\partial x^2}=0 & \quad 0 < x , 0 < t  \\ 

T(0,t)=0 & \quad 0 < t  \\ 

T(x,t)= \sqrt{\pi} & \quad 0 < t , x \rightarrow \infty \\ 

T(x,t)=\sqrt{\pi} & \quad x>0 , x \rightarrow 0 

\end{array} \right.

</math></center>

Solo queda hacer un simple cambio para que cumpla las condiciones propuestas al inicio: <math> T(x,t) = 1 - \frac{1}{\sqrt{ºpi}}T_1(x,t)</math> . Luego la solución que buscamos es:

<center><math> T(x,t) = 1 - \frac{2}{\sqrt{\pi}}\int^{\frac{x}{2\sqrt{t}}}_0e^{-z^2}dz  = 1 - erf(\frac{x}{2\sqrt{t}})  </math></center>

== Solución del calor con dato inicial <math> u_{0}(x) = 1_{[-1,1]} </math>==

==Solución fundamental de la ecuación del calor en dimensión 2==
En último lugar, vamos a estudiar la solución fundamental de la ecuación del calor en dimensión <math>2</math>. En este caso, tomamos el intervalo espacial <math>(x_1,x_2)=[-1,1]^2</math> y <math>D=1</math>. La expresión de la solución fundamental es la siguiente:

<center><math> \Phi_1 (\vec{x},t) = \frac{1}{4 \pi t}e^{-\frac{|\vec{x}|^2}{4t}}</math>.</center>

A continuación se representa para diferentes instantes de tiempo (<math> t=0.001,~t=0.01</math> y <math> t=0.1</math>):

PONER GRÁFICAS

=Referencias=

* [https://www.thermal-engineering.org/es/que-es-la-ley-de-conduccion-termica-de-fourier-definicion/ ¿Qué es la Ley de conducción térmica de Fourier? Definición. Nick Connor]