Diferencia entre revisiones de «Ecuación del calor(Grupo Eau De Parfum(EDP)))»

Revisión del 20:37 6 mar 2024

1 Introducción

2 Conocimientos previos

Antes de comenzar será de gran ayuda introducir una serie de conceptos que utilizaremos a lo largo este trabajo:

• Ley del calor de Newton: describe la tasa de pérdida o ganancia de calor de un objeto en relación con el ambiente circundante.La ley se expresa comúnmente como:

[math] \frac{d T(t)}{d t} = - k (T(t) - T_{\mathrm{a}}) [/math],

siendo [math]T(t) [/math]la temperatura del objeto en un tiempo dado, [math]T_{\mathrm{a}}[/math] la temperatura ambiente y [math]k[/math] la constante de proporcionalidad.

• Principio de conservación de la energía calorífica: establece que la variación de energía calorífica sobre un cuerpo V se debe al balance entre el calor que entra y sale del cuerpo más una producción externa.

• Ley de Fourier:intaura que el calor fluye desde regiones de alta temperatura a regiones de baja temperatura, y la cantidad de flujo de calor depende de la diferencia de temperatura y las propiedades del material.La expresión matemática de la Ley de Fourier es:

[math] \mathbf{q} = - k {\nabla} T [/math],

donde [math] \mathbf{q}[/math] es el vector de flujo de calor por unidad de superficie y [math] {\nabla} T [/math] es el gradiente del campo de temperatura en el interior del material.

• Principio del máximo: Sea [math]u \in C^{2,1}(Q_T)\cap (\overline{Q_T})[/math] tal que [math]u_t - \Delta u \leq 0[/math] en [math]Q_T[/math]. Entonces u alcanza su máximo en la frontera parabólica [math]\partial _p Q_T[/math]:

[math]\max_{\overline{Q_T}} u = \max_{\partial _p Q_T} u[/math]

3 Primer problema

Consideramos una varilla metálica de longitud 1 m , que se encuentra aislada por su superficie lateral, de manera que la conducción de calor únicamente se produce en la dirección longitudinal. La temperatura inicial de la varilla es 0ºC. En el extremo izquierdo se consigue mantener la temperatura a 0 ºC mientras que en el derecho la temperatura es siempre de 1 ºC. Además tanto la conductividad térmica como el calor específico toman un valor constante de 1.Por lo que nuestro problema a resolver queda:

Este problema se resuelve mediante el método de variación de las constantes, pero un factor clave para su implementación es que las condiciones de contorno estén homogeneizadas, es decir, que sean igual a cero. Esto se conseguirá mediante un cambio de variable, en el que es imprescindible la llamada solución estacionaria. La solución estacionaria de nuestro problema es aquella que no varía con el tiempo [math] v=\lim_{t \to \infty} u(x,t)[/math]y por tanto, su derivada respecto al tiempo es nula. Nuestro sistema queda ahora:

Resolviendo obtenemos que [math] v(x)=x [/math] es nuestra solución estacionaria, cuya expresión gráfica es:

import matplotlib.pyplot as plt
import numpy as np
# Define el rango de x entre 0 y 1
x = np.linspace(0, 1, 100)  
# Define la función  estacionaria v(x) = x
v_x = x
# Crea la gráfica
plt.plot(x, v_x)
# Etiquetas de los ejes
plt.xlabel('x')
plt.ylabel('v(x)')
# Título de la gráfica
plt.title('Gráfica de la solución estacionaria v(x)=x')
# Muestra la gráfica
plt.show()

Una vez hallada la solución estacionaria, podemos homogeneizar nuestro problema mediante el cambio de variable [math]w(x,t)= u(x,t)-v(x)[/math].Por tanto, sustituyendo en el sistema inicial obtenemos:

Nuestro sistema es ya homogéneo,por lo que aplicamos el método de separación de variables (que consiste en expresar la solución como el producto de una función que depende únicamente de la variable espacial x y una función dependiente de la variable temporal) y obtenemos como solución homogénea:

Y deshaciendo el cambio, nuestra solución al problema original es por tanto:

Mostramos gráficamente la solución, observando que se cumplen las condiciones de frontera y además según avanza el tiempo, la solución deja de depender del tiempo y alcanza un estado estacionario, determinado justamente por la solución estacionaria obtenida anteriormente, [math] \ddot u(x)=x [/math]. Cabe destacar que en [math] t=0 [/math] obtenemos las oscilaciones correspondientes a la serie de Fourier de la condición inicial, que se aproxima mediante senos.

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# Parámetros
L = 1       # Longitud de la varilla
T = 1      # Tiempo final
D=1       # Coeficiente de conductividad térmica
N = 10      # Número de términos en la serie
x_points = 100  # Puntos en la dirección x
t_points = 100  # Puntos en la dirección t

# Espacios
x = np.linspace(0, L, x_points)
t = np.linspace(0, T, t_points)
X, T = np.meshgrid(x, t)

# Solución usando los primeros N términos de la serie
def u(x, t, N):
    sum = 0
    for k in range(1, N+1):
        lambda_k = k**2 * np.pi**2 
        #Calculasmos los coeficientes C_k
        C_k = 2 * (-1)**(k) / (k * np.pi)  
        sum += C_k * np.sin(k*np.pi * x) * np.exp(-lambda_k * t) 
    return x+sum

# Calcular la solución en la malla
Z = u(X, T, N)

# Graficar
fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(111, projection='3d')
surf=ax.plot_surface(X, T, Z, cmap='viridis')

ax.set_xlabel('Posición x')
ax.set_ylabel('Tiempo t')
ax.set_zlabel('Temperatura T(x,t)')
ax.set_title('Evolución de la Temperatura en la Varilla')

fig.colorbar(surf, shrink=0.5, aspect=5)
plt.show()

4 Flujo

Para estudiar el flujo de calor saliente y entrante en ambos extremos de la varilla a lo largo del tiempo, nos debemos basar en la ley de Fourier [math] \mathbf{q} = - k u_x(x,t) [/math] Por lo que necesitamos obtener la derivada de la solución con respecto al espacio:

Ahora, el flujo en los extremos es:

Gráficamente,

import numpy as np
import matplotlib.pyplot as plt

# Constantes
 # Conductividad térmica
L = 1  # Longitud de la varilla
N = 10  # Número de términos en el sumatorio
t = np.linspace(0, 1, 1000)  # Intervalo de tiempo

# Calcular los gradientes de temperatura en los extremos como sumatorios sobre k
def flujo_at_0(t, N):
    """Gradiente de temperatura en x=0 como un sumatorio que depende de k."""
    return -1-sum(2*(-1)**(k+1) * np.exp(-k**2 * np.pi**2 * t) for k in range(1, N+1))

def flujo_at_1(t, N):
    """Gradiente de temperatura en x=L como un sumatorio que depende de k."""
    return -1-sum((2) *np.exp(-k**2 * np.pi**2 * t)  for k in range(1, N+1))

# Vectores para almacenar los valores de los flujos de temperatura
flujo_at_0_vals = np.array([flujo_at_0(time, N) for time in t])
flujo_at_1_vals = np.array([flujo_at_1(time, N) for time in t])


# Graficar
plt.figure(figsize=(10, 6))
plt.plot(t, flujo_at_0_vals, label='Flujo en x=0')
plt.plot(t, flujo_at_1_vals, label='Flujo en x=1', linestyle='--')
plt.xlabel('Tiempo')
plt.ylabel('Flujo de Calor')
plt.title('Flujo de Calor en los Extremos de la Varilla ')
plt.legend()
plt.grid(True)
plt.show()

Al comenzar con un flujo de calor nulo en el extremo izquierdo, esto indica que hay una gran diferencia de temperatura entre ambos extremos de la varilla. El extremo izquierdo está perdiendo calor rápidamente hacia la zona más fría. A medida que el calor se va transfiriendo por la varilla, el gradiente de la temperatura cambia y el extremo derecho comienza a calentarse, debido al flujo de calor que proviene del otro extremo. Este cambio conduce a que el flujo en ambos extremos comience a estar en equilibrio. Esta estabilización de flujo en -1 refleja que se ha alcanzado un estado estacionario en el sistema, indicando que cualquier calor entrante en la varilla, se compensa con el calor saliente, manteniendo una distribución de temperatura uniforme.

Este comportamiento refleja la Ley de Fourier, donde al principio la gran diferencia de temperatura induce un flujo de calor bastante significativo, que se va estabilizando llegando a un flujo de -1 en ambos extremos, alcanzando así un gradiente de temperatura constante, con la misma cantidad de calor saliendo hacia fuera desde ambos extremos.

5 Cambio de coeficiente de conductividad térmica

Hemos estudiado el problema de Cauchy para un coeficiente de difusión de 1, pero vamos a observar cómo sería para un valor cualquiera. Estableciendo que D es el coeficiente de difusión, nuestro problema es:

Debemos reescalar en tiempo, para facilitarnos los cálculos, por lo que realizamos el cambio de variable [math]\tau = \alpha t, \tau \in\mathbb{R} [/math]. Por tanto, [math] u(x,t)=v(x,\tau)[/math] y entonces [math] \alpha v_{\tau}-Dv_{xx}=0[/math] .Tomando [math] \alpha = D[/math] , obtenemos el problema

Cuya solución es

Deshaciendo el cambio de variable, nuestra solución al problema original es:

6 Cambio de condiciones

Supongamos ahora que la temperatura en el extremo derecho es [math]0 Cº[/math] y, además, inicialmente la temperatura de la barra viene dada por la función [math]u(x,0) = max\{0, 1−4|x−1/2|\}[/math]. Entonces, obtenemos este sistema:

Supongamos que, en tiempo infinito, \( t \rightarrow \infty \), se llega a una solución estacionaria y que, por tanto, ya no varía en tiempo. De modo, que tomando esta suposición llegamos al siguiente sistema con \( u(x,t) \xrightarrow{\text{t tiende a }\infty} v(x) \).

Cuya solución es la constante 0. De modo que el sistema homogeneizado, coincide con el original, pues \( w(x,t)=u(x,t)-v(x) =u(x,t)-0=u(x,t)\) Antes de exigir la condición inicial, la solución coincide con la del primer apartado, pues nos encontramos ante la misma ecuación diferencial, salvo por la función \( u(x,0)\), la cual no toma parte hasta el final de la resolución de esta suponiendo separación de variables. De modo que

[math]u(x,t) = \sum_{k=1}^{\infty} C_k e^{(-k^2\pi ^2 t )}\sin(k \pi x)[/math]

Para hallar el valor de los coeficientes [math] C_k[/math], buscamos la serie de Fourier de la función [math] max\{0, 1−4|x−1/2|\}[/math] extendida de forma impar. Esto es posible, ya que

[math] \int^1_0 (max\{0, 1−4|x−1/2|\})^2 dx=0.1667 \lt \infty[/math]

El cálculo de los coeficientes se ha realizado numéricamente aproximando por la fórmula del trapecio. Los coeficientes se pueden consultar ejecutando el siguiente código en Python

import numpy as np 
def f(x):   # Defino función que se quiere aproximar
    return np.maximum(0, 1 - 4 * np.abs(x - 1/2))

def fourier_coefficients(f, n, x):
    coefficients = np.zeros(n) # Guardo los coeficientes en una lista
    for k in range(1,n+1):     # Compiezo en k=1, pues sin(0)=0
        coefficients[k-1] = 2 * np.trapz(f(x) * np.sin(k * np.pi * x), x)
    return coefficients

fourier_coefficients(f, 10, np.linspace(0, 1, 1000))

De modo que ya se ha calculado la solución de este nuevo problema. Visualmente, responde a la siguiente gráfica, resultado de ejecutar este código:

import sympy as sp
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# Parámetros
L = 1       # Longitud de la varilla
T = 1       # Tiempo final
N = 10      # Número de términos en la serie
# Espacios
x = np.linspace(0, L, 1000) # Puntos en la dirección x
t = np.linspace(0, T, 1000) # Puntos en la dirección t
X, T = np.meshgrid(x, t)    # Generar una cuadrícula de puntos

def f(x):    # Defino mi función condición inicial
    return np.maximum(0, 1 - 4 * np.abs(x - 1/2))

def fourier_coefficients(f, n, x): # Coeficientes, recuperado del anterior trabajo
    coefficients = np.zeros(n)
    for k in range(1,n+1):         # Comenzamos en 1, pues sen(0)=0
        coefficients[k-1] = 2 * np.trapz(f(x) * np.sin(k * np.pi * x), x)
    return coefficients.tolist()   # Convertir a lista
coef=fourier_coefficients(f, N, x)

def u(x, t, n): # Definir la solución en términos de la serie de Fourier con N términos
    sum = 0                 # Almacenamos la suma (resultado) hasta el término 10
    for k in range(1, n+1): # Para cada k calculo el lambda y sumo su coeficiente respectivo
        lambda_k = k* np.pi
        sum += coef[k-1] * np.sin(lambda_k * x) * np.exp(-lambda_k**2 * t)
    return sum

Z = u(X, T, N) # Calcular los valores de la función para cada par (x, t)

# Graficar la función en 3D
fig = plt.figure(figsize=(12, 8))  # Cambio aquí para hacer la figura más grande
ax = fig.add_subplot(111, projection='3d')
surf=ax.plot_surface(X, T, Z,cmap='viridis')

ax.set_xlabel('Posición x [m]')
ax.set_ylabel('Tiempo t [s]')
ax.set_zlabel('Temperatura T(x,t) [ºC]')
ax.set_title('Evolución de la Temperatura en la Varilla')
fig.colorbar(surf, shrink=0.5, aspect=5)  #Añado una guía de colores
surf.set_clim(0, 1)                       # Le introduzco cota inf/sup a los colores
plt.show()

​ Vemos como la temperatura sufre una gran variación en tiempos cercanos a 0, y rápidamente tiende a su estado estacionario, [math]v(x)=0[/math]. Para poder visualizarlo mejor mostraremos una gráfica restringida a [math]t \in [0,0.05][/math] Al principio, u es nulo fuera del origen. Una vez que [math]t\gt0[/math] u se hace positivo en todas partes, lo que significa que la temperatura se difunde casi instantáneamente a lo largo del eje x, resultando en una velocidad de propagación infinita. \ref Terminamos el estudio de esta ecuación de difusión con el estudio del flujo, la inversa de la variación de la variable x de la solución. Para ello derivamos la solución que hemos obtenido. En términos de [math]C_k[/math], puesto que no son dependientes de x, obtenemos: [math]\frac{\partial}{\partial x}u(x,t)=\sum_{k=1}^{\infty} C_k e^{(-k^2\pi ^2 t )}\cos(k \pi x)\frac{1}{k \pi}[/math] como flujo de la ecuación homogénea. Adjuntamos el código con el cual se ha calculado la expresión de forma numérica y su respectiva representación gráfica:

import sympy as sp
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# Parámetros
L = 1       # Longitud de la varilla
T = 1       # Tiempo final
N = 10      # Número de términos en la serie
# Espacios
x = np.linspace(0, L, 1000) # Puntos en la dirección x
t = np.linspace(0, T, 1000) # Puntos en la dirección t
X, T = np.meshgrid(x, t)

def f(x):   # Defino la funcion de la ci para hallar los coef
    return np.maximum(0, 1 - 4 * np.abs(x - 1/2))

def fourier_coefficients(f, n, x): # Hallo los coef de Fourier  
    coefficients = np.zeros(n)
    for k in range(1,n+1): 
        coefficients[k-1] = 2 * np.trapz(f(x) * np.sin(k * np.pi * x), x)
    return coefficients.tolist()   # Convertir a lista
coef=fourier_coefficients(f, N, x)

def u(x, t, n):  # Definir la derivada en términos de la serie de Fourier con N términos
    sum = 0
    for k in range(1, n+1):
        lambda_k = k* np.pi
        sum += coef[k-1] *np.cos(k * np.pi * x)/(k * np.pi) * np.exp(-lambda_k**2 * t)
    return sum

# Graficamos
plt.figure(figsize=(8, 6))
plt.plot(t, u(0,t,N), label='Flujo en x=0', color='blue')  # Cálculo del flujo en x=0
plt.plot(t, u(1,t,N), label='Flujo en x=1', color='red')   # Cálculo del flujo en x=1
plt.xlabel('Tiempo [s]')
plt.ylabel(r'Flujo de Calor [$\frac{^\circ C}{m^2 \cdot s}$]')
plt.title('Flujo de Calor en los Extremos de la Varilla ')
plt.legend()
plt.grid(True)
plt.show()