简介：这份Python代码包面向复合材料力学方向的学生与工程师，围绕经典层压理论（CLT）实现复合层定义、层压板铺层、应力应变分布计算与失效准则判定，覆盖杨氏模量、刚度矩阵、强度校核等核心环节，适合课程设计、科研验证或工程预分析。压缩包共63个文件，以26个Python脚本为核心计算模块，辅以21个pyc编译文件与10张PNG结果图，另有2份PDF参考手册（含Nastran手册）、LaTeX排版文件、辅助脚本及README说明文档，整体仅2.9MB，下载后即可快速上手。该资源已有307人学习浏览，在复合材料计算场景中具备一定参考价值。通过运行脚本可查看每一层应力应变分布图、失效步骤与层合板强度校核结果，代码内还包含Chamis细观力学模型手工计算示例，便于对照理解CLT完整求解流程。
1. 复合层压材料性能分析：从弹性常数到失效判据的完整计算链路
做复合材料结构设计的人，迟早会遇到一个尴尬：手算层压板的等效刚度可以，但一旦涉及逐层应力、失效指数和强度比，就不得不依赖商业软件，而很多商业软件的黑盒交互方式让你根本不知道结果是怎么来的。复合层压材料性能的完整计算，本质上是三条线的交汇：单层板的宏观弹性常数（杨氏模量、剪切模量、泊松比）、经典层合板理论（CLT）下的刚度矩阵组装、以及失效准则对逐层应力状态的判读。这三条线串起来，就是一段不超过三百行的 Python 代码能覆盖的完整链路。这篇文章要解决的就是：如何用 Python 从零搭起这套分析工具，覆盖从材料输入、偏轴刚度转换、ABD 矩阵求解到 Tsai-Wu 与 Hashin 准则判别的全过程，并给出能直接复制运行的代码和参数选取依据。
2. 复合层压材料性能的理论底座：杨氏模量与刚度矩阵的数学结构
2.1 单层板的正轴刚度：为什么杨氏模量不是唯一输入
要分析复合层压材料性能，第一步不是写代码，而是认清单层板的力学描述方式。单向复合材料板在材料主方向上有五个独立弹性常数：纵向杨氏模量 E1、横向杨氏模量 E2、面内剪切模量 G12，以及主泊松比 ν12 和次泊松比 ν21。很多人误以为有了 E1、E2 就够，实际上在二维应力状态下，正轴柔度矩阵 S 的完整形态是：
S = [[1/E1, -ν12/E1, 0], [-ν21/E2, 1/E2, 0], [0, 0, 1/G12]]其中 ν21 不是独立参数，由ν21 = ν12 * E2 / E1 约束。这是复合材料与各向同性材料最大的区别：刚度矩阵非对角项由两个不同模量和两个泊松比联合决定，任何一个参数的测量误差都会在后续应力计算中被放大。所以做性能分析前，必须校验输入数据的物理合理性，常见手段是检查正定条件——刚度矩阵必须正定，即 ν12 * ν21 < 1。
实际工程中，E1 由纤维主导，E2 和 G12 由基体主导。碳纤维/环氧体系的典型值范围是 E1 = 120-180 GPa，E2 = 8-12 GPa，G12 = 4-7 GPa，ν12 = 0.25-0.35。如果输入数据落在这个范围外，先怀疑数据来源而不是程序。我做材料参数校核时，习惯先把五个常数代入正定检查，再进入后续计算，这能挡掉一大部分脏数据。
2.2 偏轴刚度转换：层合板每一层的刚度方向不同
复合层压材料性能计算的核心难点不在正轴，而在偏轴。每个单层按特定铺层角铺放，其材料主方向与层合板参考坐标系之间存在夹角 θ。此时需要把正轴刚度矩阵 Q 旋转到偏轴方向。刚度矩阵的转换公式为：
Q_bar = T_inv * Q * T_inv^T （张量转换形式）展开写就是工程上常用的 Q11_bar、Q22_bar、Q12_bar、Q66_bar 以及耦合项 Q16_bar、Q26_bar。其中 Q16_bar 和 Q26_bar 的引入意味着偏轴层同时存在拉剪耦合，这是复合材料层合板区别于各向同性板的重要特征。对于对称均衡铺层，这两个耦合项会在层合板层面被抵消，但对非对称铺层，它们会进入 ABD 矩阵的 B 子矩阵，导致弯曲-拉伸耦合。
Python 实现时，我一般用角度转弧度后直接按三角函数表达式计算转换矩阵，而不是先构造四阶张量再做缩并。虽然张量方法通用性更好，但工程分析中只需要二维情况，显式表达式更直观也更容易调试。下面给出偏轴刚度转换的具体代码。
import numpy as np def rotate_stiffness(E1, E2, G12, nu12, theta_deg): """计算偏轴刚度矩阵 Q_bar (3x3)。 参数对应单层板正轴弹性常数，theta_deg 为铺层角（度）。 """ nu21 = nu12 * E2 / E1 # 正定检查 if nu12 * nu21 >= 1.0: raise ValueError("泊松比乘积 >= 1，刚度矩阵不正定，请检查输入") # 正轴刚度矩阵 Q denom = 1.0 - nu12 * nu21 Q11 = E1 / denom Q22 = E2 / denom Q12 = nu12 * E2 / denom Q66 = G12 theta = np.radians(theta_deg) c, s = np.cos(theta), np.sin(theta) c2, s2 = c*c, s*s c4, s4 = c2*c2, s2*s2 # 偏轴刚度显式计算（经典 CLT 表达） Q11_bar = Q11*c4 + 2*(Q12 + 2*Q66)*s2*c2 + Q22*s4 Q22_bar = Q11*s4 + 2*(Q12 + 2*Q66)*s2*c2 + Q22*c4 Q12_bar = (Q11 + Q22 - 4*Q66)*s2*c2 + Q12*(c4 + s4) Q66_bar = (Q11 + Q22 - 2*Q12 - 2*Q66)*s2*c2 + Q66*(c4 + s4) Q16_bar = (Q11 - Q12 - 2*Q66)*s*c2*c - (Q22 - Q12 - 2*Q66)*s2*s*c Q26_bar = (Q11 - Q12 - 2*Q66)*s2*s*c - (Q22 - Q12 - 2*Q66)*s*c2*c return np.array([ [Q11_bar, Q12_bar, Q16_bar], [Q12_bar, Q22_bar, Q26_bar], [Q16_bar, Q26_bar, Q66_bar] ])这段代码的关键点在于 Q16_bar 和 Q26_bar 的计算，它们的符号直接取决于铺层角的正负，而这两个耦合项在后续 ABD 矩阵组装中是判别层合板是否发生拉剪耦合的依据。代码里我显式做了正定检查，因为不少工程数据表格里的泊松比是近似值，直接把近似值丢进公式可能得到负刚度，这在物理上不可能。
2.3 层合板 ABD 矩阵组装：刚度矩阵从单层到整体的升维
有了每一层的偏轴刚度 Q_bar，下一步是沿着厚度方向积分，组装层合板的面内刚度矩阵 A、耦合刚度矩阵 B 和弯曲刚度矩阵 D。这是复合层压材料性能分析中最关键的一步，也是容易出错的地方。
三个子矩阵的定义是：
A_ij = Σ (Q_bar_ij)_k * (z_k - z_{k-1}) B_ij = 1/2 * Σ (Q_bar_ij)_k * (z_k^2 - z_{k-1}^2) D_ij = 1/3 * Σ (Q_bar_ij)_k * (z_k^3 - z_{k-1}^3)其中 z_k 是第 k 层底面到层合板中面的距离。注意所有 z 坐标都是相对于中面定义的，即使层合板不对称，中面也是计算基准。代码实现时，先按铺层顺序逐层累加厚度确定 z 坐标，再循环组装。
def assemble_abd(layers): """输入层合板铺层信息，输出 ABD 矩阵。 layers: list of dict，每个 dict 含 E1, E2, G12, nu12, thickness(mm), theta(deg) """ total_thickness = sum(l['thickness'] for l in layers) z0 = -total_thickness / 2.0 A = np.zeros((3, 3)) B = np.zeros((3, 3)) D = np.zeros((3, 3)) z_prev = z0 for l in layers: Q_bar = rotate_stiffness(l['E1'], l['E2'], l['G12'], l['nu12'], l['theta']) z_curr = z_prev + l['thickness'] # 厚度坐标积分 A += Q_bar * (z_curr - z_prev) B += 0.5 * Q_bar * (z_curr**2 - z_prev**2) D += (1/3.0) * Q_bar * (z_curr**3 - z_prev**3) z_prev = z_curr return A, B, D组装完成后，6x6 的 ABD 矩阵把层合板的面内力和弯矩与中面应变和曲率关联起来：
[N] = [A B] [epsilon_0] [M] [B D] [kappa_0]对于对称铺层，B 矩阵为零，面内问题和弯曲问题解耦，这也是大多数实际结构采用对称铺层的原因。判断对称性有个快速方法：铺层角序列关于中面对称，且对应层材料相同。如果你发现 B 矩阵非零导致面内拉伸伴随弯曲变形，先检查铺层序列是否按对称排列。
3. 用 Python 实现强度和失效准则应用：Tsai-Wu 与 Hashin 的可执行方案
3.1 逐层应力恢复：失效判据的前提条件
失效准则应用的前提是得到每一层的应力状态。ABD 矩阵给出的是中面应变和曲率，每一层的真实应力需要先由中面应变加弯曲应变合成该层的总应变，再乘以该层的偏轴刚度矩阵还原为应力。
具体流程是：已知外力 N 和 M，先求中面应变 epsilon_0 和曲率 kappa_0，然后第 k 层的应变为 epsilon_k = epsilon_0 + z_k * kappa_0，最终应力为 sigma_k = Q_bar_k * epsilon_k。代码实现如下。
def layer_stresses(A, B, D, layers, N, M): """计算每一层上下表面的应力，单位为 MPa。 N: 面内力向量 [Nx, Ny, Nxy] N/mm M: 弯矩向量 [Mx, My, Mxy] N*mm/mm """ abd = np.block([[A, B], [B, D]]) strain_curv = np.linalg.solve(abd, np.concatenate([N, M])) eps0 = strain_curv[:3] kappa = strain_curv[3:] stresses = [] z_prev = -sum(l['thickness'] for l in layers) / 2.0 for l in layers: Q_bar = rotate_stiffness(l['E1'], l['E2'], l['G12'], l['nu12'], l['theta']) z_curr = z_prev + l['thickness'] # 上下表面的应变与应力 eps_top = eps0 + z_curr * kappa eps_bottom = eps0 + z_prev * kappa sigma_top = Q_bar @ eps_top sigma_bottom = Q_bar @ eps_bottom stresses.append({layer['theta']: (sigma_bottom, sigma_top, z_prev, z_curr)}) z_prev = z_curr return stresses注意这里的应力是偏轴坐标系下的值，也就是沿层合板参考方向 x、y 和 xy 方向的应力分量。失效准则应用时，需要把这些应力转换回每一层的材料主方向，即正轴坐标系下的 σ1、σ2、τ12。转换方法是把偏轴应力逆旋转回正轴，这本质上也是坐标变换，只不过应用对象是应力张量而不是刚度矩阵。
def rotate_stress_to_principal(sigma_xy, theta_deg): """把偏轴应力转换到材料主方向。 sigma_xy: [sigma_x, sigma_y, tau_xy] 单位 MPa """ th = np.radians(theta_deg) c, s = np.cos(th), np.sin(th) T = np.array([ [c*c, s*s, 2*s*c], [s*s, c*c, -2*s*c], [-s*c, s*c, c*c - s*s] ]) return T @ sigma_xy这里有一个常见的坑：很多人直接把偏轴应力代入失效准则，得到的结果完全错误。纤维方向的拉伸强度通常比横向高一个数量级，如果应力方向不转换，横向应力被误判为纵向应力，失效指数会严重失真。所以失效准则应用的第一步永远是坐标变换，而不是选公式。
3.2 Tsai-Wu 张量准则：参数标定与失效指数计算
Tsai-Wu 准则在多轴应力状态下有良好的适用性。其失效判据为张量多项式：
F_ij * sigma_i * sigma_j + F_i * sigma_i = 1对于二维正轴应力状态展开后，表达式为：
F11*σ1² + F22*σ2² + F66*τ12² + 2*F12*σ1*σ2 + F1*σ1 + F2*σ2 = 1其中系数由强度参数标定：F11 = 1/(XtXc)，F22 = 1/(YtYc)，F1 = 1/Xt - 1/Xc，F2 = 1/Yt - 1/Yc，F66 = 1/S²。F12 是交互项，一般取 -1/2 * sqrt(F11*F22) 或取零值保守估计。Tsai-Wu 准则的输出是失效指数 FI，FI < 1 安全，FI = 1 临界，FI > 1 失效。由于表达式是二次型，FI 不线性正比于载荷，所以工程上通常用反推法求失效载荷——把载荷乘放大系数 iter 直到 FI 逼近 1。
def tsai_wu_fi(sigma1, sigma2, tau12, Xt, Xc, Yt, Yc, S): """Tsai-Wu 失效指数计算。所有强度参数单位 MPa，压强度取正数。 """ F11 = 1.0 / (Xt * Xc) F22 = 1.0 / (Yt * Yc) F1 = 1.0 / Xt - 1.0 / Xc F2 = 1.0 / Yt - 1.0 / Yc F66 = 1.0 / (S * S) # F12 交互项，采用各向同性近似折减 F12 = -0.5 * np.sqrt(F11 * F22) fi = (F11*sigma1**2 + F22*sigma2**2 + F66*tau12**2 + 2*F12*sigma1*sigma2 + F1*sigma1 + F2*sigma2) return fi参数标定是这个准则应用中最容易出错的地方。压强度 Xc、Yc 在公式里必须取正数，但有些文献里压强度写成负值，直接代入会得到错误的线性项系数。另外，如果只有拉伸强度而没有压缩强度数据，F1 和 F2 的标定是不完整的，此时宁可不考虑线性项也不能乱猜。交互项 F12 的取值对失效包络面形状影响很大，在双轴拉压工况下结果差异明显，我一般建议先取零，再用双轴实验数据修正。
3.3 Hashin 准则：区分纤维失效与基体失效的判别逻辑
Hashin 准则的价值在于区分失效模式，而不是只给一个失效指数。它把失效分成纤维拉伸、纤维压缩、基体拉伸、基体压缩四种模式，每种模式有独立的判据。二维状态下常用的是纤维拉伸与压缩、基体拉伸与压缩四个公式。
def hashin_failure(sigma1, sigma2, tau12, Xt, Xc, Yt, Yc, S): """Hashin 准则逐模式判别。 返回 dict，包含各失效模式的指数和触发布尔值。 """ result = {} # 纤维拉伸（sigma1 >= 0） if sigma1 >= 0: ff = (sigma1 / Xt)**2 + (tau12 / S)**2 mode = 'fiber_tension' else: ff = (sigma1 / Xc)**2 mode = 'fiber_compression' result[mode] = ff # 基体拉伸（sigma2 >= 0） if sigma2 >= 0: mt = (sigma2 / Yt)**2 + (tau12 / S)**2 mode = 'matrix_tension' else: # 基体压缩考虑剪应力贡献 mt = ((sigma2 / (2*S))**2 + ((Yc / (2*S))**2 - 1) * (sigma2 / Yc) + (tau12 / S)**2) mode = 'matrix_compression' result[mode] = mt result['failed'] = {k: v >= 1.0 for k, v in result.items()} return resultHashin 准则应用中的关键判断是模式优先级：如果纤维拉伸指数和基体拉伸指数同时超限，结构设计上应从纤维失效开始处理，因为纤维断裂通常是灾难性的，而基体开裂可能只意味着继续承载能力下降但结构未解体。工程上，基体拉伸失效指数在 1.0-1.5 之间时，很多设计规范允许带伤工作，但纤维失效指数超过 1.0 就认为结构达到承载极限。这个差异源于两种失效模式的后果严重程度不同。
4. 完整案例实操：从铺层设计到失效准则应用的一站式 Python 代码
4.1 问题定义与铺层参数输入
以一个典型碳纤维层合板为例：铺层为 [0/90/±45]s，共 8 层，单层厚度 0.125 mm，总厚度 1.0 mm。材料为 T300/环氧体系，弹性常数为 E1 = 135 GPa，E2 = 9.5 GPa，G12 = 5.2 GPa，ν12 = 0.31。强度参数为 Xt = 1500 MPa，Xc = 1200 MPa，Yt = 45 MPa，Yc = 180 MPa，S = 75 MPa。受载状态为 Nx = 100 N/mm，Ny = 0，Nxy = 0，即单向拉伸。
这是一个标准的面内载荷问题。由于铺层对称且均衡，B 矩阵为零，但仍按完整流程计算以验证代码正确性。把所有参数组织成数据结构，方便后续修改铺层角度或载荷。
material = { 'E1': 135e3, 'E2': 9.5e3, 'G12': 5.2e3, 'nu12': 0.31, 'Xt': 1500.0, 'Xc': 1200.0, 'Yt': 45.0, 'Yc': 180.0, 'S': 75.0 } layer_thickness = 0.125 layers = [] for theta in [0, 90, 45, -45, -45, 45, 90, 0]: layers.append({ 'E1': material['E1'], 'E2': material['E2'], 'G12': material['G12'], 'nu12': material['nu12'], 'thickness': layer_thickness, 'theta': theta }) N = np.array([100.0, 0.0, 0.0]) # N/mm M = np.array([0.0, 0.0, 0.0]) # N*mm/mm参数选取说明：E1 和 E2 的量级差异是复合材料各向异性的直接体现，数值上差了 14 倍，这意味着在 90 度铺层中，载荷主要由基体承担，应力水平会显著高于 0 度层。强度参数中 Yt 只有 45 MPa，是整套参数中的最薄弱环节，失效大概率会从基体拉伸模式触发。
4.2 刚度矩阵计算与结果验证
调用前两章的功能函数，完成 ABD 矩阵组装和逐层应力计算。为了验证正确性，可以做一个中间检查：在相同载荷下，[0/90/±45]s 层合板的等效面内刚度应该介于所有层全是 0 度和全 90 度之间。更直观的验证手段是计算等效工程常数：Ex = (A11*A22 - A12²) / (A22 * t_total)。
A, B, D = assemble_abd(layers) # 等效面内工程常数验证 t_total = sum(l['thickness'] for l in layers) Ex = (A[0,0]*A[1,1] - A[0,1]**2) / (A[1,1] * t_total) Ey = (A[0,0]*A[1,1] - A[0,1]**2) / (A[0,0] * t_total) Gxy = A[2,2] / t_total nu_xy = A[0,1] / A[1,1] print(f"等效纵向模量 Ex = {Ex/1e3:.1f} GPa") print(f"等效横向模量 Ey = {Ey/1e3:.1f} GPa") print(f"等效剪切模量 Gxy = {Gxy/1e3:.1f} GPa") stresses = layer_stresses(A, B, D, layers, N, M) for layer in stresses: for theta, (sig_b, sig_t, _, _) in layer.items(): sig_principal_b = rotate_stress_to_principal(sig_b, theta) sig_principal_t = rotate_stress_to_principal(sig_t, theta) print(f"铺层角 {theta:>3}°: 下表面正轴应力 σ1={sig_principal_b[0]:.2f}, " f"σ2={sig_principal_b[1]:.2f}, τ12={sig_principal_b[2]:.2f} MPa")输出结果中值得关注的是 90 度层的横向应力 σ2。由于 90 度层纤维方向垂直于载荷方向，载荷大部分由基体传递，σ2 会明显高于其他铺层。如果 σ2 已经接近 Yt = 45 MPa，说明层合板在该载荷下已经接近基体拉伸失效。这是复合层压材料性能分析的典型结论：层合板的强度下限由最薄弱的铺层方向和失效模式决定，而不是简单地平均分配。
4.3 失效指数计算与失效模式解读
对每一层的正轴应力分别计算 Tsai-Wu 失效指数和 Hashin 四种模式指数。由于是纯面内拉伸，上下表面应力符号相同，只取上表面结果即可。
print("\n失效指数评估（上表面）:") for layer in stresses: for theta, (_, sig_t, _, _) in layer.items(): sp = rotate_stress_to_principal(sig_t, theta) fi_tsai = tsai_wu_fi(sp[0], sp[1], sp[2], material['Xt'], material['Xc'], material['Yt'], material['Yc'], material['S']) h = hashin_failure(sp[0], sp[1], sp[2], material['Xt'], material['Xc'], material['Yt'], material['Yc'], material['S']) failed_modes = [k for k, v in h['failed'].items() if v] print(f"θ={theta:>3}° | Tsai-Wu FI={fi_tsai:.3f} | " f"Hashin={ {k: round(v,3) for k,v in h.items() if k!='failed'} }") if failed_modes: print(f" → 触发失效模式: {failed_modes}")分析逻辑上，先看 Tsai-Wu 的全局失效指数，再看 Hashin 的逐模式指数定位失效源头。如果 Tsai-Wu FI 小于 1 但某层 Hashin 基体拉伸指数大于 1，说明该层正在承担超过其横向强度的应力，层合板虽未整体失效，但已经出现局部基体开裂。这种跨准则的交叉验证比单一准则可靠得多，也是工程审查时最有说服力的分析路径。
4.4 强度比与安全裕度计算：失效准则应用的反向思考
失效准则应用不仅用于验算某一载荷下是否失效，更常用的是反推强度比——即层合板能承载的最大载荷倍数。这类似于有限元分析中的安全裕度计算。定义强度比 R = 极限载荷 / 当前载荷，则各层应力按 R 比例缩放。对于 Tsai-Wu 准则，由于表达式中同时包含应力的二次项和线性项，求 R 需要解一元二次方程：
A*R² + B*R = 1 其中 A = F11*σ1² + F22*σ2² + F66*τ12² + 2*F12*σ1*σ2, B = F1*σ1 + F2*σ2 R = (-B + sqrt(B² + 4A)) / (2A)def tsai_wu_strength_ratio(sig1, sig2, tau12, Xt, Xc, Yt, Yc, S): """计算 Tsai-Wu 强度比 R，表示极限载荷与当前载荷的倍数关系。""" F11 = 1.0 / (Xt * Xc) F22 = 1.0 / (Yt * Yc) F1 = 1.0 / Xt - 1.0 / Xc F2 = 1.0 / Yt - 1.0 / Yc F66 = 1.0 / (S * S) F12 = -0.5 * np.sqrt(F11 * F22) A = (F11*sig1**2 + F22*sig2**2 + F66*tau12**2 + 2*F12*sig1*sig2) B = F1*sig1 + F2*sig2 R = (-B + np.sqrt(B**2 + 4*A)) / (2*A) return R # 在所有层中找最小强度比，即最危险层 min_R = float('inf') weakest_layer = None for layer in stresses: for theta, (_, sig_t, _, _) in layer.items(): sp = rotate_stress_to_principal(sig_t, theta) R = tsai_wu_strength_ratio(sp[0], sp[1], sp[2], material['Xt'], material['Xc'], material['Yt'], material['Yc'], material['S']) if R < min_R: min_R = R weakest_layer = theta print(f"最小强度比 R_min = {min_R:.3f}，出现在 {weakest_layer}° 层")强度比的意义在于给出了直观的载荷裕度。R = 2.0 意味着当前载荷还可以翻倍才会达到 Tsai-Wu 预测的失效点。在工程设计中，飞机结构一般要求 R ≥ 1.5，汽车部件要求 R ≥ 1.2 左右。如果最薄弱层的 R 不满足要求，调整铺层角度比增加层数更高效——因为失效往往由横向应力主导，改变载荷传力路径比堆厚度更经济。
5. 参数灵敏度与常见误区：复合层压材料性能分析中的隐藏陷阱
5.1 泊松比与剪切模量对失效判据的影响程度
复合层压材料性能分析中，多数工程师会把精力放在 E1 和 E2 的精度上，但实际对结果影响更大的是 G12 和 ν12。原因在于：偏轴刚度矩阵中的耦合项 Q16_bar 和 Q26_bar 包含了 G12 的贡献，而层间剪应力的大小直接取决于这两个耦合项。当铺层角为 45 度时，G12 每偏差 10%，层间剪应力会偏差接近 7%-8%，这个误差传递率在失效准则应用中是显著不可忽略的。
ν12 的影响则更隐蔽。ν12 同时出现在柔度矩阵的非对角项和正定检查中，如果 ν12 偏大，可能导致 A 矩阵的等效泊松比 νxy 偏大，进而影响弯曲分析的曲率分布。建议做分析前先对 G12 做 ±5% 的扰动测试，观察失效指数变化幅度。如果变化超过 10%，说明当前分析结果对剪切模量敏感，需要在材料测试环节重点控制 G12 的测量精度。
5.2 失效准则选用策略与保守性排序
不同失效准则应用在同一工况下会给出不同结果。经验上的保守程度排序大致为：Tsai-Wu（F12=0）> Hashin 基体拉伸 > Tsai-Wu（F12 取经验值）> 最大应力准则。这个排序的物理背景是 F12 的取值改变了失效包络面的曲率方向，F12=0 时包络面内凹，判据更严。
选准则时不应只看保守性，也要看失效模式的物理含义。Hashin 准则能区分纤维和基体失效，适合做损伤容限评估；Tsai-Wu 适合做全局强度校核和初步设计筛选。如果目标是写论文或做失效机理研究，应选 Hashin 或 Puck 等具有物理失效模式的准则；如果目标是做工程校核报告，Tsai-Wu 加 1.5 安全系数是行业普遍接受的做法。
5.3 层合板理论适用边界：什么时候结果不可信
经典层合板理论的假设是每个单层处于平面应力状态且层间完美粘结。这意味着三点限制：一是层合板的面内尺寸至少是厚度的 10 倍以上，否则自由边效应导致层间应力不可忽略；二是载荷不能引起局部屈曲，否则 ABD 矩阵的线性刚度关系失效；三是每一层内应力沿厚度均匀，对于厚度方向剪切变形显著的厚板，CLT 的精度会下降。
当分析对象是层合板自由边附近区域，或结构存在开孔、缺口时，纯 CLT 的结果只具有参考意义。此时应把 Python 代码算出的远场应力作为边界条件输入有限元模型，用三维实体单元或高阶层合单元复核局部应力。常见做法是先用这里提供的代码快速扫参数、筛铺层，再用有限元验证关键工况。这种两级分析流程能兼顾速度和精度。
6. 最后一层：铺层优化中的失效准则应用技巧与验证方法
铺层优化是复合层压材料性能分析的终极应用场景。给定载荷条件，如何在最小重量下找到最优铺层顺序，本质上是一个带约束的优化问题。优化变量是每一层的角度和厚度，约束条件包括对称性、均衡性、最小铺层比例，以及所有层的失效指数小于 1。这里给出一个基于暴力枚举的扫参方法，适合铺层数不超过 12 层的设计空间。
import itertools def optimize_layup(material, candidate_angles, n_layers, N, M): """枚举铺层组合，找最小总厚度下满足 Tsai-Wu 失效约束的方案。 仅演示逻辑，完整代码需按实际需求调整。 """ layer_thickness = 0.125 best = None for n in range(1, n_layers + 1): for combo in itertools.product(candidate_angles, repeat=n): # 强制对称铺层，只枚举一半角度 if n % 2 == 0: full_angles = list(combo) + list(reversed(combo)) else: full_angles = list(combo) + list(reversed(combo[:-1])) layers = [{ 'E1': material['E1'], 'E2': material['E2'], 'G12': material['G12'], 'nu12': material['nu12'], 'thickness': layer_thickness, 'theta': th } for th in full_angles] A, B, D = assemble_abd(layers) stresses = layer_stresses(A, B, D, layers, N, M) max_fi = 0.0 for layer in stresses: for theta, (_, sig_t, _, _) in layer.items(): sp = rotate_stress_to_principal(sig_t, theta) fi = tsai_wu_fi(sp[0], sp[1], sp[2], material['Xt'], material['Xc'], material['Yt'], material['Yc'], material['S']) max_fi = max(max_fi, fi) if max_fi <= 1.0: thickness = n * layer_thickness if best is None or thickness < best['thickness']: best = {'angles': full_angles, 'thickness': thickness, 'max_fi': max_fi} break # 当前层数已有可行解，不必继续同层数枚举 return best这个枚举方案效率不高，但作为验证思路足够：它展示了失效准则应用在优化中的角色——作为约束函数判断候选铺层是否可行。真正的工程优化会改用遗传算法或梯度法缩小搜索空间，但约束判断的核心逻辑不变。每评估一个铺层，就要走一遍刚度组装和失效计算的完整链路。
验证优化结果的方法有两种。第一种是自洽性验证：把优化出的铺层重新代入 Tsai-Wu 和 Hashin 两种准则，确认最小失效指数对应的铺层和失效模式是否与优化过程中的记录一致。第二种是交叉验证：用有限元软件建立同样的层合板模型，施加相同载荷，对比 Python 计算的逐层应力和失效指数。这两个结果之间允许有 5% 以内的偏差，主要来源于有限元模型的网格离散误差和 CLT 的平面应力假设。
最后给你一个具体的验证脚本思路：对同一铺层，手动计算 0 度方向拉伸刚度 A11/t，再通过施加单位应变反求应力，看是否满足 σ = Q_bar * ε 的基本关系。这种最小验证能在一分钟内确认整套代码链路没有符号错误，是每次修改代码后应该做的回归测试。我在调整失效准则参数或修改铺层逻辑后，都会先跑这段验证再进入正式分析。
本文还有配套的精品资源，点击获取
