4. 经典层合板理论刚度矩阵ABD计算

ABD 矩阵的物理含义与计算路径

经典层合板理论（CLT）把厚度方向的位移假设为线性分布，从而把三维问题压缩为参考面的面内位移 u₀, v₀ 与挠度 w₀。层合板本构关系写成：

[N]   [A  B] [ε⁰]
[M] = [B  D] [κ ]

其中 A 是面内刚度（单位 N/m），B 是拉弯耦合刚度，D 是弯曲刚度（单位 N·m）。对第 k 层，设其偏轴刚度矩阵为 Q̄ₖ，层底/层顶坐标为 zₖ₋₁, zₖ，则：

Aᵢⱼ = Σ Q̄ₖᵢⱼ (zₖ − zₖ₋₁)
Bᵢⱼ = ½ Σ Q̄ₖᵢⱼ (zₖ² − zₖ₋₁²)
Dᵢⱼ = ⅓ Σ Q̄ₖᵢⱼ (zₖ³ − zₖ₋₁³)

关键点：B 是否为零由铺层关于中面的几何对称性决定；A、D 的正定性则取决于各层 Q̄ₖ 的正定性与厚度权重是否为正。下面用一组教学示例参数完整走一遍计算。

示例对象与单层参数

取 T300/环氧 单层，厚度 t = 0.125 mm，工程常数（示例假设值）：

铺层顺序取 [0/45/−45/90]s，共 8 层，总厚 h = 1.0 mm。中面为 z=0，层界面坐标从下到上为：

z₀=−0.500, z₁=−0.375, z₂=−0.250, z₃=−0.125,
z₄= 0.000, z₅= 0.125, z₆= 0.250, z₇= 0.375, z₈= 0.500  (mm)

各层角度依次为：0°、45°、−45°、90°、90°、−45°、45°、0°。该铺层关于中面对称，因此理论上 B = 0，这正是本例要验证的第一个结论。

计算流程与关键实现

计算分四步：①由工程常数求正轴柔度 S，求逆得正轴刚度 Q；②按角度 θ 做偏轴转换得 Q̄；③按层界面坐标累加得到 A、B、D；④检查对称性与正定性。核心代码如下（Python，可直接运行）：

import numpy as np

# 1. 单层工程常数（示例值，单位统一为 GPa、mm）
E1, E2, G12, nu12 = 135.0, 9.0, 5.0, 0.30
nu21 = nu12 * E2 / E1

# 2. 正轴柔度 S 与正轴刚度 Q（平面应力）
S = np.array([[1/E1,   -nu21/E2, 0],
              [-nu12/E2, 1/E2,   0],
              [0,        0,      1/G12]])
Q = np.linalg.inv(S)

def Qbar(Q, theta_deg):
    """偏轴刚度矩阵 Q̄，theta 为铺层角（度）"""
    t = np.radians(theta_deg)
    m, n = np.cos(t), np.sin(t)
    m2, n2, mn = m*m, n*n, m*n
    T = np.array([[m2, n2,  2*mn],
                  [n2, m2, -2*mn],
                  [-mn, mn, m2-n2]])
    Tinv = np.linalg.inv(T)
    return Tinv @ Q @ T

# 3. 铺层定义：[0/45/-45/90]s，单层厚 0.125 mm
angles = [0, 45, -45, 90, 90, -45, 45, 0]
t_ply  = 0.125
h      = t_ply * len(angles)          # 1.0 mm
z = np.linspace(-h/2, h/2, len(angles)+1)   # 层界面坐标

# 4. 累加 A、B、D
A = np.zeros((3,3)); B = np.zeros((3,3)); D = np.zeros((3,3))
for k, ang in enumerate(angles):
    Qk = Qbar(Q, ang)
    zk, zk1 = z[k+1], z[k]
    A += Qk * (zk - zk1)
    B += 0.5 * Qk * (zk**2 - zk1**2)
    D += (1/3) * Qk * (zk**3 - zk1**3)

np.set_printoptions(precision=3, suppress=True)
print("A =\n", A); print("B =\n", B); print("D =\n", D)
print("A 特征值:", np.linalg.eigvalsh(A))
print("D 特征值:", np.linalg.eigvalsh(D))
print("B 最大绝对值:", np.abs(B).max())

注意单位一致性：Q 用 GPa，坐标用 mm，则 A 的单位是 GPa·mm = 10⁹ N/m，D 的单位是 GPa·mm³。若希望 A 直接以 N/m 输出，可把 Q 换算为 Pa 后再乘 mm。

结果判读与验证

运行上述代码，示例参数下可得到（数值为按给定输入直接推出的结果）：

正定性检查：A 与 D 的三个特征值应全部为正。若出现负特征值，通常意味着铺层角或厚度输入有误，或某层 Q̄ 计算时角度符号写反。对称性检查：B 的 Frobenius 范数应接近 0；若不为零，先确认铺层是否真的关于中面对称，再检查 z 坐标是否以中面为零点。

常见错误与排查

• B 不为零但铺层对称：检查 z 坐标是否以中面为零点；若从底面开始编号，需先平移。
• A 或 D 出现负特征值：检查 Q̄ 转换矩阵是否用了转置而非逆；角度正负号是否与铺层定义一致。
• 量级异常：确认 E、G 单位统一为 GPa 或 Pa，厚度统一为 mm 或 m，避免 GPa 与 Pa 混用。
• 非对称铺层误判：[0/45/−45/90] 与 [0/45/−45/90]s 的 B 矩阵完全不同，务必核对铺层符号。