简介以SIMP方法为主线的这份资料聚焦人工髋关节的轻量化与耐疲劳设计面向生物力学、机械工程领域的研究者与工程技术开发者尤其适合已在柔度最小化代码上工作、希望进一步引入疲劳约束的读者。内容以单份PDF文档呈现压缩包约1018KB、共1个PDF文件系统梳理了从CAD模型导入GMSH、载荷与边界条件定义、网格划分到结构拓扑优化的完整流程并给出材料Ti6Al4V下基于S-N曲线的疲劳寿命建模思路与参考文献。文档进一步说明动态循环载荷的建立与雨流计数法识别应力循环的方法涵盖60kg人体行走、跑步两种工况下的多级峰值循环载荷以及目标函数柔度与体积、最大应力、最大疲劳损伤等约束的数值输出和绘图需求。目前已有75人学习下载读者可据此复现疲劳约束的耦合逻辑并借助文献参考攻克建模难点。1. 从柔度最小到疲劳约束人工髋关节柄拓扑优化的下一步一个 60 kg 的成年人慢跑时髋关节接触力的峰值可以冲到体重的五到六倍也就是三千多牛顿而这个载荷会在十年里重复上千万次。很多做髋关节拓扑优化的项目第一版跑出来的结构在静力学上非常漂亮——柔度降了三成体积减了四成——但拿去做疲劳评估时柄颈根部和干骺端内侧的应力幅直接落在 Ti6Al4V 的有限寿命区S-N 曲线上一查寿命只有几十万次循环。问题不在优化器而在优化列式里压根没有疲劳这一项。手上这套代码已经把 SIMP 的柔度最小化、体积约束、GMSH 网格、STL 后处理全部封装成了流水线输入 .step、输出 .stl 全自动。现在要做的是在这条流水线里插进两件东西一个是随时间变化的动态载荷谱一个是能把应力历史换算成损伤值并反馈给优化器的疲劳约束。适合谁看做过 SIMP 但没碰过疲劳约束的攻城狮以及正在用 Python 搭生物力学优化流程的研究人员。2. SIMP 密度插值与疲劳损伤约束的数学耦合2.1 从柔度最小到多约束列式的改写原来的优化问题只有一个约束体积分数 V/V₀ ≤ f。现在要扩成三约束min c(x) U^T K(x) U s.t. V(x)/V0 0.4 sigma_max(x) sigma_allow D(x) D_crit x_min x_e 1, e 1..ND(x)是整个结构的疲劳损伤指标D_crit是允许值。工程上一般取D_crit 1Miner 准则的临界点但髋关节这种不可维修的植入体我建议压到 0.30.5相当于留了 23 倍的安全裕度。这一步最先要确认的是你手上的两篇参考文献到底把损伤定义在单元级还是结构级。如果约束写成max_e D_e D_crit那是单元最坏值约束对优化器的梯度很不友好因为最大值函数不可微如果写成sum_e D_e D_crit那实际上是在约束总损伤量物理意义会偏。主流做法是用 p-norm 把单元损伤聚合成一个光滑函数import numpy as np def pnorm_aggregate(values, p8.0): 把逐单元损伤聚合成可微的标量约束函数。 p 越大越接近 max()但梯度越尖锐工程上 p 取 6~12 比较稳。 values 需为非负加 1e-12 防止 0**p 出现梯度断点。 v np.maximum(np.asarray(values, dtypefloat), 1e-12) return np.sum(v ** p) ** (1.0 / p)这个聚合函数的逻辑是当某一单元的损伤远大于其它单元时它会被 p 次幂放大求和后开 p 次方结果逼近该最大值。参数 p 是调节旋钮——p6 时曲线比较平缓优化器好走但可能压不住局部热点p12 时更贴近真实最大值但灵敏度矩阵里会出现很大的数值需要同步把约束的缩放因子调小。2.2 应力松弛SIMP 里最容易踩的坑SIMP 的刚度插值是E_e E_min x_e^p (E_0 - E_min)p 通常取 3。但如果你直接用同样的插值去算应力会得到一个很反直觉的结果低密度单元x_e 接近 0的应力会趋近于它的真实承载应力除以一个极小的刚度导致sigma_e x_e^p * sigma_e^0这种写法把中间密度区域的应力压得过低优化器就学会了“用一堆半密度单元骗过疲劳约束”。正确的做法是给应力单独一套松弛指数物理量插值形式指数取值说明弹性模量E_e x_e^p · E₀p 3经典 SIMP用于刚度阵组装单元应力σ_e x_e^q · σ_e⁰q 0.5应力松弛避免低密度区伪高应力材料密度ρ_e x_e · ρ₀线性用于质量/体积分数统计疲劳损伤D_e(σ_e)由 S-N 反解损伤是应力的非线性函数不额外松弛q 0.5是 Bruggi 那套松弛方案里最常用的取值。改起来很简单在提取单元应力的地方加一行# x_phys 是过滤投影之后的物理密度 # sigma_vm 是标准有限元解出的单元 von Mises 应力 sigma_relaxed (x_phys ** 0.5) * sigma_vm注意x_phys必须是经过密度过滤和 Heaviside 投影之后的物理密度不是设计变量本身。如果你在优化循环里用的是x而不是x_phys疲劳约束会对过滤半径特别敏感改个 rmin 就要重新调参数。2.3 S-N 曲线与 Goodman 平均应力修正Ti6Al4V 的 S-N 曲线用 Basquin 形式σ_a σ_f · (2N_f)^b反解寿命N_f 0.5 · (σ_a / σ_f)^(1/b)但问题是髋关节的载荷谱是变幅的雨流计数出来的每个循环都带自己的均值 σ_m。必须先把均值修正到对称循环再用 S-N 曲线否则寿命会高估几倍。Goodman 关系式σ_ar σ_a / (1 - σ_m / σ_u)下面是完整的寿命与损伤计算模块import numpy as np # Ti6Al4V 参数占位值实际请以手中两篇文献的取值替换 SIGMA_F 2030.0 # 疲劳强度系数 sigma_fMPa B_EXP -0.104 # 疲劳强度指数 b UTS 950.0 # 抗拉强度 sigma_uMPa def fatigue_life(sigma_a, sigma_m0.0, utsUTS): Goodman 修正 Basquin 反解疲劳寿命返回许用循环数 Nf。 sa np.asarray(sigma_a, dtypefloat) sm np.asarray(sigma_m, dtypefloat) # 平均应力不能超过抗拉强度否则一次性静载就断了 denom np.clip(1.0 - sm / uts, 1e-6, None) sa_eq sa / denom Nf 0.5 * (np.abs(sa_eq) / SIGMA_F) ** (1.0 / B_EXP) # 寿命夹到 [1, 1e12]防止 log 溢出和下界为 0 return np.clip(Nf, 1.0, 1e12) def miner_damage(cycles, utsUTS): Miner 线性累积损伤。cycles [(幅值, 均值, 计数), ...] D 0.0 for amp, mean, cnt in cycles: D cnt / fatigue_life(amp, mean, uts) return DB_EXP是个负数这一点很关键(sa_eq / SIGMA_F)通常小于 1负指数会把它翻到大于 1得到合理的寿命量级。如果误写成正数寿命会小于 1整个约束直接失效。np.clip那一行是为了防止σ_m接近σ_u时除零同时给Nf加上下界——损伤计算里除以 0 会让整个优化循环直接崩掉。2.4 损伤对设计变量的灵敏度链式分解优化器要能用必须给出∂D/∂x_e。按链式法则拆∂D/∂x_e Σ_i (∂D/∂σ_i) · (∂σ_i/∂x_e)后半项从有限元解里拿前半项对 Basquin 加 Goodman 求导后是一个闭式表达式。实际调试时我一般先用有限差分校一遍解析灵敏度误差超过 1% 就说明松弛指数或者物理密度映射那一步有问题。这一步不做后面 MMA 的收敛曲线会抖得没法看。3. 行走/跑步变幅载荷谱构建与雨流计数落地3.1 一个载荷周期的转折点怎么排你手上已经定了四个峰值1107.9 N、1507.9 N、1808.9 N、3433.5 N。这四个值不是四个独立工况而是同一个周期里四个加载阶段的峰值。行走和跑步的区别主要体现在两处一是峰值序列的整体缩放系数二是单位时间内的循环次数步频。要把这四个峰值变成雨流计数能吃的输入得先展开成峰谷交替的时程。缺口处的谷值文献没给常见的处理有两种取前一状态的稳态值或者取峰值的一个固定比例比如 0.10.15 倍。我用第二种import numpy as np def build_cycle(peaks, valley_ratio0.12): 把峰值序列展开成峰谷交替的一维时程。 peaks: 一个周期内的峰值序列单位 N valley_ratio: 峰谷之间的回落比例缺失实验数据时用经验值 返回: [p1, v1, p2, v2, ..., pn, vn] pts [] for p in peaks: pts.append(float(p)) pts.append(float(p) * valley_ratio) return np.asarray(pts, dtypefloat) walk_peaks np.array([1107.9, 1507.9, 1808.9, 3433.5]) run_peaks walk_peaks * 1.35 # 跑步整体放大系数按你的文献取值 sig_walk build_cycle(walk_peaks) sig_run build_cycle(run_peaks)valley_ratio这个参数对损伤结果的影响比想象中大。它决定了每个小循环的幅值而损伤对幅值是指数敏感的——比例从 0.12 改到 0.2累积损伤可能翻一倍。所以这个值要么来自你手上文献里的实测载荷曲线要么在结果里明确标注为假设值。跑步的缩放系数 1.35 也只是占位同样要按文献来。3.2 峰谷提取与四点雨流计数原始时程里相邻点大多是单调的先做极值点压缩把不需要的中间点去掉然后再做雨流计数。ASTM E1049-85 的四点法是标准做法import numpy as np def extract_extrema(series): 提取峰谷点只保留斜率变号的位置 s np.asarray(series, dtypefloat) d np.diff(s) idx np.where(np.diff(np.sign(d)) ! 0)[0] 1 return np.concatenate(([s[0]], s[idx], [s[-1]])) def rainflow_count(ext): 四点雨流计数。ext 为峰谷交替序列。 返回 [(幅值, 均值, 计数)]计数 1.0 为全循环0.5 为半循环。 S list(ext) cycles [] while True: found False for i in range(len(S) - 3): r1 abs(S[i1] - S[i]) r2 abs(S[i2] - S[i1]) r3 abs(S[i3] - S[i2]) # 中间段同时小于左右两段该段构成一个完整循环 if r2 r1 and r2 r3: amp r2 / 2.0 mean (S[i1] S[i2]) / 2.0 cycles.append((amp, mean, 1.0)) del S[i1:i3] found True break if not found: break # 剩下的点两两配对当作半循环处理 for i in range(len(S) - 1): amp abs(S[i1] - S[i]) / 2.0 mean (S[i1] S[i]) / 2.0 cycles.append((amp, mean, 0.5)) return cycles逻辑说明外层while True反复扫描只要发现某三段满足「中间段幅度小于等于左右两段」就把它抽出来作为一个全循环并把中间两个点从序列里删掉然后重新开始扫描。这个「抽掉再扫」的过程就是雨流法的灵魂——不删除的话会重复计数。found标志用来判断这一轮有没有抽到循环抽不到就退出剩下的序列按半循环配对。参数说明amp是幅值等于载荷范围的一半mean是均值等于两端点的算术平均。这两个值后面要分别送进 Goodman 公式和 S-N 曲线。全循环计数记 1.0半循环记 0.5最后累加的时候直接乘进去就行。3.3 从载荷幅值到单元应力幅值这里有个省钱的办法。因为静力线弹性下应力与载荷严格成正比所以不需要对每个载荷峰值都跑一遍有限元。跑一次单位载荷1 N的静态分析得到每个单元的应力影响系数σ_unit然后# sigma_unit: 单位载荷下的单元应力场形状 (n_elem,) def stress_cycles(sigma_unit, cycles, scale1.0): 把载荷循环谱映射成单元应力循环谱 out [] for amp, mean, cnt in cycles: out.append((sigma_unit * amp * scale, sigma_unit * mean * scale, cnt)) return out # 逐单元累积损伤 def element_damage(sigma_unit): D_e np.zeros_like(sigma_unit) for mode_sig in (sig_walk, sig_run): cyc rainflow_count(extract_extrema(mode_sig)) for amp, mean, cnt in cyc: sa sigma_unit * amp sm sigma_unit * mean D_e cnt / fatigue_life(sa, sm) return D_e传入的sigma_unit记得先用前面说的x_phys ** 0.5做一次松弛否则低密度单元的损伤会被系统性地低估。两种工况的循环次数步频 × 服役年限通过cnt的倍数体现出来不用改代码结构。4. 迭代过程的目标/约束数值采集与 Excel 落盘4.1 在优化循环里埋一个记录钩子每次迭代要落四组数柔度、体积分数、最大应力、最大疲劳损伤。麻烦的是scipy.optimize.minimize的callback只传设计变量如果在那儿重新算一遍目标函数和约束等于每次迭代多跑两遍有限元耗时直接翻倍。正确做法是把目标函数和约束函数里的关键中间量缓存下来callback只负责取数落盘import numpy as np import pandas as pd class IterLogger: 缓存式迭代记录器目标/约束函数写缓存callback 只读缓存 def __init__(self, xlsx_pathhip_opt_log.xlsx, sheetiter): self.path xlsx_path self.sheet sheet self.cache {} self.it 0 def update(self, key, value): # 目标函数和约束函数里调用key 用字符串标识 self.cache[key] float(value) def dump(self, xNone): self.it 1 row { iter: self.it, compliance: self.cache.get(compliance), vol_frac: self.cache.get(vol_frac), sigma_max: self.cache.get(sigma_max), damage_max: self.cache.get(damage_max), } if x is not None: row[x_mean] float(np.mean(x)) self._append(row) def _append(self, row): df pd.DataFrame([row]) try: with pd.ExcelWriter(self.path, engineopenpyxl, modea, if_sheet_existsoverlay) as w: if self.sheet in w.book.sheetnames: ws w.book[self.sheet] startrow ws.max_row header False else: startrow 0 header True df.to_excel(w, sheet_nameself.sheet, indexFalse, startrowstartrow, headerheader) except FileNotFoundError: df.to_excel(self.path, sheet_nameself.sheet, indexFalse)update在目标函数里写柔度在体积约束函数里写体积分数在应力约束里写sigma_max在损伤约束里写damage_max。dump由callback调用只做一次字典组装和一次 Excel 追加。注意modea和if_sheet_existsoverlay这两个参数需要 pandas ≥ 1.3 和 openpyxl ≥ 3.0。startrowws.max_row是关键——openpyxl 的行号从 1 开始而 pandas 的startrow从 0 开始所以新数据的起始行正好等于当前最大行号不会覆盖也不会留空行。header只在第一次创建 sheet 时为 True。如果迭代上千次每轮都开一次 Excel 会明显拖慢速度。常见的折中做法是加一个计数器每 10 轮或者每 5 秒 flush 一次def dump_throttled(self, xNone, every10): self.it 1 if self.it % every 0: self._append({...}) # 组装当前记录或者更稳妥把行攒在内存的 list 里注册一个atexit回调做最终落盘这样中途崩了也不会丢数据前提是捕获异常后仍然执行 atexit。4.2 各列的含义与判读口径落盘的表头不是随便定的每一列对应优化器里的一个具体量列名来源量纲判读要点iter计数器—从 1 开始便于和 Origin 的横轴对齐complianceU^T K UN·mm单调下降是正常反弹说明步长太大vol_fracsum(x_phys)/N—必须 ≤ 目标体积分数超了说明该轮不可行sigma_maxp-norm 聚合 σ_eMPa要和sigma_allow手动比对一次damage_maxp-norm 聚合 D_e—目标值 ≤D_crit注意这是聚合值不是真实最大值x_meanmean(x)—长期停在 0.5 附近说明惩罚不足或过滤半径过大damage_max那一列特别容易看错。p-norm 聚合出来的数值比真实最大值略小p 越大越接近所以约束满足不代表真实最坏单元也满足。稳妥的做法是每 20 轮额外算一次真实的max_e D_e写进另一列作为独立校核。5. 疲劳约束下的收敛排错与 Origin 出图技巧5.1 约束震荡的三个常见原因疲劳约束一加进去最常见的就是柔度曲线不再单调下降而是和损伤曲线交替震荡。按我踩过的顺序大概三个原因第一是p值太大。p-norm 的梯度在 p12 以上时对少数单元极其敏感优化器一步就把热点单元的密度拉满下一步又因为体积约束被迫压回去。把 p 从 12 降到 8再给设计变量加一个 0.02 的移动限震荡通常就压住了。第二是灵敏度没做一致性检查。疲劳损伤经过 Goodman 修正和 Basquin 反解之后导数的解析式和有限差分经常对不上尤其是σ_m接近σ_u的那部分单元。花十分钟跑一遍有限差分比调三天参数值。第三是滤波半径和网格尺寸的比例不对。rmin小于 1.5 倍单元尺寸时疲劳约束会在单元尺度上形成高频振荡看起来像棋盘格但不完全是。把rmin提到 23 倍单元尺寸同时把投影的beta从 1 缓慢升到 16震荡一般会消失。5.2 用 Origin 画迭代曲线的一个实用技巧Excel 里导出的数据直接拖进 Origin 有个坑柔度的量级是 10³10⁴损伤的量级是 10⁻³10⁰放在同一张图里损伤基本贴地。不要用双 Y 轴硬叠而是把柔度归一化到第一轮的值C_norm(i) C(i) / C(1) D_norm(i) D(i) / D_crit两列都变成无量纲的 01.x 量级横轴用迭代次数两条曲线叠在一张图上看得很清楚。再加一条y 1.0的水平参考线表示损伤约束边界约束从上方逼近参考线还是穿过参考线一眼就能判断优化是不是真的收敛了。如果要把 STL 结果和曲线对应起来我一般会在 Excel 里多存一列stage取值init、mid、finalOrigin 里用这一列做数据筛选器切换时同时刷新曲线和旁边的结构云图截图。这样汇报的时候不用来回切文件。本文还有配套的精品资源点击获取