编辑代码

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

# 模拟参数
NUM_PARTICLES = 50  # 颗粒数量
SIM_TIME = 10.0  # 总模拟时间
DT = 0.02  # 时间步长
WINDOW_SIZE = 10.0  # 模拟区域尺寸
AIR_VISCOSITY = 1.8e-5  # 空气动力粘度(Pa·s)
PARTICLE_DENSITY = 1500  # 颗粒密度(kg/m³)
AIR_DENSITY = 1.2  # 空气密度(kg/m³)


class Particle:
    def __init__(self):
        self.pos = np.random.rand(2) * WINDOW_SIZE  # 初始随机位置
        self.vel = np.zeros(2)  # 初始速度
        self.diameter = np.random.uniform(0.02, 0.05)  # 随机粒径
        self.mass = PARTICLE_DENSITY * np.pi * (self.diameter ** 3) / 6

    def apply_force(self, force):
        self.vel += force / self.mass * DT


class FlowField:
    def __init__(self):
        # 简单的二维流场（示例：中心区域风速较高）
        self.wind_speed = 2.0  # 基准风速(m/s)

    def get_velocity(self, pos):
        # 创建简单的流场（可修改为更复杂的模式）
        x, y = pos
        center = WINDOW_SIZE / 2
        factor = 1 - abs(x - center) / center  # 中心区域风速高
        return np.array([self.wind_speed * factor, 0])


def update(frame):
    global particles, flow
    plt.cla()

    # 更新每个粒子
    for p in particles:
        # 流场作用力
        air_vel = flow.get_velocity(p.pos)
        rel_vel = air_vel - p.vel

        # 计算阻力（简化模型）
        Re = np.linalg.norm(rel_vel) * p.diameter * AIR_DENSITY / AIR_VISCOSITY
        Cd = 24 / Re + 6 / (1 + np.sqrt(Re)) + 0.4  # 标准阻力系数公式
        drag_force = 0.5 * AIR_DENSITY * np.pi * (p.diameter ** 2) / 4 * Cd * rel_vel

        # 重力
        gravity = np.array([0, -9.81 * p.mass])

        # 合力应用
        p.apply_force(drag_force + gravity)

        # 更新位置
        p.pos += p.vel * DT

        # 边界处理
        p.pos = np.clip(p.pos, 0, WINDOW_SIZE)
        if p.pos[1] <= 0 or p.pos[1] >= WINDOW_SIZE:
            p.vel[1] *= -0.5  # 弹性碰撞

    # 绘制
    plt.scatter([p.pos[0] for p in particles],
                [p.pos[1] for p in particles],
                s=[p.diameter * 1000 for p in particles],
                c='blue', alpha=0.6)

    # 绘制流场箭头
    x = np.linspace(0, WINDOW_SIZE, 6)
    y = np.linspace(0, WINDOW_SIZE, 6)
    X, Y = np.meshgrid(x, y)
    U = np.array([flow.get_velocity([xi, yi])[0] for xi, yi in zip(X.flatten(), Y.flatten())])
    V = np.array([flow.get_velocity([xi, yi])[1] for xi, yi in zip(X.flatten(), Y.flatten())])
    plt.quiver(X, Y, U.reshape(X.shape), V.reshape(Y.shape), color='red')

    plt.xlim(0, WINDOW_SIZE)
    plt.ylim(0, WINDOW_SIZE)
    plt.title(f"Gas-Solid Flow Simulation (t={frame * DT:.2f}s)")
    plt.xlabel("X Position")
    plt.ylabel("Y Position")


# 初始化
particles = [Particle() for _ in range(NUM_PARTICLES)]
flow = FlowField()

# 创建动画
fig = plt.figure(figsize=(8, 6))
ani = FuncAnimation(fig, update, frames=int(SIM_TIME / DT), interval=50)
plt.show()