1956：有限元素法 FEM — 把飛機翅膀切成小碎片#

副標題：剛度矩陣拼起來，複雜形狀也能算。

案發現場#

1950 年代波音工程師盯著三角翼發愁。噴射客機要更快更輕，翅膀應力卻算不出來。彈性力學方程漂亮，解只存在於矩形與圓形。真實翅膀變厚度、帶孔、有加強肋，解析解投降。有限差分要規則網格，遇斜邊圓孔就歪七扭八。被害人是 707 研發時程：算不出應力，不敢試飛。主角是 Turner、Clough、Martin、Topp 四位工程師。他們只想讓工程師在電腦上算出答案。 1956 年論文沒有「有限元素」一詞，卻藏著全部真相。 Clough 1960 年才命名 FEM ，但 1956 年是案發年。

偵查過程#

線索一：分片多項式的逆向思維#

傳統思路找全域函數滿足微分方程，難如登天。 FEM 反其道：切成小三角形，每片用簡單多項式。三角內部位移設線性函數，係數由三節點位移決定。形函數在己節點取一，別節點取零，支集只在相鄰單元。圓孔用小三角形逼近，複雜形狀就此被馴服。連續性只要求邊上位移連接，應變可逐片常數。工程師直覺戰勝數學家完美主義。

線索二：剛度矩陣的拼裝術#

每三角形單元有自己的小剛度矩陣 \(k_e\) 。它連起節點力與位移，由虛功原理在小域上積分得到。全域剛度 \(K\) 由所有 \(k_e\) 按節點編號組裝。規則是對號入座：共享節點貢獻直接相加。最終得全域系統，行內可記為 \(Ku = f\) 的形式。其中 \(K\) 為剛度矩陣， \(u\) 為位移， \(f\) 為載荷。矩陣 \(K\) 稀疏對稱正定，迭代法正好大展身手。固定節點划行划列，邊界處理直白。

線索三：1956 論文的證據鏈#

論文處理平面應力，三角單元大獲成功。懸臂樑與帶孔板計算對上實驗，誤差工程可忍。 | 單元數 | 自由度 | 最大應力誤差 | |--------|--------|--------------| | 粗網格 | 數十 | 約一成 | | 中網格 | 數百 | 數個百分點 | | 細網格 | 上千 | 工程可用 | 審稿人質疑變分基礎，工程師用數據反駁。數學家後補 Sobolev 空間與誤差估計。 NASA 資助的 NASTRAN 於 1968 年問世，FEM 推向世界。弱形式只要求加權積分成立，Galerkin 代入即得矩陣方程。

結案報告#

1956 年 FEM 破解複雜形狀偏微分方程之謎。連續體切碎再用矩陣拼回，化無限維為有限維。飛機、反應爐、橋樑從此有了數值保證。 NASTRAN、ANSYS、ABAQUS 皆本案直系子孫。它與有限體積、無網格法三分天下，仍是固體模擬之王。

證據與工具#

全域系統獨立公式如下，為核心證物。

\[ Ku = f \]

其中 \(K\) 為全域剛度， \(u\) 為位移向量， \(f\) 為節點力。單元剛度獨立公式如下，積分只在小單元上。

\[ k_e = \int_{\Omega_e} B^T D B \, d\Omega \]

其中 \(B\) 為應變位移矩陣， \(D\) 為本構矩陣。 | 符號 | 意義 | 性質 | |------|------|------| | \(K\) | 全域剛度 | 稀疏對稱正定 | | \(u\) | 節點位移 | 待求未知數 | | \(f\) | 節點力 | 含體力與邊界力 | 一維桿練習：切 \(n\) 段組三對角 \(K\) ，一端固定一端加力。加密網格，觀察端部位移收斂到解析解。因果鏈：傅立葉連續體方程在此被切碎，1977 年 SPH 另闢無網格路。數學家後補先驗誤差估計，收斂階數與多項式次數掛鉤。自適應網格在應力集中處自動加密，省下大把自由度。本案檔案編號 FEM-1956 ，狀態已結案，網格仍在細化。等參單元把曲邊映射到母單元，圓孔逼近更光滑。商業軟體 NASTRAN 1968 年問世，工程界全面換裝。

補充：程式實作#

對應程式：1956-fem_1d_poisson.py

理論呼應#

本文斷言連續體可切碎再拼回，化無限維偏微分方程為有限維矩陣方程。

程式以一維卜瓦松方程示範線性元，單元剛度按節點編號組裝為全域三對角矩陣。

全域系統正是正文證物 \(K u = f\) ，其中 \(K\) 為剛度而 \(u\) 為位移向量。

解析解為 \(x(1-x) / 2\) ，加密網格即見數值解逐點靠攏，正好呼應收斂表格。

碎片早已按編號歸位，剛度矩陣親自指認兇手就是離散誤差，案情就此收網。

如何執行#

在本文所屬目錄下執行以下指令即可重現全部數字。

python3 _code/1956-fem_1d_poisson.py

程式只依賴 numpy ，已固定亂數種子，任何機器重跑都應一致。

預設切分為 N = 50 段，邊界施加齊次狄利克雷條件，載荷為常數一。

實測結果抄錄#

以下數字為本次實際執行所得，未經修飾，直接抄錄關鍵行。

N=50, h=0.0200

max_err=5.412337e-16

VERIFY: max_err=0.0000000000

驗證行顯示最大誤差遠小於千分之一，程式最後印出 ALL CHECKS PASS 。

由此可見一維線性元在此光滑問題上幾乎重現了解析解的二次曲線。

動手實驗#

讀者可把 N 改為 10 或 200 ，觀察最大誤差如何隨網格加密而二次下降。

完整程式如下：

# -*- coding: utf-8 -*-
"""1956 一維 Poisson P1 有限元素 (對應 wiki: 計算模擬學 / 有限元素法)
只用 numpy,固定種子。-u''=1,u(0)=u(1)=0,N=50 線性元。
驗證:與解析解 x(1-x)/2 的最大誤差 < 1e-3。
"""
import numpy as np

np.random.seed(0)
N = 50
h = 1.0 / N
n_in = N - 1  # 內點數

# 剛度矩陣 (1/h)*tridiag(-1,2,-1),載重全為 h (f=1)
diag = np.full(n_in, 2.0 / h)
off = np.full(n_in - 1, -1.0 / h)
A = np.diag(diag) + np.diag(off, 1) + np.diag(off, -1)
F = np.full(n_in, h * 1.0)

u_in = np.linalg.solve(A, F)
x_in = np.linspace(h, 1.0 - h, n_in)
u_exact = x_in * (1 - x_in) / 2.0
err = float(np.max(np.abs(u_in - u_exact)))
ok = err < 1e-3
print(f"N={N}, h={h:.4f}, max_err={err:.6e} (驗證 <1e-3: {'PASS' if ok else 'FAIL'})")
print(f"VERIFY: max_err={err:.10f}")
assert ok, "FEM 誤差過大"
print("ALL CHECKS PASS")