Cookbook / 工程 / 04
60 分钟数值教学模型已执行的网页阅读版
下载并在 VibeIt 中运行所有教程

组装并验证梁有限元模型

组装 Hermite 梁单元、施加固定端，并将节点挠度与反力和闭式参考比较。

工程 · 第 04 课🔋 离线核心分析⏱ 60 分钟数值教学模型📦 numpy · pandas · matplotlib · IPython

目标与运行方法

本课适合具备 Python 基础的本科生。先阅读每一步的方法与图注，再修改参数。AI 练习为可选部分，不需要 AI 账号即可完成核心课程。数据写在 notebook metadata 中，不需要另下载数据文件。请保持 notebook 已保存到当前工作文件夹。

学习任务是理解本课的研究问题、执行透明的计算、校验结果，并说明结论的边界。

环境与参数

In [1] · Python
from pathlib import Path
import base64, gzip, hashlib, json, platform
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from IPython.display import display, HTML
from IPython.display import display as _display
def display(value):
    if isinstance(value,pd.DataFrame):
        table=value.to_html(max_rows=12,escape=True)
        style="<style>.vibeit-readable-table{max-width:100%;overflow-x:auto}.vibeit-readable-table table{width:max-content;max-width:none;border-collapse:collapse}.vibeit-readable-table td,.vibeit-readable-table th{white-space:nowrap!important;word-break:normal!important;overflow-wrap:normal!important}</style>"
        return _display(HTML(style+'<div class="vibeit-readable-table">'+table+'</div>'))
    return _display(value)
RECIPE_ID='en04-beam-fem'
OUTPUT_DIR=Path.cwd()/"results"/RECIPE_ID
plt.rcParams.update({"figure.dpi":120,"font.size":10,"axes.spines.top":False,"axes.spines.right":False,
                      "axes.prop_cycle":plt.cycler(color=["#18785f","#416aa6","#bc6541","#9974a5"])})
print("Python:",platform.python_version(),"; NumPy:",np.__version__,"; Pandas:",pd.__version__)
print("Working folder:",Path.cwd())
print("Core data mode: embedded snapshot; no network requests")

Python: 3.13.13 ; NumPy: 2.5.3 ; Pandas: 3.0.6
Working folder: /Users/zhiluo/Developer/vibeit-cookbook/cookbook/downloads
Core data mode: embedded snapshot; no network requests

读取内嵌快照

In [2] · Python
def load_snapshot(recipe_id, expected_sha256, snapshot_key):
    """Read embedded data from this notebook, including a renamed copy."""
    for path in sorted(Path.cwd().glob("*.ipynb")):
        try:
            notebook = json.loads(path.read_text(encoding="utf-8"))
            metadata = notebook.get("metadata", {}).get("vibeit_cookbook", {})
            blob = metadata.get("snapshots", {}).get(snapshot_key)
            if metadata.get("recipe_id") != recipe_id or not blob:
                continue
            if blob.get("sha256") != expected_sha256:
                continue
            raw = gzip.decompress(base64.b64decode(blob["gzip_base64"], validate=True))
            if hashlib.sha256(raw).hexdigest() != expected_sha256:
                raise ValueError("Embedded snapshot checksum mismatch: " + snapshot_key)
            return json.loads(raw)
        except (OSError, json.JSONDecodeError):
            continue
    raise FileNotFoundError(
        "Save/open the complete .ipynb in the active working folder. "
        "The embedded snapshot could not be located for " + recipe_id)

SNAPSHOT_HASHES={'models': '48fd2e6e706650fe2ee0c0ad80b6d794a546a1dec78dd9d370ac09db3cc3adbf'}
def snapshot(key):
    return load_snapshot(RECIPE_ID,SNAPSHOT_HASHES[key],key)
print("Embedded snapshots:",list(SNAPSHOT_HASHES))

Embedded snapshots: ['models']

方法函数

In [3] · Python
def cantilever_fem(elements,length_m,young_Pa,second_moment_m4,tip_force_N):
    """Euler–Bernoulli Hermite beam elements: transverse displacement and rotation DOFs."""
    if not isinstance(elements,(int,np.integer)) or elements<1 or min(length_m,young_Pa,second_moment_m4)<=0:
        raise ValueError('Need positive beam geometry/material and an integer element count')
    nodes=elements+1;ell=length_m/elements;stiffness=np.zeros((2*nodes,2*nodes));load=np.zeros(2*nodes)
    local=young_Pa*second_moment_m4/ell**3*np.array([[12,6*ell,-12,6*ell],[6*ell,4*ell**2,-6*ell,2*ell**2],[-12,-6*ell,12,-6*ell],[6*ell,2*ell**2,-6*ell,4*ell**2]])
    for element in range(elements):
        indices=np.arange(2*element,2*element+4);stiffness[np.ix_(indices,indices)]+=local
    load[-2]=tip_force_N;displacement=np.zeros(2*nodes)
    displacement[2:]=np.linalg.solve(stiffness[2:,2:],load[2:])
    reaction=stiffness@displacement-load
    return dict(x=np.linspace(0,length_m,nodes),displacement_m=displacement[::2],rotation_rad=displacement[1::2],stiffness=stiffness,load=load,all_DOFs=displacement,reaction=reaction)

print("Method ready: cantilever_fem")

Method ready: cantilever_fem

数据与模型边界

建筑记录包含来源中单栋建筑的真实十分钟能耗及环境观测。机场天气字段含来源插值，不是同位置校准参考。梁、电路和 PID 例子是明确的数值设计模型；计算响应不能证明硬件性能或实际验收。

分步分析

1. 组装材料、几何及边界条件

每个 Hermite 梁节点有横向位移 m 与转角 rad。固定端将根部两自由度置零。名义 E、I 为教学参数，模型假定细长线弹性梁及小挠度。

In [4] · Python
p=snapshot("models")["beam"]
beam=cantilever_fem(10,p["length_m"],p["young_modulus_Pa"],p["second_moment_m4"],p["tip_force_N"])
dense=np.linspace(0,p["length_m"],501)
analytic=p["tip_force_N"]*dense**2*(3*p["length_m"]-dense)/(6*p["young_modulus_Pa"]*p["second_moment_m4"])
fig,ax=plt.subplots(figsize=(8,4));ax.plot(dense,1000*analytic,label="Closed tip-load solution")
ax.plot(beam["x"],1000*beam["displacement_m"],"o",label="Finite-element nodes")
ax.set(xlabel="Beam coordinate (m)",ylabel="Transverse deflection (mm)",title="Declared cantilever tip-load model");ax.legend();plt.tight_layout();plt.show()
Computed figure from this code cell

2. 检查组装后的连接结构

固定端约束去除前的刚度稀疏模式反映相邻单元耦合。位移与转角自由度单位不同，原始条件数依赖缩放，不是材料性质。

In [5] · Python
fig,ax=plt.subplots(figsize=(6,5));ax.spy(beam["stiffness"],markersize=3)
ax.set(xlabel="DOF column",ylabel="DOF row",title="Assembled beam stiffness connectivity");plt.tight_layout();plt.show()
Computed figure from this code cell

3. 检查此可精确表示载荷的不同网格

闭式挠度为三次函数，匹配 Hermite 基，因此节点端载解即便粗网格也可在舍入精度内正确。这是特殊正确性基准，不能证明任意载荷、几何或动力问题只需一单元。

In [6] · Python
records=[]
reference_tip=p["tip_force_N"]*p["length_m"]**3/(3*p["young_modulus_Pa"]*p["second_moment_m4"])
for count in [1,2,4,8,16]:
 trial=cantilever_fem(count,p["length_m"],p["young_modulus_Pa"],p["second_moment_m4"],p["tip_force_N"])
 records.append({"elements":count,"tip_deflection_m":trial["displacement_m"][-1],"reference_m":reference_tip,"relative_error":abs(trial["displacement_m"][-1]/reference_tip-1)})
mesh=pd.DataFrame(records);display(mesh)
fig,ax=plt.subplots(figsize=(7,4));ax.plot(mesh.elements,1000*mesh.tip_deflection_m,"o-")
ax.axhline(1000*reference_tip,color="gray",linestyle="--")
ax.set(xlabel="Elements",ylabel="Tip deflection (mm)",title="Special cubic benchmark across meshes");plt.tight_layout();plt.show()
elements tip_deflection_m reference_m relative_error
0 1 0.00162 0.00162 2.220446e-16
1 2 0.00162 0.00162 2.220446e-15
2 4 0.00162 0.00162 2.375877e-14
3 8 0.00162 0.00162 1.620926e-14
4 16 0.00162 0.00162 2.161826e-12
Computed figure from this code cell

4. 核对平衡反力并导出

按声明符号约定，根部反力平衡端力，根部反力矩平衡 F·L，内部残差应为零。这是模型平衡检查，不是加载试验验收。

In [7] · Python
assert np.allclose(beam["stiffness"],beam["stiffness"].T)
assert np.max(abs(beam["reaction"][2:]))<1e-6
assert np.isclose(beam["reaction"][0],-p["tip_force_N"],rtol=1e-8)
assert np.isclose(beam["reaction"][1],-p["tip_force_N"]*p["length_m"],rtol=1e-8)
assert mesh.relative_error.max()<1e-8
OUTPUT_DIR.mkdir(parents=True,exist_ok=True)
pd.DataFrame({"x_m":beam["x"],"deflection_m":beam["displacement_m"],"rotation_rad":beam["rotation_rad"]}).to_csv(OUTPUT_DIR/"beam-nodes.csv",index=False)
mesh.to_csv(OUTPUT_DIR/"cubic-benchmark-meshes.csv",index=False)
print("Assembly, reaction and closed-solution checks passed; exports:",OUTPUT_DIR)

Assembly, reaction and closed-solution checks passed; exports: /Users/zhiluo/Developer/vibeit-cookbook/cookbook/downloads/results/en04-beam-fem

自主练习与参考答案

  1. 为何三次端载基准看似不依赖网格？
  2. 位移与转角自由度为何不可互换？

参考解释。 该特殊载荷的闭式挠度为三次，Hermite 单元可表示。位移单位 m、转角 rad，刚度耦合不同物理量。

先保留当前 notebook 副本，再修改参数。把新图与原图一起比较，并记录改变了哪些输入、哪些计算和哪些结论。

可选的 AI 编程练习

提示词 1

仅用 NumPy、Matplotlib 增加解析转角，并比较全部节点转角。

提示词 2

使用同样的包分别改变 E、I，检查逆刚度缩放，保留力/力矩平衡。

人工验收。 固定端自由度为零，反力平衡 F 与 F·L，节点参考在数值容差内一致。

检查修改后的源代码，再手动运行受影响单元及后续校验。要求助手解释修改的目的、使用的包和数据来源，并用实际输出核对解释。

可选联网扩展

In [8] · Python
RUN_ONLINE_UPDATE=False
if RUN_ONLINE_UPDATE:
    import requests
    try:
        response=requests.get('https://docs.scipy.org/doc/scipy/tutorial/linalg.html',timeout=20)
        response.raise_for_status()
        print("Source:",response.url)
        print(response.text[:700])
    except requests.RequestException as error:
        print("Online extension unavailable:",type(error).__name__,str(error)[:160])
else:
    print("Online extension skipped; embedded core data unchanged")

Online extension skipped; embedded core data unchanged

来源、快照与下一步

  • IAPWS saturation reference release — IAPWS SR1-86(1992): unrestricted publication allowed in all countries (cover page); retrieved 2026-10-09. SHA-256 a10d38d73339b0713c86cfac88f9adcba137202167ccdfa1cd58238c438361d7. Reference equations and numerical coefficients transcribed from release pages 2–3; pedagogical model values are not new experiments.

  • Authored numerical teaching model parameters — Original teaching parameters/code MIT; reference equation coefficients attributed to IAPWS SR1-86(1992); retrieved 2026-10-09. SHA-256 4350c4305a83c772fa295a8be0ef86a247f2bfaae930a081720d3cd0c753abe0. Declared parameters and seed; generated trajectories and calibration observations must be labelled simulated.

下一步。 在保留本课的数据检查、参数记录和结果校验之后，将相同方法用于自己的公开或已获授权数据。记录研究设计与方法限制，不把示例结果当成新实验的验证。完整数据许可与生成说明见 cookbook 网页。