每日算法 — 使用java实现桥梁建造师:桁架有限元分析与受力平衡求解

桥梁建造师(Bridge Constructor)是一款把结构工程师的日常伪装成游戏的益智作品:预算有限,却要在峡谷上架一座能让卡车安全通过、自己又不被压塌的桥。玩家每次加一根钢梁,都要在”省材料”和”别塌”之间赌博。

这篇文章用 Java 实现它的计算内核:桁架结构有限元分析。核心是把整座桥离散成节点与杆件,组装全局刚度矩阵,解一个线性方程组得到每个节点的位移,再由位移反推每根杆件的轴力,最后与材料强度比较判定安全与破坏。

一、物理模型:为什么桁架可以只用轴力

真实的梁会受到弯矩,但桁架(truss)有一个简化假设:所有杆件只在两端受力,节点视为理想铰接,不传递弯矩。因此每根杆件只承受沿自身轴线方向的拉/压轴力(axial force),整座桥的计算就退化成纯粹的”节点力平衡”。

模型定义如下:

  • 节点:带 2 个自由度(x、y 位移),承重物固定在若干节点上
  • 杆件:连接两个节点,材料参数为杨氏模量 E、截面积 A、长度 L
  • 边界条件:桥两端锚固在地面,位移为 0

杆件轴向刚度是关键量:

k = E * A / L

它决定了”拉长单位长度需要多少力”。刚度越大越不形变,成本也越高——这正是游戏里钢材与木材的核心差异。

二、单元刚度矩阵与坐标变换

一根与水平方向夹角为 θ 的杆件,其局部刚度矩阵(沿轴向)需要旋转到全局 x-y 坐标系。标准结果是 4×4 矩阵:

K_e = k * [ c²   cs  -c²  -cs ]
          [ cs   s²  -cs  -s² ]
          [-c²  -cs   c²   cs ]
          [-cs  -s²   cs   s² ]

其中 c = cos θ = dx/L,s = sin θ = dy/L。

public final class Node {
    public final int id;
    public double x, y;        // 初始坐标(几何不变)
    public double dx, dy;      // 待求位移
    public double fx, fy;      // 外力(含自重与车辆)
    public boolean fixed;      // 锚点,位移强制为 0

    public Node(int id, double x, double y, boolean fixed) {
        this.id = id; this.x = x; this.y = y; this.fixed = fixed;
    }
}

public final class Member {
    public final Node a, b;
    public final double E;     // 杨氏模量
    public final double A;     // 截面积
    public final double length;
    public double force;       // 求解后得到的轴力:正为拉,负为压
    public final double maxTension;    // 抗拉极限
    public final double maxCompression;// 抗压极限

    public Member(Node a, Node b, double E, double A,
                  double maxTension, double maxCompression) {
        this.a = a; this.b = b; this.E = E; this.A = A;
        this.length = Math.hypot(b.x - a.x, b.y - a.y);
        this.maxTension = maxTension;
        this.maxCompression = maxCompression;
    }

    public double stiffness() {
        return E * A / length;   // 轴向刚度 k
    }

    /** 方向余弦:把轴向变形投影到 x / y */
    public double[] directionCosines() {
        return new double[]{ (b.x - a.x) / length, (b.y - a.y) / length };
    }
}

组装时按节点编号映射到全局矩阵的下标:节点 i 占用行/列 2i 与 2i+1。

public final class GlobalAssembler {
    /** 返回 2n x 2n 的全局刚度矩阵 */
    public static double[][] assemble(java.util.List<Node> nodes,
                                      java.util.List<Member> members) {
        int n = nodes.size() * 2;
        double[][] K = new double[n][n];

        for (Member m : members) {
            double k = m.stiffness();
            double[] cs = m.directionCosines();
            double c = cs[0], s = cs[1];
            double c2 = c * c, s2 = s * s, cs2 = c * s;
            double[][] ke = {
                { c2,  cs2, -c2, -cs2 },
                { cs2,  s2, -cs2, -s2 },
                { -c2, -cs2, c2,  cs2 },
                { -cs2, -s2, cs2,  s2 }
            };
            int[] idx = { 2 * m.a.id, 2 * m.a.id + 1, 2 * m.b.id, 2 * m.b.id + 1 };
            for (int i = 0; i < 4; i++) {
                for (int j = 0; j < 4; j++) {
                    K[idx[i]][idx[j]] += k * ke[i][j];
                }
            }
        }
        return K;
    }
}

由于 K 是对称半正定矩阵,锚点节点对应行列会造成奇异。标准做法是”消元法处理边界条件”:把固定自由度对应的行与列删掉(或把行列置为单位、右端置零),只解剩余自由节点。

三、核心求解:带主元的高斯消元

方程组形式为 K · d = f,其中 d 是待求位移向量,f 是节点外力向量(自重 + 车辆载荷)。桁架规模通常几十到上百个自由度,直接用部分主元高斯消元即可,O(n³) 在游戏规模下毫无压力。

public final class LinearSolver {
    /** 部分主元高斯消元,原地求解 A x = b;A 会被修改 */
    public static double[] solve(double[][] A, double[] b) {
        int n = b.length;
        for (int col = 0; col < n; col++) {
            // 选主元:取该列绝对值最大的行,避免除以极小数
            int pivot = col;
            for (int r = col + 1; r < n; r++) {
                if (Math.abs(A[r][col]) > Math.abs(A[pivot][col])) pivot = r;
            }
            if (Math.abs(A[pivot][col]) < 1e-12) {
                throw new IllegalStateException("刚度矩阵奇异:结构存在机构自由度(未约束)");
            }
            if (pivot != col) {
                double[] tmp = A[pivot]; A[pivot] = A[col]; A[col] = tmp;
                double t = b[pivot]; b[pivot] = b[col]; b[col] = t;
            }
            // 消去下方各行
            for (int r = col + 1; r < n; r++) {
                double factor = A[r][col] / A[col][col];
                if (factor == 0) continue;
                for (int c = col; c < n; c++) A[r][c] -= factor * A[col][c];
                b[r] -= factor * b[col];
            }
        }
        // 回代
        double[] x = new double[n];
        for (int i = n - 1; i >= 0; i--) {
            double sum = b[i];
            for (int j = i + 1; j < n; j++) sum -= A[i][j] * x[j];
            x[i] = sum / A[i][i];
        }
        return x;
    }
}

那句”刚度矩阵奇异”的异常其实是游戏里最有价值的信息:它说明结构存在机构自由度(mechanism),也就是某一部分可以自由摆动而没有杆件抵抗。桥梁建造师里表现为”桥面一晃就散”,本质是缺少三角形支撑。

四、从位移反推轴力与安全判定

位移求出后,杆件的轴向变形量是两端沿杆方向的相对位移投影:

Δ = (d_b - d_a) · n
N = k · Δ

N > 0 为拉力,N < 0 为压力。判据是 |N| 与材料极限比较。

public final class TrussAnalyzer {
    /** 求解整座桥,返回是否安全 */
    public static boolean analyze(java.util.List<Node> nodes,
                                  java.util.List<Member> members) {

        double[][] K = GlobalAssembler.assemble(nodes, members);

        // 收集自由自由度,建立"全局下标 -> 压缩下标"映射
        int[] map = new int[nodes.size() * 2];
        java.util.Arrays.fill(map, -1);
        int free = 0;
        for (int i = 0; i < nodes.size(); i++) {
            Node nd = nodes.get(i);
            if (!nd.fixed) { map[2 * i] = free++; map[2 * i + 1] = free++; }
        }

        double[][] Kf = new double[free][free];
        double[] ff = new double[free];
        for (int i = 0; i < nodes.size() * 2; i++) {
            if (map[i] < 0) continue;
            ff[map[i]] = (i % 2 == 0) ? nodes.get(i / 2).fx : nodes.get(i / 2).fy;
            for (int j = 0; j < nodes.size() * 2; j++) {
                if (map[j] < 0) continue;
                Kf[map[i]][map[j]] = K[i][j];
            }
        }

        double[] u = LinearSolver.solve(Kf, ff);

        // 回填位移
        for (int i = 0; i < nodes.size(); i++) {
            Node nd = nodes.get(i);
            if (map[2 * i] >= 0) {
                nd.dx = u[map[2 * i]];
                nd.dy = u[map[2 * i + 1]];
            } else {
                nd.dx = 0; nd.dy = 0;
            }
        }

        // 计算每根杆件轴力
        boolean safe = true;
        for (Member m : members) {
            double[] cs = m.directionCosines();
            double c = cs[0], s = cs[1];
            double delta = (m.b.dx - m.a.dx) * c + (m.b.dy - m.a.dy) * s;
            m.force = m.stiffness() * delta;

            if (m.force > m.maxTension) {
                System.out.printf("杆件 #%d-%d 拉断:轴力 %.1f > 极限 %.1f%n",
                        m.a.id, m.b.id, m.force, m.maxTension);
                safe = false;
            } else if (-m.force > m.maxCompression) {
                System.out.printf("杆件 #%d-%d 压溃:轴力 %.1f > 极限 %.1f%n",
                        m.a.id, m.b.id, -m.force, m.maxCompression);
                safe = false;
            }
        }
        return safe;
    }
}

注意压力的极限通常远小于拉力:细长杆受压会屈曲(buckling),而不是单纯的材料强度问题。用欧拉屈曲公式估算临界力更贴近物理:

P_cr = π² E I / L²

I 是截面惯性矩。这就是桥梁建造师里”长杆做压杆必失败”的根源——L² 在分母,长度翻倍临界力降到 1/4。

五、优化:动载荷与材料预算

游戏的真正难度在于”卡车移动”。车辆不是静态载荷,而是一个沿桥面移动的集中力。做法是对位置做采样,逐个工况求解,取最坏情况。

public final class LoadCaseRunner {
    /**
     * 让重为 truckWeight 的车辆沿桥面节点序列移动,
     * 每个位置作为一个工况求解,返回全过程中最危险杆件的安全系数
     */
    public static double worstCaseSafety(java.util.List<Node> nodes,
                                         java.util.List<Member> members,
                                         java.util.List<Node> deck,
                                         double truckWeight) {
        double worst = Double.MAX_VALUE;

        for (int pos = 0; pos < deck.size(); pos++) {
            // 清空外力
            for (Node nd : nodes) { nd.fx = 0; nd.fy = 0; }

            // 自重:每个非锚点节点承受自身重力(简化为均布)
            for (Node nd : nodes) if (!nd.fixed) nd.fy += 10.0;

            // 车轮载荷分配到相邻两个桥面节点(线性插值)
            Node cur = deck.get(pos);
            cur.fy += truckWeight * 0.5;
            if (pos + 1 < deck.size()) deck.get(pos + 1).fy += truckWeight * 0.5;
            else cur.fy += truckWeight * 0.5;

            TrussAnalyzer.analyze(nodes, members);

            for (Member m : members) {
                double used = (m.force >= 0)
                        ? m.force / m.maxTension
                        : -m.force / m.maxCompression;
                double safety = 1.0 - used;   // 越小越危险
                worst = Math.min(worst, safety);
            }
        }
        return worst;
    }
}

复杂度上,若桥面有 P 个位置、自由度 N,总代价 O(P·N³)。P=20、N=100 时约 2×10⁷ 次浮点运算,Java 单线程可以在几十毫秒内跑完,完全可以做成”玩家点击测试”时实时反馈。

如果要进一步提速,可以利用刚度矩阵在几何固定时不变化这一点——预先做一次 LU 分解,每个载荷工况只需 O(N²) 的前代回代。这是结构分析里最经典的优化,能把整个测试从 O(P·N³) 降到 O(N³ + P·N²)。

六、建桥策略的算法视角

把求解器跑通之后,剩下的就是”如何在预算内建出安全桥”。这其实是一个带约束的拓扑优化问题,几个可落地的启发式:

  1. 先搭三角形再补横杆。静定桁架(三角形递推)不会因单根杆件失效而整体崩溃,算法上的原因是刚度矩阵满秩且无机构自由度。
  2. 长杆尽量做成受拉。受拉不受屈曲限制,同样长度下承载力高一个数量级。这是”下承式拱桥”和”斜拉”结构的算法依据。
  3. 把锚点分布开。锚点少意味着反力集中,锚固节点附近的杆件轴力会爆表;分散锚点相当于把边界反力摊平。
  4. 用安全系数热力图指导删杆。跑完全部工况后,每根杆件的利用率就是它的”价值”指标,接近 0 的杆件可以删掉换预算,这是最朴素的拓扑优化。

桥梁建造师的乐趣在于”用直觉猜结构”,而有限元分析把这份直觉变成了可计算的数字:刚度矩阵告诉你哪里软,轴力告诉你哪里危险,屈曲公式告诉你为什么细长杆不能受压。理解这三件事之后,你在游戏里省下的每一块钢板,都是算出来的。