Modelo não isotérmico da variação da densidade atmosférica conforme o aumento da altitude

A atmosfera terrestre é dividida em camadas que variam em temperatura, densidade, pressão e altitude. Hoje vamos focar em uma dessas variáveis: a densidade do ar; ela é definida como a razão entre massa e volume $(\frac{Kg}{m^2})$:

$$\rho = \frac{m}{V}$$

Quanto maior a altitude e mais distante da superfície, menor será a densidade. Mas por que?

Podemos atribuir à essa queda três fatores:

1. À medida que aumentamos a altura, há menos “ar” acima de nós, pressionando-nos, diminuindo a pressão e, assim, a densidade;
2. Pela gravidade diminuir conforme aumenta-se a altitude, o peso efetivo do ar decai, diminuindo também a densidade;
3. Em maiores altitudes, a temperatura do ar fica mais fria, sendo assim, menos denso.

Com isso, nota-se que não é simples compreender a divisão atmosférica, visto que há tantos fatores envolvidos. Portanto, a partir de um modelo no qual a atmosfera é uma única camada isotérmica, ou seja, não possui variação de temperatura, e em que há uma diminuição exponencial da densidade com a altitude, vamos complexificar nosso estudo.

A função simplificada original seria:

$$\rho = \rho_0 e^{\frac{-h}{H}}$$

Para $\rho_0 = 1225 \frac{Kg}{m^3}$ como a densidade ao nível do mar. Vamos criar um código que modifique esse pensamento inicial e mostre-nos a diminuição da densidade considerando as mudanças de temperatura!

Primeiro estabelecemos duas constantes, a gravidade e a constante dos gases. Elas serão necessárias para calcular a densidade com a equação de estado dos gases ideais:

$$\rho = \frac{P}{RT}$$

g0 = 9.80665      # gravidade (m/s²)
R = 287.05        # constante dos gases (J/kg·K)

Então estabelecemos as camadas atmosféricas até 80 km, conforme o padrão internacional. Cada camada terá altitude base (m), temperatura na base, um gradiente de temperatura (K/m ou K/km) e a altitude final. Também vamos assumir um intervalo para a altitude, com 800 valores entre 0 e 80 km.

# Camadas: (base_alt, base_temp, lapse_rate, top_alt)
layers = [
   (0,     288.15, -0.0065, 11000),   # Troposfera
   (11000, 216.65,  0.0,    20000),   # Tropopausa
   (20000, 216.65,  0.001,  32000),   # Estratosfera 1
   (32000, 228.65,  0.0028, 47000),   # Estratosfera 2
   (47000, 270.65,  0.0,    51000),   # Estratopausa
   (51000, 270.65, -0.0028, 71000),   # Mesosfera 1
   (71000, 214.65, -0.002,  80000),   # Mesosfera 2
]


# Altitudes e inicialização
altitudes = np.linspace(0, 80000, 800)
densities = np.zeros_like(altitudes)

Definimos as variáveis iniciais da primeira camada, que serão atualizados ao que aumentamos a altitude, e criamos um loop para cada camada da atmosfera. Esse loop será responsável por calcular a densidade naquela faixa, dependendo se há um gradiente de temperatura por km ou não.

# Condições iniciais
p_base = 101325  # Pa
T_base = 288.15  # K
h_base = 0       # m


alt_pointer = 0
# Loop sobre as camadas
for layer in layers: #Iteramos
   h_start, T_start, L, h_end = layer
   while alt_pointer < len(altitudes) and altitudes[alt_pointer] <= h_end: # Dentro de cada camada
       h = altitudes[alt_pointer]
       delta_h = h - h_base
# Cálculo da densidade a depender se naquela camada há um gradiente ou não
       if L == 0: 
           # Isotérmica
           T = T_base
           p = p_base * np.exp(-g0 * delta_h / (R * T))
       else:
           # Com gradiente
           T = T_base + L * delta_h
           p = p_base * (T / T_base) ** (-g0 / (R * L))


       rho = p / (R * T) # Para calcular a densidade

Os dados T_base e p_base são atualizados no fim de cada camada, a fim de que sirvam como base para a próxima faixa. Então, com o modelo estabelecido, agora vamos plotar essas relações em um gráfico, com faixas coloridas indicando as camadas da atmosfera estudadas:

Dados simulados, não são baseados em medidas experimentais. Trata-se apenas de uma representação sobre a proporcionalidade direta entre as variáveis.

Logo, podemos ver como é o decaimento da densidade nesse modelo mais complexo. A queda é mais sutil e sensível às características de cada camada. No entanto, a diminuição ainda sim é observada e comprovada.

Código utilizado para a simulação:

import numpy as np
import matplotlib.pyplot as plt
import imageio
import os


# Constantes
g0 = 9.80665      # gravidade (m/s²)
R = 287.05        # constante dos gases (J/kg·K)


# Camadas: (base_alt, base_temp, lapse_rate, top_alt)
layers = [
   (0,     288.15, -0.0065, 11000),   # Troposfera
   (11000, 216.65,  0.0,    20000),   # Tropopausa
   (20000, 216.65,  0.001,  32000),   # Estratosfera 1
   (32000, 228.65,  0.0028, 47000),   # Estratosfera 2
   (47000, 270.65,  0.0,    51000),   # Estratopausa
   (51000, 270.65, -0.0028, 71000),   # Mesosfera 1
   (71000, 214.65, -0.002,  80000),   # Mesosfera 2
]


# Altitudes e inicialização
altitudes = np.linspace(0, 80000, 800)
densities = np.zeros_like(altitudes)


# Condições iniciais
p_base = 101325  # Pa
T_base = 288.15  # K
h_base = 0       # m


alt_pointer = 0


for layer in layers:
   h_start, T_start, L, h_end = layer
   while alt_pointer < len(altitudes) and altitudes[alt_pointer] <= h_end:
       h = altitudes[alt_pointer]
       delta_h = h - h_base


       if L == 0:
           # Isotérmica
           T = T_base
           p = p_base * np.exp(-g0 * delta_h / (R * T))
       else:
           # Com gradiente
           T = T_base + L * delta_h
           p = p_base * (T / T_base) ** (-g0 / (R * L))


       rho = p / (R * T)
       densities[alt_pointer] = rho
       alt_pointer += 1


   # Atualiza base para próxima camada
   h_base = h_end
   delta_h_layer = h_base - h_start
   if L == 0:
       T_base = T_start
       p_base = p_base * np.exp(-g0 * (h_base - h_start) / (R * T_base))
   else:
       T_base = T_start + L * delta_h_layer
       p_base = p_base * (T_base / T_start) ** (-g0 / (R * L))


# Camadas para exibir no gráfico
visual_layers = [
   {"name": "Troposfera",    "start": 0,   "end": 11,  "color": "#d0f0c0"},
   {"name": "Tropopausa",    "start": 11,  "end": 20,  "color": "#96dee9"},
   {"name": "Estratosfera",  "start": 20,  "end": 50,  "color": "#ffcc9a"},
   {"name": "Estratopausa",  "start": 50,  "end": 51,  "color": "#ff9a9a"},
   {"name": "Mesosfera",     "start": 51,  "end": 80,  "color": "#d79aff"},
]


# Criar diretório para os frames
if not os.path.exists("frames"):
   os.makedirs("frames")


# Criar os frames
for i in range(1, len(altitudes) + 1):
   plt.figure(figsize=(10, 6))


   for layer in visual_layers:
       plt.axvspan(layer["start"], layer["end"], color=layer["color"], alpha=0.8, zorder=1)
       plt.text(
           (layer["start"] + layer["end"]) / 2,
           0.3,
           layer["name"],
           ha='center',
           va='center',
           fontsize=9,
           rotation=90,
           color='black',
           zorder=2
       )


   plt.plot(altitudes[:i] / 1000, densities[:i], label='Densidade Atmosférica', color='black', zorder=10)
   plt.xlabel('Altitude (km)')
   plt.ylabel('Densidade (kg/m³)')
   plt.title('Densidade Atmosférica vs. Altitude (Modelo Estratificado)')
   plt.yscale('log')
   plt.grid(True, which='both', linestyle='--', linewidth=0.5)
   plt.ylim(1e-5, 1.5)
   plt.xlim(0, 80)
   plt.legend()


   filename = f"frames/frame_{i:03d}.png"
   plt.savefig(filename)
   plt.close()


# Criar GIF
frames = [imageio.imread(f"frames/frame_{i:03d}.png") for i in range(1, len(altitudes) + 1)]
imageio.mimsave("atmospheric_density_ni.gif", frames, fps=10)


print("GIF com modelo atmosférico corrigido criado com sucesso!")