简介这份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 GPaE2 8-12 GPaG12 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/XcF2 1/Yt - 1/YcF66 1/S²。F12 是交互项一般取 -1/2 * sqrt(F11*F22) 或取零值保守估计。Tsai-Wu 准则的输出是失效指数 FIFI 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 GPaE2 9.5 GPaG12 5.2 GPaν12 0.31。强度参数为 Xt 1500 MPaXc 1200 MPaYt 45 MPaYc 180 MPaS 75 MPa。受载状态为 Nx 100 N/mmNy 0Nxy 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} | fHashin{ {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-WuF120 Hashin 基体拉伸 Tsai-WuF12 取经验值 最大应力准则。这个排序的物理背景是 F12 的取值改变了失效包络面的曲率方向F120 时包络面内凹判据更严。选准则时不应只看保守性也要看失效模式的物理含义。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, repeatn): # 强制对称铺层只枚举一半角度 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 * ε 的基本关系。这种最小验证能在一分钟内确认整套代码链路没有符号错误是每次修改代码后应该做的回归测试。我在调整失效准则参数或修改铺层逻辑后都会先跑这段验证再进入正式分析。本文还有配套的精品资源点击获取