Oscilaciones extremas y valores en la ecuación de calor en un anillo cerrado 1D

Estoy tratando de crear un modelo (muy) simple de difusión de calor a lo largo de una barra cerrada (por ejemplo, un anillo) con una condición inicial en la que un segmento está más caliente que el resto.

Estoy usando la siguiente fórmula:

T X t + 1 = T X t + Δ t α T X 1 t + T X 1 t 2 T X 1 t ( Δ X ) 2

con T siendo la temperatura, X posición, y α = λ ρ C siendo el coeficiente de transmisión dado por la conductividad térmica λ , densidad ρ y capacidad térmica C .

Me di cuenta de que cuando α acercarse a uno el sistema se vuelve inestable, con oscilaciones crecientes que producen temperaturas extremas por encima y por debajo del rango de temperatura inicial.

Aquí está el código R de una simulación para una barra larga de 100 segmentos con un segmento a 10 °C y el resto a 1 °C. El alfa está configurado para tener una conductividad térmica de una décima parte de ρ C . Esta configuración no produce inestabilidad y los resultados parecen realistas:

    library(ggplot2)
    
        v <- c(rep(1, 49), 10, rep(1, 50))
        
        out <- data.frame()
        
        for (i in 1:20) {
            d2T <- c(v[2:length(v)], v[1]) + c(v[length(v)], v[1:(length(v) - 1)]) - 2 * v
            
            # dx = dt = 1
            d2X <- 1
            dT <- 1
            
            out <- bind_rows(out, data.frame(
                val = v,
                x = 1:length(v),
                grad = d2T,
                iter = i
            ))

            v <- v + .1 * d2T/d2X
            
        }
        
        
        print(ggplot(out, aes(x, val)) +
            geom_line(alpha = .5, show.legend = F) +
            geom_point(aes(color = grad)) +
            scale_color_viridis_c(option = 'A') +
            facet_wrap(~ iter) +
            theme_minimal()) +
            labs(x = 'x', y = 'T', color = 'd2T/d2x') 

ingrese la descripción de la imagen aquí

Pero si α se establece en .5 o más, observamos tendencias inestables.

ingrese la descripción de la imagen aquí

El modelo se vuelve totalmente irreal como α se acerca a 1, con valores mucho más allá del rango de temperatura inicial de 1-10°C. Esto es α = .9 :

ingrese la descripción de la imagen aquí

Supongo que estropeé el modelo o es un efecto secundario de la pobre discretización de un proceso continuo. ¿Qué estoy haciendo mal?

No estás haciendo nada malo. El esquema explícito es inestable en valores superiores a 2. Para que su problema desaparezca, debe cambiar a un esquema implícito de tiempo como Euler inverso.

Respuestas (1)

El criterio de estabilidad de Von Neumann o Fourier falla para el esquema explícito de diferencias finitas en este contexto . el paso del tiempo Δ t debe ser más pequeño que ( Δ X ) 2 / 2 α para suprimir inestabilidades, donde Δ X corresponde a la discretización espacial y α es la difusividad térmica. Esta es la razón por la cual los pasos de unidad de tiempo y distancia tienden a fallar con α 0.5 .

Si está limitado en este aspecto, considere usar un método implícito de solución de diferencias finitas, que siempre es estable (pero más intensivo desde el punto de vista computacional).

¡Gracias! ¡No tenía ni idea! si por casualidad tienes alguna referencia sobre una solucion implementada te lo agradeceria, sino vere que encuentro!
incluso mejor si hay una solución exacta
Tal vez las Matemáticas de la difusión de Crank .
Mientras buscaba soluciones más apropiadas, parcheé el modelo para prohibir d T / d t mayor que ( T X 1 + T X + T X + 1 ) / 3 , es decir, la temperatura en un segmento no puede ir más allá de la temperatura media en su triplete de segmentos. Esto estabiliza la dinámica, eliminando tanto las oscilaciones como los valores extremos. Sé que es un truco, pero me preguntaba si se aproxima a algunos principios más generales (o en realidad es solo un truco feo).
lo siento, quise decir un d T / d t mayor que T X ( T X 1 + T X + T X + 1 ) / 3 en valor absoluto