傾度風のシミュレーション(定常状態への収束・発散を調べる)

keidohu-sim.py

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

# 定数の定義
dt = 1  # タイムステップ
num_steps = 200  # ステップ数
pressure_gradient_force = np.array([0, 1.1])  # 気圧傾度力の大きさ（仮定）
coriolis_parameter = 1e-1  # コリオリ力の係数（仮定）
friction_coefficient = 0.0  # 摩擦係数

# 初期条件
initial_velocity = np.array([10.0, 0.0])  # 浮動小数点数として初期化
velocity = initial_velocity.copy()

# 図の設定
rng = 15
fig, ax = plt.subplots()
ax.set_xlim(-rng, rng)
ax.set_ylim(-rng, rng)

# ベクトルの初期プロット
velocity_arrow = ax.quiver(0, 0, velocity[0], velocity[1], color="black", angles="xy", scale_units="xy", scale=1)
pgf_arrow = ax.quiver(0, 0, pressure_gradient_force[0], pressure_gradient_force[1], color="red", angles="xy", scale_units="xy", scale=1)
coriolis_arrow = ax.quiver(0, 0, 0, 0, color="blue", angles="xy", scale_units="xy", scale=1)
net_force_arrow = ax.quiver(0, 0, 0, 0, color="purple", angles="xy", scale_units="xy", scale=1)

# ステップ数表示用テキスト
step_text = ax.text(-9, 9, "", fontsize=12)


def update(step):
    global velocity
    if step == 0:
        velocity[:] = initial_velocity  # ループ開始時に初期条件にリセット

    coriolis_force = coriolis_parameter * np.array([velocity[1], -velocity[0]])  # 北半球のコリオリ力（右向き）
    friction_force = 0
    friction_force = -friction_coefficient * velocity
    net_force = pressure_gradient_force + coriolis_force + friction_force
    velocity += net_force * dt

    # ベクトルの更新
    velocity_arrow.set_UVC(velocity[0], velocity[1])
    coriolis_arrow.set_UVC(coriolis_force[0], coriolis_force[1])
    net_force_arrow.set_UVC(net_force[0], net_force[1])

    # ステップ数の更新
    step_text.set_text(f"Step: {step}")

    return velocity_arrow, pgf_arrow, coriolis_arrow, net_force_arrow, step_text


ani = animation.FuncAnimation(fig, update, frames=num_steps, interval=100, blit=True)

plt.show()

はじめから気圧傾度力とコリオリ力が釣り合っていれば、その状態を保つが、そうでないとき、定常状態に収束するのか......?
摩擦力を導入する(摩擦係数を適当に設定する)と、適当なところで釣り合いやすくなる。