当板壳理论遇上Python，这个题目听起来多少有点标题党，但我今天想聊的确实是件特别“手撕”的活：不依赖ANSYS的APDL分层模拟，不靠ABAQUS的UMAT子程序，只用Python把功能梯度板（FGM板）的自由振动分析从头到尾写出来。这些年做力学数值分析，我见过太多人卡在同一个地方——FGM材料的弹性模量和密度沿厚度方向连续变化，而商业软件偏偏只擅长按“层”定义材料，你只能把板剖成几十层给每层单独赋材料，建模烦、收敛慢、改个梯度指数又得重来。与其这样，不如把板壳理论推到前台，用Python直接求解。
如果你正准备做FGM板的模态分析、参数扫描，或者你只是个想拿力学问题练手Python数值方法的工程师、研究生，这篇文章应该能给你一条完整通路：先聊材料模型和板理论怎么选，再推导FSDT五自由度控制方程，然后借Navier解写一个几十行代码就能跑的模态求解器，最后用DQ法把边界条件从简支扩展到固支、悬臂。代码我会给出核心部分，所有结论都能自己复现。
先说个重要前提：整篇文章只依赖numpy、scipy、matplotlib三个库，Python 3.9到3.12实测都能跑。环境卡脖子的同学，终端敲一行pip install numpy scipy matplotlib就能开工。
1. 功能梯度板在算振动前，得先把“梯度”这件事说清楚
1.1 功能梯度材料的三张“配方表”：P-FGM、S-FGM与Mori-Tanaka
功能梯度材料的概念并不复杂：某种陶瓷和某种金属，在厚度方向上从一侧连续过渡到另一侧。拿最常见的ZrO2/Al体系举例，顶部是纯陶瓷（耐高温、刚度高），底部是纯金属（韧性好、抗断裂），中间没有明显界面。这个“没有明显界面”是FGM区别于传统层合板的核心，也是数值模拟的难点——材料参数随厚度坐标z连续变化，建模时必须把这种连续变化写进刚度积分里。
描述变化规律的“配方”主要有三类，我最常用的是P-FGM（幂律分布）：
[ E(z) = (E_c - E_m)\left(\frac{z}{h} + \frac12\right)^p + E_m ]
其中 (p) 是梯度指数，(E_c) 是陶瓷的弹性模量，(E_m) 是金属的模量，z从-h/2到h/2。p=0时整块板全是陶瓷，p趋向无穷大时趋近金属，p越大材料越“软”、越偏金属。密度(\rho(z))按同样形式变化，泊松比在大部分文献里假设为常数，这也是本文采用的近似。
S-FGM（Sigmoid分布）适合“上下表面都是陶瓷、中间是金属”的夹心式梯度设计，它能避免幂律在中面附近的突变；Mori-Tanaka模型则基于细观力学，在陶瓷体积分数较高时对等效剪切模量的预测更准。对振动频率而言，P-FGM是性价比最高的选择，原因很实际：公式简单、物理趋势直观、文献对比数据多。S-FGM和Mori-Tanaka的差别主要体现在定量数值上，定性规律一致。
1.2 板理论选择的AB面：CPT、FSDT与TSDT到底差在哪
板理论本质上是对三维弹性问题做降维处理。经典薄板理论（CPT）基于Kirchhoff假设：中面法线变形后仍垂直于中面，横向剪切变形被直接忽略。这个假设在薄板里成立，但一到中厚板就会把结构算“硬”，频率偏高。一阶剪切变形理论（FSDT）做了修正：法线在变形后不再垂直中面，引入两个独立的转角自由度，允许横向剪切。
FSDT相比高阶理论（TSDT，比如Reddy三阶理论）的实现成本低不少，精度对工程初步分析足够。三阶理论因为要满足上下表面剪应力为零的条件，位移场复杂度显著上升，刚度矩阵和惯性矩阵会多出一批高阶项，调试成本翻倍。我的建议是：除非你的板特别厚（a/h小于5）或者需要精确的层间应力，否则FSDT是手写代码的第一选择。下表是它们的粗略对比：
| 理论 | 自由度 | 适用厚跨比 | 剪切变形 | 代码难度 |
|---|---|---|---|---|
| CPT | 1 (w) | a/h > 20 | 忽略 | 低 |
| FSDT | 5 (u,v,w,φx,φy) | a/h > 5 | 常剪切修正 | 中 |
| TSDT | 5+高阶项 | a/h任意 | 高阶精确 | 高 |
1.3 为什么FGM板比均质板更需要关注剪切变形
这一点容易被忽略。FGM板是两种材料“混”出来的，金属一侧剪切模量低，而陶瓷一侧刚度高，整个截面的等效剪切刚度其实比按算术平均估计的更低。厚度越厚、陶瓷金属模量比越大，剪切效应越明显。我做参数扫描时发现，当a/h从20降到5，FSDT和CPT算出的第一阶频率差异可以超过10%，高阶模态差异还会更大。如果你用CPT去算中厚FGM板，得到“漂亮但偏危险”的频率值，这在工程上是不可接受的。
2. 从虚功原理到五自由度矩阵：FSDT控制方程的推导与离散
2.1 位移场假设与应变-几何关系
FSDT的位移场是：
[ \begin{aligned} u(x,y,z) &= u_0(x,y) + z\phi_x(x,y) \ v(x,y,z) &= v_0(x,y) + z\phi_y(x,y) \ w(x,y,z) &= w_0(x,y) \end{aligned} ]
这里的 (u_0、v_0、w_0) 是中面的三个位移，(\phi_x、\phi_y) 是中面法线的转角。横向剪切应变的表达式是：
[ \gamma_{xz} = w_{,x} + \phi_x, \quad \gamma_{yz} = w_{,y} + \phi_y ]
注意这个表述里，(\phi_x) 不是单纯的法线转角，而是包含了剪切角度的量。这也是FSDT和经典薄板理论最本质的差别。
2.2 本构关系与厚度方向的梯度积分
对各向同性材料，平面应力本构为：
[ \begin{bmatrix} \sigma_x \ \sigma_y \ \tau_{xy} \end{bmatrix} = \frac{E(z)}{1-\nu^2} \begin{bmatrix} 1 & \nu & 0 \ \nu & 1 & 0 \ 0 & 0 & \frac{1-\nu}{2} \end{bmatrix} \begin{bmatrix} \varepsilon_x \ \varepsilon_y \ \gamma_{xy} \end{bmatrix} ]
把位移场代入几何方程，再沿厚度积分，就能得到广义力与广义应变的关系：
[ \begin{bmatrix} \mathbf{N} \ \mathbf{M} \end{bmatrix}
\begin{bmatrix} \mathbf{A} & \mathbf{B} \ \mathbf{B} & \mathbf{D} \end{bmatrix} \begin{bmatrix} \boldsymbol{\varepsilon}^0 \ \boldsymbol{\kappa} \end{bmatrix} ]
其中 (\mathbf{A}, \mathbf{B}, \mathbf{D}) 分别是拉伸、拉弯耦合和弯曲刚度，它们的元素全都是对厚度方向积分：
[ (A_{ij}, B_{ij}, D_{ij}) = \int_{-h/2}^{h/2} Q_{ij}(z) (1, z, z^2) dz ]
横向剪切项是：
[ Q_x = \kappa_s A_s \gamma_{xz}, \quad Q_y = \kappa_s A_s \gamma_{yz}, \quad A_s = \int_{-h/2}^{h/2} G(z) dz ]
这里的 (\kappa_s) 是剪切修正系数，一般取5/6，但我会在后面专门指出它在FGM板里的坑。积分计算我的经验是用Gauss-Legendre公式，50个高斯点已经能让刚度系数收敛到小数点后六位以上。
2.3 运动方程与Navier级数代入
FSDT的五个运动方程是从哈密顿原理推来的，形式上是三个平动方程加两个转动方程。对四边简支（SSSS）矩形板，可以用Navier法获得解析解：把五个广义位移展开成双三角级数：
[ \begin{aligned} u_0 &= U_{mn} \cos(\alpha x) \sin(\beta y) \ v_0 &= V_{mn} \sin(\alpha x) \cos(\beta y) \ w_0 &= W_{mn} \sin(\alpha x) \sin(\beta y) \ \phi_x &= X_{mn} \cos(\alpha x) \sin(\beta y) \ \phi_y &= Y_{mn} \sin(\alpha x) \cos(\beta y) \end{aligned} ]
其中 (\alpha = m\pi/a, \beta = n\pi/b)。这套展开天然满足简支边界：(w=0)，(M_x=M_y=0)。代入运动方程后，偏微分方程就退化成代数特征值问题，得到一个5×5的对称矩阵。
2.4 从连续方程到广义特征值问题
最终的五自由度方程是：
[ \mathbf{K}{mn} \mathbf{d}{mn} = \omega^2 \mathbf{M}{mn} \mathbf{d}{mn} ]
其中 (\mathbf{d}_{mn} = [U, V, W, X, Y]^T)。刚度矩阵元素可以通过代换直接写出，例如：
[ \begin{aligned} K_{11} &= A_{11}\alpha^2 + A_{66}\beta^2 \ K_{12} &= (A_{12}+A_{66})\alpha\beta \ K_{14} &= B_{11}\alpha^2 + B_{66}\beta^2 \ K_{15} &= (B_{12}+B_{66})\alpha\beta \ K_{33} &= A_s\alpha^2 + A_s\beta^2 \ K_{34} &= A_s\alpha \ K_{44} &= D_{11}\alpha^2 + D_{66}\beta^2 + A_s \end{aligned} ]
惯性矩阵由 ((I_0, I_1, I_2)) 组成，其中 (I_1 = \int \rho(z) z dz) 在FGM板里通常不为零，它对应着面内位移和转动的惯性耦合。很多简化代码把(I_1)顺手设成0，这在梯度材料里会引入不可忽视的误差。
3. Navier解落地：用Python写出简支FGM板模态求解器
3.1 代码骨架与数据流
整个求解器我按四步组织：材料模块（计算梯度分布）、积分模块（计算A/B/D/As和I0/I1/I2）、组装模块（构造5×5矩阵）、求解模块（eigh求解特征值）。这样无论后面做参数扫描还是换积分规则，动一个模块就行。
3.2 核心实现：材料梯度与厚度积分
import numpy as np from scipy.linalg import eigh # 材料与几何参数（ZrO2/Al体系） E_c, nu_c, rho_c = 151e9, 0.3, 3000.0 # 陶瓷 E_m, nu_m, rho_m = 70e9, 0.3, 2707.0 # 金属 p = 1.0 # 梯度指数 a, b, h = 1.0, 1.0, 0.1 # 矩形板边长与厚度 kappa_s = 5.0 / 6.0 # 剪切修正系数（均质板默认值） def E_z(z): return (E_c - E_m) * (z/h + 0.5)**p + E_m def rho_z(z): return (rho_c - rho_m) * (z/h + 0.5)**p + rho_m # 厚度方向 Gauss-Legendre 积分 N_GP = 50 xi, wi = np.polynomial.legendre.leggauss(N_GP) z_pts = h / 2 * xi E_pts = E_z(z_pts) rho_pts = rho_z(z_pts) G_pts = E_pts / (2 * (1 + nu_c)) # 这里假设 nu_c = nu_m Q11 = E_pts / (1 - nu_c**2) Q12 = nu_c * E_pts / (1 - nu_c**2) Q66 = G_pts scale = h / 2 A11 = scale * np.sum(Q11 * wi) A12 = scale * np.sum(Q12 * wi) A66 = scale * np.sum(Q66 * wi) B11 = scale * np.sum(Q11 * z_pts * wi) B12 = scale * np.sum(Q12 * z_pts * wi) B66 = scale * np.sum(Q66 * z_pts * wi) D11 = scale * np.sum(Q11 * z_pts**2 * wi) D12 = scale * np.sum(Q12 * z_pts**2 * wi) D66 = scale * np.sum(Q66 * z_pts**2 * wi) As = kappa_s * scale * np.sum(G_pts * wi) I0 = scale * np.sum(rho_pts * wi) I1 = scale * np.sum(rho_pts * z_pts * wi) I2 = scale * np.sum(rho_pts * z_pts**2 * wi)这段代码里最容易错的是scale = h/2这个系数。Gauss-Legendre积分默认区间在[-1,1]，节点 z_pts 已经换算到 [-h/2, h/2]，所以积分权重也要乘 h/2，漏掉这个系数会让刚度和惯量同时差一个数量级，频率误差会直接放大到结果不可信。
3.3 组装5×5矩阵并求解
def navier_freq(m, n): alpha = m * np.pi / a beta = n * np.pi / b K = np.zeros((5, 5)) K[0, 0] = A11 * alpha**2 + A66 * beta**2 K[0, 1] = (A12 + A66) * alpha * beta K[0, 3] = B11 * alpha**2 + B66 * beta**2 K[0, 4] = (B12 + B66) * alpha * beta K[1, 1] = A66 * alpha**2 + A11 * beta**2 K[1, 3] = (B12 + B66) * alpha * beta K[1, 4] = B66 * alpha**2 + B11 * beta**2 K[2, 2] = As * (alpha**2 + beta**2) K[2, 3] = As * alpha K[2, 4] = As * beta K[3, 3] = D11 * alpha**2 + D66 * beta**2 + As K[3, 4] = (D12 + D66) * alpha * beta K[4, 4] = D66 * alpha**2 + D11 * beta**2 + As # 对称填充 K = K + K.T - np.diag(np.diag(K)) M = np.zeros((5, 5)) M[0, 0] = I0 M[1, 1] = I0 M[2, 2] = I0 M[3, 3] = I2 M[4, 4] = I2 M[0, 3] = M[3, 0] = I1 M[1, 4] = M[4, 1] = I1 w2, vec = eigh(K, M) return np.sqrt(np.maximum(w2, 0)) / (2 * np.pi), vec freqs, modes = navier_freq(1, 1) print("前五阶频率(Hz)：", freqs[:5])注意我在对称填充时用了K = K + K.T - np.diag(np.diag(K))，因为上面只填了上三角。这种做法比手动逐个补对称项更不容易抄错。如果结果里出现 (w_2) 和 (w_3) 不相等，那基本可以断定矩阵某个交叉项写错了。
3.4 怎么确定代码算出来的数是对的
第一件事：把梯度指数设成 p=1e-8，厚度 h 改成 0.001，跑出来的基频应当非常接近经典薄板简支方板的解 ( \bar{\omega} \approx 19.7392 )。这里的无量纲定义是 (\bar{\omega} = \omega a^2 \sqrt{\rho_m h / D_m})，其中 (D_m) 是金属材料的弯曲刚度。如果这个数不对，说明积分、组装或者边界条件展开里有bug，先修好再往下走。
第二件事：把梯度指数设回p=1，算完看一眼趋势。因为陶瓷模量高于金属，p增大意味着陶瓷占比下降，板整体变软，基频应该单调下降并逐步趋近金属板的值。如果出现先升后降的诡异曲线，多半是B矩阵或I1的符号有问题，而不是物理规律变了。
第三件事：有条件的话，用一个40层均匀分层的3D实体有限元模型交叉验证。每层给不同的等效模量，层数越多越接近连续梯度。我的经验是，只要层数超过30层，分层有限元结果和FSDT的差值会小于1%，这个误差主要来自分层近似，不是你的代码。
3.5 振型可视化
拿到特征向量后，可视化是检验“解是否像样”的最直观手段。把特征向量按Navier基函数叠加，可以还原出整个面上的挠度场：
import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D N_x, N_y = 41, 41 x = np.linspace(0, a, N_x) y = np.linspace(0, b, N_y) X, Y = np.meshgrid(x, y) # 以第一阶模态的 w 分量振型为例 mode = modes[:, 0] W_amp = mode[2] # w 的幅值 Z = W_amp * np.sin(np.pi * X / a) * np.sin(np.pi * Y / b) fig = plt.figure(figsize=(8, 5)) ax = fig.add_subplot(111, projection='3d') ax.plot_surface(X, Y, Z, cmap='viridis') ax.set_xlabel('x'); ax.set_ylabel('y'); ax.set_zlabel('w') plt.show()如果画出来的第一阶模态是半个正弦波、中间最大四周为零，那求解器基本可以进入参数扫描环节了。如果出现锯齿形或者边界不归零，回查边界条件展开项是否满足简支条件。
4. 不满足于简支：用DQ法把手伸向任意边界条件
4.1 为什么Navier解走到这里就到头了
Navier解能处理的边界条件非常有限，它要求对边的边界条件配对满足特定三角展开，最常见的就是四边简支。工程上固支方板、悬臂板、四边自由板才是常态，这些情况用解析法几乎无解。传统的处理方法是用有限元，但有限元网格加密、收敛验证、后处理一套流程下来，节奏很慢。这时候可以上微分求积法（DQ法）：用少量全局节点把偏微分方程直接离散成代数方程，对光滑问题收敛速度飞快。
4.2 DQ加权系数的生成方法
DQ法的基本思想很简单：函数在某个节点上的导数，用所有节点函数值的加权和来近似。关键是加权系数怎么算。经典做法是在Chebyshev-Gauss-Lobatto（CGL）节点上构造一阶加权系数矩阵，然后用矩阵乘法得到二阶系数矩阵：
def dq_weights(x): N = x.size - 1 W1 = np.zeros((N + 1, N + 1)) for i in range(N + 1): for j in range(N + 1): if i == j: continue num = 1.0 den = 1.0 for k in range(N + 1): if k == i or k == j: continue num *= (x[i] - x[k]) den *= (x[j] - x[k]) W1[i, j] = num / ((x[i] - x[j]) * den) for i in range(N + 1): W1[i, i] = -np.sum(W1[i, :]) W2 = W1 @ W1 return W1, W2 N = 15 x_node = 0.5 - 0.5 * np.cos(np.pi * np.arange(N + 1) / N) W1, W2 = dq_weights(x_node)四边简支方板的双调和算子在二维张量积网格上可以写成：
I = np.eye(N + 1) L = np.kron(W2, I) + 2 * np.kron(W1, W1) + np.kron(I, W2)这个 L 是离散的双调和算子，对应薄板振动方程 (D\nabla^4 w = \rho h \omega^2 w)。拿到L之后，把边界条件通过行变换施加进去，就变成一个标准的广义特征值问题：(L_{mod} \mathbf{w} = \lambda \mathbf{w})。
4.3 边界条件的三种处理手法
DQ法施加边界条件比有限元更讲究，常见套路有三种：
- 直接替换法：把边界点对应行替换成单位向量（w=0）或一阶差分近似（w_x=0），实现最简单，容易破坏刚度矩阵的对称性，但小节点数下仍可用。
- 消去法：把边界节点自由度从未知量中消去，再求解内部节点。精度高，但代码要处理索引映射，稍繁琐。
- 方程约束法：把边界条件当作约束方程，用拉格朗日乘子或者罚函数引入。通用性强，适合复杂边界组合。
我做薄板DQ时更倾向于消去法，它的稳定性最好。以四边固支为例，固支条件是 (w=0) 且 (w_{,x}=0)（或 (w_{,y}=0)）。边界节点的位移直接置零，边界导数的约束则通过把边界点的相邻内部点值代入一阶加权系数方程，把边界点的二阶导数约束化掉。这一步写起来要小心索引，但逻辑非常机械。
4.4 从薄板DQ到FSDT板DQ的扩展思路
你可能会问：上面这套只是薄板CPT的DQ，和FGM板的FSDT有什么关系？关系在于组装思路完全相同。FSDT的DQ离散并不是构建一个大双调和算子，而是把每个高斯点的5个自由度 ((u,v,w,\phi_x,\phi_y)) 全部拉直成一个长向量，每个控制方程在内部点上写成关于这个长向量的代数方程，边界条件则对相应的边界自由度做约束。最后得到的仍然是一个广义特征值问题。
我自己实测下来，采用 CGL 节点时，(N=13) 到 (N=17) 个节点就能让前五阶频率收敛到小数点后三位。这是DQ法最大的好处：节点少、不用划分网格、改边界条件只改几行约束代码。如果你需要处理非矩形域或者变厚度板，建议直接转向有限元或等几何分析，DQ在规则域上的优势会更明显。
5. 参数扫描里藏着工程答案：梯度指数、厚跨比与频率的定量关系
5.1 梯度指数 p 对基频的单调性
把上一章的求解器包进一个循环，扫描梯度指数 (p)。这是最简单的参数扫描代码：
ps = [0.0, 0.2, 0.5, 1.0, 2.0, 5.0, 10.0] for p_val in ps: p = p_val # 重新计算 E_z, rho_z, 再进行厚度积分 # （需要把第3章的积分代码放进函数里，或者用闭包） freq, _ = navier_freq(1, 1) print(p_val, freq[0])结果会呈现清晰的单调下降趋势：p=0时板是纯陶瓷，刚度最大，频率最高；p增大意味着金属比例增加，板变软，频率下降；当p超过5，结果会非常接近纯金属板，后续再增大p对频率的影响已经很小。这个趋势很有用，工程上想减重又不想牺牲太多刚度，p取1到2附近往往是最佳区间。
5.2 厚跨比 a/h：什么时候不能再用薄板理论
把 a/h 从100扫到5，FSDT和CPT的差异是一条上升的曲线。在 a/h=100 时两者几乎重合，差异不到0.1%；a/h=20 时差异约1%到2%，尚在工程容差内；到 a/h=10 时差异开始明显，能达到5%以上；等到 a/h=5，差异可能超过10%。这意味着，如果你用薄板理论去算一块中厚FGM板的频率，结果会偏大，也就是偏不安全。
所以我的建议是：a/h 大于20时，可以用CPT快速估算；低于20，老老实实上FSDT；低于5，最好再往TSDT或三维实体单元走。这个判断标准放在FGM板上尤其重要，因为梯度材料让剪切刚度进一步降低，实际修正量比均质板更大。
5.3 高阶模态的剪切效应与多阶校核
很多人做模态分析只看第一阶频率，但FGM板的剪切效应对高阶模态的影响更严重。从FSDT的刚度矩阵可以看出，剪切项（带As的项）在模态阶次升高时，由于 (\alpha^2+\beta^2) 增大，对频率的修正占比也在变大。如果你的结构在工作频率范围内可能激发出第三、第四阶模态，只按基频校核会明显低估风险。
我在实际项目里一般会扫描前五阶，然后看每阶对应的振型是弯曲主导还是扭转主导。对于矩形板，正方形板会出现 (w_2=w_3) 的重频现象，这是对称性导致的，不是错误。如果参数扫描时发现相邻模态频率曲线有交叉，也要留意是否有振型交换，这对后续的响应分析和优化迭代非常关键。
6. 这五个坑，我替你们踩了一遍
6.1 剪切修正系数5/6是均质板的答案，不是FGM的
FSDT最大的软肋就是剪切修正系数。经典值5/6是从均质各向同性板推导来的，但FGM板厚度方向剪切模量在变化，严格的修正系数并不是一个常数。如果直接沿用5/6，中厚FGM板的频率会偏低一些。学术文献里有人专门推导过FGM板的剪切修正系数，结果通常落在2/3到5/6之间，具体值依赖p和材料模量比。
工程上的务实做法是：承认这个不确定性，在结果对比中注明使用的是5/6修正；如果你的项目对精度要求高，建议升级到TSDT，或者用三维实体模型交叉验证。用力学的话说，这不是代码bug，是模型误差。
6.2 厚度方向的积分别懒，别用低阶梯形法
积分方法决定刚度系数的精度。我见过有人为了省事用10层梯形法做厚度积分，结果p大于2时频率误差能到百分之几。原因在于FGM的指数变化在厚度方向上分布不均，线性近似跟不上。用Gauss-Legendre积分之后，哪怕用30个高斯点，积分误差就能压到机器精度以下。代码里设N_GP=50完全不浪费，积分才花几微秒，特征值求解才是主要开销。
6.3 特征值求解出现“虚频”先别慌
我在写第一版求解器时，遇到过 eigh 解出来的特征值有负值，取根号后一片NaN。排查下来是B矩阵没有对称填充导致K不对称，负特征值就冒出来了。所以如果你发现np.sqrt(w2)里出现NaN，第一件事不是检查物理模型，而是检查K矩阵是否对称。另外，用scipy.linalg.eigh而不是eig，因为前者能利用对称性，数值稳定性更好，速度也更快。对5×5的小矩阵可能看不出差别，但DQ法组装出来的大矩阵，差别会非常明显。
6.4 无量纲频率的“口径”问题
文献里的无量纲频率至少有三种常见定义：(\omega a^2\sqrt{\rho_c h/D_c})、(\omega h\sqrt{\rho_c/E_c})、(\omega a^2/h \sqrt{\rho_c/E_c})。这三种定义数量级差别很大，直接对比很容易把结果搞错。我的做法是在所有代码里先输出绝对频率（Hz），最后需要对比文献时再显式换算，换算式写清楚，而不是在代码里隐式转换。这样虽然多写两行，但避免了“对不上文献时根本不知道是自己算错还是无量纲定义不同”的尴尬。
6.5 单位制统一与Python环境的小事
弹模用Pa、密度用kg/m³、长度用m、频率输出Hz，这套SI单位制必须自始至终统一。我踩过最无语的坑是把密度写成了g/cm³，频率差了30多倍。环境方面，真不用追求最新Python版本，3.10或者3.11都很稳，装好numpy、scipy、matplotlib三个包就够。如果你用的是Anaconda，直接conda install numpy scipy matplotlib就行，别去折腾复杂的环境配置，把力气留在调矩阵上。
数值分析这个活，很多时候最花时间的不是推导，而是复盘“结果为什么不对”。功能梯度板问题不算新，但用Python把它从材料模型一路拆到频率输出，整个过程既能帮你复习板壳理论，又给你留了一整套可以随意扩展的代码框架：改梯度分布函数就能模拟S-FGM，加一个线性阻尼项就能做复模态，把边界条件换成固支就能直接进入工程实践。这套代码我现在还在用，每次加新功能都比重新建模省事太多。
