📖 第17课：FHP模型

FHP 三角形网格 Navier-Stokes 各向同性

🔺 三角形网格——物理仿真的突破

1986年，Frisch、Hasslacher和Pomeau提出了FHP模型——使用三角形网格（六角形邻域）的Lattice Gas。这个看似简单的改变解决了HPP模型的各向异性问题，使得Lattice Gas能够真正模拟Navier-Stokes方程。

FHP模型的核心改进：

六角形网格（6个速度方向 + 静止粒子），自然各向同性
6个运动方向：间隔60°，等间距
FHP-I: 6速度（无静止粒子）
FHP-II: 7速度（含静止粒子）
FHP-III: 7速度 + 额外碰撞规则

六角形网格的对称群是D₆（六角对称群），足以保证宏观应力张量的各向同性！

FHP碰撞规则

二体碰撞（最基本）：

两个对头粒子碰撞后旋转60°

方向0+方向3 → 方向1+方向4 或方向2+方向5（两种结果，随机选择）

每种二体碰撞有2种可能结果，概率各1/2

三体碰撞：

三个间隔120°的粒子同时碰撞

方向0+方向2+方向4 → 方向1+方向3+方向5（旋转60°）

反之亦然

⚡ FHP模型的Verilog实现

FHP-II Lattice Gas引擎

// ============================================================================
// lattice_gas_fhp.v - FHP-II模型Lattice Gas引擎
// 7速度方向：6个运动方向(0-5,间隔60°) + 1个静止粒子(6)
// 三角形网格用"奇偶行偏移"编码
// ============================================================================
module lattice_gas_fhp #(
    parameter WIDTH  = 64,
    parameter HEIGHT = 64
)(
    input  wire              clk,
    input  wire              rst_n,
    input  wire              enable,
    input  wire              init,
    output wire [6:0]        grid_out [0:WIDTH*HEIGHT-1],
    output wire [31:0]       step_count,
    output wire [31:0]       density,      // 宏观密度
    output wire signed [31:0] velocity_x,  // 宏观速度x
    output wire signed [31:0] velocity_y   // 宏观速度y
);

    // ---- 网格存储：7位/格点 ----
    // bit[5:0] = 6个运动方向粒子
    // bit[6]   = 静止粒子
    // 方向编码: 0=右(0°) 1=右上(60°) 2=左上(120°)
    //           3=左(180°) 4=左下(240°) 5=右下(300°)
    reg [6:0] grid [0:WIDTH*HEIGHT-1];
    reg [6:0] grid_nxt [0:WIDTH*HEIGHT-1];
    reg [31:0] steps;

    // ---- 伪随机数（用于碰撞选择的随机性） ----
    reg [15:0] rng;
    always @(posedge clk or negedge rst_n) begin
        if (!rst_n) rng <= 16'hCAFE;
        else rng <= {rng[14:0], rng[15] ^ rng[14] ^ rng[12] ^ rng[3]};
    end

    // ---- 碰撞规则 ----
    function [6:0] collide;
        input [6:0] state;
        input rand_bit;
        begin
            collide = state;  // 默认不变
            // 二体碰撞：对头粒子
            case (state[5:0])
                6'b100001: collide = rand_bit ? 6'b010010 : 6'b001100;  // 0+3
                6'b010010: collide = rand_bit ? 6'b100001 : 6'b001100;  // 1+4
                6'b001100: collide = rand_bit ? 6'b100001 : 6'b010010;  // 2+5
                default: collide = state[5:0];
            endcase
            collide[6] = state[6];  // 静止粒子不变
        end
    endfunction

    // ---- 邻居计算（三角形网格） ----
    // 偶数行和奇数行偏移不同
    integer idx, x, y;
    reg is_odd_row;
    reg [6:0] after_col;

    always @(*) begin
        for (idx = 0; idx < WIDTH*HEIGHT; idx = idx + 1) begin
            x = idx % WIDTH;
            y = idx / HEIGHT;
            is_odd_row = y[0];

            // 步骤1：碰撞
            after_col = collide(grid[idx], rng[0]);

            // 步骤2：传播（根据方向移动到邻居）
            // 方向0(右): 从(x-1,y)接收
            // 方向1(右上): 从左下邻居接收
            // 方向2(左上): 从左下邻居接收
            // 方向3(左): 从(x+1,y)接收
            // 方向4(左下): 从右上邻居接收
            // 方向5(右下): 从右上邻居接收
            // 方向6(静止): 本格点

            // 简化实现：环形边界
            grid_nxt[idx][0] = grid[(y)*WIDTH + (x==0?WIDTH-1:x-1)][0];  // 右向来自左边
            grid_nxt[idx][3] = grid[(y)*WIDTH + (x==WIDTH-1?0:x+1)][3];  // 左向来自右边
            // 对角方向需要考虑奇偶行
            if (is_odd_row) begin
                grid_nxt[idx][1] = (y==0)?1'b0 : grid[(y-1)*WIDTH+x][1];
                grid_nxt[idx][2] = (y==0)?1'b0 : grid[(y-1)*WIDTH+(x==0?WIDTH-1:x-1)][2];
                grid_nxt[idx][4] = (y==HEIGHT-1)?1'b0 : grid[(y+1)*WIDTH+x][4];
                grid_nxt[idx][5] = (y==HEIGHT-1)?1'b0 : grid[(y+1)*WIDTH+(x==0?WIDTH-1:x-1)][5];
            end else begin
                grid_nxt[idx][1] = (y==0)?1'b0 : grid[(y-1)*WIDTH+(x==WIDTH-1?0:x+1)][1];
                grid_nxt[idx][2] = (y==0)?1'b0 : grid[(y-1)*WIDTH+x][2];
                grid_nxt[idx][4] = (y==HEIGHT-1)?1'b0 : grid[(y+1)*WIDTH+(x==WIDTH-1?0:x+1)][4];
                grid_nxt[idx][5] = (y==HEIGHT-1)?1'b0 : grid[(y+1)*WIDTH+x][5];
            end
            grid_nxt[idx][6] = after_col[6];  // 静止粒子不移动
        end
    end

    // 状态更新
    always @(posedge clk or negedge rst_n) begin
        if (!rst_n) begin
            steps <= 32'd0;
            for (idx = 0; idx < WIDTH*HEIGHT; idx = idx + 1)
                grid[idx] <= 7'd0;
        end else if (init) begin
            steps <= 32'd0;
            for (idx = 0; idx < WIDTH*HEIGHT; idx = idx + 1)
                grid[idx] <= 7'd0;
        end else if (enable) begin
            for (idx = 0; idx < WIDTH*HEIGHT; idx = idx + 1)
                grid[idx] <= grid_nxt[idx];
            steps <= steps + 32'd1;
        end
    end

    // 宏观量计算
    reg [31:0] dens; reg signed [31:0] vx, vy;
    always @(*) begin
        dens = 0; vx = 0; vy = 0;
        for (idx = 0; idx < WIDTH*HEIGHT; idx = idx + 1) begin
            dens = dens + grid[idx][0]+grid[idx][1]+grid[idx][2]+
                   grid[idx][3]+grid[idx][4]+grid[idx][5]+grid[idx][6];
            // x方向动量: 方向0(+) 方向3(-) 方向1,5(+0.5) 方向2,4(-0.5)
            vx = vx + grid[idx][0] - grid[idx][3] +
                  (grid[idx][1]+grid[idx][5])/2 - (grid[idx][2]+grid[idx][4])/2;
            // y方向动量: 方向1,2(+) 方向4,5(-)
            vy = vy + (grid[idx][1]+grid[idx][2]) - (grid[idx][4]+grid[idx][5]);
        end
    end

    genvar gi;
    generate
        for (gi = 0; gi < WIDTH*HEIGHT; gi = gi + 1)
            assign grid_out[gi] = grid[gi];
    endgenerate
    assign step_count = steps;
    assign density = dens;
    assign velocity_x = vx;
    assign velocity_y = vy;

endmodule

📐 从微观到宏观——统计力学推导

Chapman-Enskog展开：

1. 定义分布函数 f_i(x,t) = 格点(x,t)方向i有粒子的概率

2. 写出Boltzmann方程：f_i(t+1, x+c_i) = f_i(t,x) + Ω_i(f)

3. 其中 Ω_i 是碰撞算子

4. 对小Knudsen数展开：f = f⁰ + εf¹ + ε²f² + ...

5. 零阶：平衡分布 → Euler方程

6. 一阶：非平衡修正 → Navier-Stokes方程！

关键结果：FHP模型的宏观方程为

∂ρ/∂t + ∇·(ρu) = 0（质量守恒）

∂(ρu)/∂t + ∇·(ρuu) = -∇p + ν∇²u（Navier-Stokes！）

其中 ν = (1/12)(1/(2τ) - 1) 是运动粘度，τ与碰撞频率有关

🏋️ 练习

练习17.1：在FHP模型中模拟Poiseuille流——两个平行壁面之间的压力驱动流。测量中心线速度分布，与解析解对比。

练习17.2：模拟圆柱绕流——在均匀流中放置圆形障碍物，观察尾流是否形成卡门涡街。

练习17.3：测量FHP模型的运动粘度ν，通过观察扩散速率。与理论值对比。

练习17.4：比较HPP和FHP在同一初始条件下的行为差异。FHP是否解决了各向异性问题？

练习17.5（挑战）：实现Lattice Boltzmann Method（LBM）——FHP的"浮点数版本"，用浮点分布函数替代布尔粒子。分析精度提升。

🏆 成就解锁

🏅 流体力学师

你已经理解了Lattice Gas从微观到宏观的统计力学推导——简单碰撞规则如何涌现出Navier-Stokes方程！

✅ 实现FHP-II模型
✅ 理解三角形网格的各向同性
✅ 掌握碰撞-传播两步算法
✅ 了解Chapman-Enskog展开

📐 深入分析：FHP模型的粘度计算

运动粘度的理论推导：

FHP-I模型的运动粘度：

ν = (1/12) × (1/(2τ) - 1/2)

其中τ是平均碰撞间隔时间

对于稀薄气体（低密度极限）：τ → ∞, ν → 1/24

对于稠密气体（高密度）：τ减小，ν增大

粘度的硬件测量方法：

1. 初始化一个正弦速度场 v_x(y) = A·sin(ky)

2. 测量振幅衰减率 A(t) = A₀·exp(-νk²t)

3. 从衰减率反推ν

这个方法在硬件中可以直接实现——只需要一个正弦波初始条件和振幅检测器

补充实现

// 粘度测量模块 - 通过衰减率测量FHP的运动粘度
module viscosity_meter #(
    parameter WIDTH = 128,
    parameter HEIGHT = 128,
    parameter WAVELENGTH = 32
)(
    input  wire             clk,
    input  wire             rst_n,
    input  wire [6:0]       grid [0:WIDTH*HEIGHT-1],
    output wire [15:0]      viscosity,
    output wire [15:0]      amplitude,
    output wire             measurement_valid
);
    // 测量特定y坐标处的x方向动量
    localparam SAMPLE_Y = HEIGHT / 2;
    localparam K = 2 * 31416 / WAVELENGTH;  // 波数
    
    reg signed [31:0] vx_sum;
    integer x;
    always @(*) begin
        vx_sum = 32'd0;
        for (x = 0; x < WIDTH; x = x + 1) begin
            // x方向动量 = 方向0粒子数 - 方向3粒子数
            vx_sum = vx_sum + grid[SAMPLE_Y*WIDTH+x][0] - grid[SAMPLE_Y*WIDTH+x][3];
        end
    end
    
    assign amplitude = (vx_sum[31:16] < 0) ? -vx_sum[15:0] : vx_sum[15:0];
    assign viscosity = 16'd0;  // 需要多步测量计算
    assign measurement_valid = 1'b1;
endmodule

性能与优化分析

FHP的雷诺数限制：

雷诺数 Re = UL/ν，其中U是特征速度，L是特征长度

在FHP中，U最大为1（每步移动1格），L最大为网格宽度

对于128×128网格：Re_max ≈ 128/0.04 ≈ 3200

这足以观察到层流→湍流转捩！

更高的Re需要更大的网格或更低的粘度（更多碰撞规则）

📖 扩展阅读

Wolfram, S. A New Kind of Science, Wolfram Media, 2002.
Toffoli, T. & Margolus, N. Cellular Automata Machines, MIT Press, 1987.
Ilachinski, A. Cellular Automata: A Discrete Universe, World Scientific, 2001.

🔬 补充专题：实验方法论

在进行CA实验时，科学的方法论至关重要。以下是一些通用的实验指导原则：

实验设计三要素：

控制变量：每次只改变一个参数，保持其他不变
可重复性：记录所有参数，确保实验可以精确复现
统计显著性：多次运行取平均，避免偶然结果误导

参数扫描方法

系统性地扫描参数空间是理解CA行为的关键技术：

参数	范围	步长	测量指标
规则号	0-255	1	种群密度/周期
初始密度	0.1-0.9	0.1	收敛时间
网格大小	16-256	×2	有限尺寸效应
边界条件	环形/固定/镜像	离散	边界效应

💡 实验记录模板：每次实验记录以下信息——日期、规则号、初始条件、网格大小、边界条件、运行步数、关键观察、数据文件路径。这是科学研究的良好习惯，也使得结果可以追溯和验证。

数据可视化技巧

CA实验产生的数据通常是高维的（空间+时间+状态）。有效的可视化对于理解至关重要：

时空图：横轴=空间位置，纵轴=时间，颜色=状态——最常用的CA可视化方式
相图：x轴=参数A，y轴=参数B，颜色=测量指标——用于参数空间扫描
频谱图：对种群时间序列做FFT，识别周期性
直方图：统计状态/密度的分布

自相关分析：

时间自相关：R(τ) = ⟨n(t)·n(t+τ)⟩ / ⟨n²⟩

空间自相关：C(r) = ⟨n(x)·n(x+r)⟩ / ⟨n²⟩

如果R(τ)以周期T振荡 → 系统有周期T的行为

如果C(r)幂律衰减 → 系统处于临界状态

◀ 第16课 📜 目录第18课 ▶

元胞自动机课程 · 从Conway到Langton到Lattice Gas