★ FEATURED ARTICLE 简介围绕 LS-DYNA 非线性动力学计算中出现的负体积报错这份文档面向已有一定显式分析基础、正在做撞击、大变形或接触仿真的工程人员与研究生帮助其定位计算中断的成因并形成可操作的排查路径。内容从材料参数与状态方程、人工时间步长、模型约束、网格质量等角度梳理负体积的典型来源并给出调整 CONTROL_TIMESTEP 第二个参数如由 0.9 改为 0.7、检查边界条件与接触、修复单元形状等处理思路也提及修改材料后出现速度无穷大的连带现象。资源包共 1 个 docx 文件大小约 32KB篇幅精炼适合作为桌面速查清单或调试记录。已有 1221 人学习下载说明该问题在显式分析实践中较为常见。读者可据此建立“先查约束与网格、再校材料与状态方程、最后调时间步长”的排错顺序减少盲目试算提升模型收敛与结果可靠性。1. 负体积不是网格的错它表示“这个单元已经被压缩到体积为负”压溃算例跑到十几毫秒d3plot 还在正常输出屏幕突然刷出一行Error negative volume in solid element # 38412接着时间步长从 1e-7 掉到 1e-11运气好一点直接 terminate运气差一点节点速度窜到 1e20整场算例变成一堆 NaN。负体积这个报错看着像网格质量问题实际上它只是求解器在告诉你某个单元的雅可比已经反号体积算成了负数再往下算的每一个结果都不可信。这份笔记解决的就是这条从负体积到速度无穷大的链路——材料参数、状态方程、时间步长因子、约束和网格一条条拆开给出可以照着改的卡和可以照着看的输出。它适合手上有 LS-DYNA 算例、正在被负体积反复打断的人也适合刚接手冲击、压溃、跌落这类大变形仿真、还没建立起排查顺序的人。2. 材料参数与状态方程负体积的第一嫌疑怎么逐条排除大变形算例一出负体积很多人的第一反应是回去改网格把局部加密一遍再算一次结果报错换了个单元号接着来。真正排在前面的嫌疑是材料卡和状态方程因为这两处一旦填错整个模型的波速、体积模量和压力响应从一开始就是错的网格划得再漂亮也挡不住。2.1 从 *MAT 卡到状态方程哪几个参数真的决定体积响应体积响应这件事本质上由三个量控制密度 ρ、体积模量 K、以及材料在高压下的状态方程。K 由弹性模量 E 和泊松比 ν 共同决定关系是 K E / (3(1-2ν))。这个式子值得盯一眼当 ν 逼近 0.5 时分母趋近于 0K 发散单元几乎不可压体积应变被强行压到极小压力却要顶着巨大的刚度往上冲数值上非常容易出现压力振荡和负体积。橡胶、泡沫、土壤这几类材料是重灾区。常见的错误做法是把泊松比直接写成 0.5觉得“橡胶本来就不可压”。更稳妥的填法是 0.49 到 0.499或者干脆换用为不可压材料准备的模型比如泡沫类的可压扁泡沫模型它用一条应力-应变曲线描述压实过程不依赖泊松比去撑体积刚度。第二类错误出在状态方程上。带状态方程的材料压力-体积关系由 EOS 参数控制声速 C、Slope 参数 S1/S2/S3、Grüneisen 系数 GAMAO任何一个填错量级μ压缩比为负时算出来的压力符号就是反的单元会被“吸”着往外膨胀然后下一步又反过来来回几下体积就负了。下面这张是常见的弹塑性材料卡字段顺序在大多数发行版里都一样*MAT_PLASTIC_KINEMATIC $# mid ro e pr sigy etan beta 1 7.85e-9 210000.0 0.3 355.0 800.0 0.0这段用的是 mm-ton-s-MPa 单位制密度 7.85e-9 ton/mm³弹性模量 210000 MPa屈服 355 MPa切线模量 800 MPabeta 取 0 表示随动硬化。参数里最容易被忽略的是 etan它决定了塑性段的硬化斜率如果填成和 E 一个量级材料在塑性段几乎不变形应力全部转成压力单元被压穿只是时间问题。我一般会把 etan 控制在 E 的 1% 到 5% 之间先跑一个小模型确认应力应变曲线形态对不对。带状态方程的材料卡要成对出现材料卡管偏应力EOS 卡管压力*MAT_JOHNSON_COOK $# mid ro e pr 3 7.85e-9 210000.0 0.3 $# a b n c m tm tr 355.0 500.0 0.36 0.014 1.0 1793.0 293.0 *EOS_GRUNEISEN $# eos c s1 s2 s3 gamao a 3 4569.0 1.49 0.0 0.0 2.17 0.46C 是材料声速单位跟长度单位一致mm/sS1 是冲击绝热线的斜率GAMAO 是 Grüneisen 系数。这里的坑在于有人把 C 填成 4.569m/s 单位制下的值挪到 mm 里声速小了三个数量级临界时间步长就大了三个数量级第一步算下去单元直接翻面。C 的取值可以用体积声速粗算sqrt(K/ρ) 或者 sqrt(E(1-ν)/((1ν)(1-2ν)ρ))算出来是什么量级填进去就该是什么量级。2.2 一页核对表把密度、模量、泊松比设在同一套单位制里单位制混用是最隐蔽的一类错误因为它不会报错模型照常读入只有结果莫名其妙。下面这张表是常用的三套单位制对照做交叉检查时可以直接对着看物理量mm-ton-s-MPam-kg-s-Pacm-g-μs-Mbar长度mmmcm质量tonkgg时间ssμs力NN1e7 N应力MPaPaMbar钢密度7.85e-978507.85钢弹性模量2100002.1e112.1钢屈服强度3553.55e80.00355判断自己用的是哪套最直接的办法是看密度7.85e-9 是 mm-ton-s7850 是 m-kg-s7.85 是 cm-g-μs。这三者差着六个数量级一旦混进来一个模型的波速就会整体错位。手工翻卡容易漏我一般会用一个脚本把所有材料卡的密度扫一遍先看量级分布import re mat_head re.compile(r^\*MAT_[A-Z_], re.I) def scan_density(path): 扫描 k 文件里每张材料卡的第一个数据行取出密度字段 with open(path, encodingutf-8, errorsignore) as f: lines f.read().splitlines() result [] i 0 while i len(lines): if mat_head.match(lines[i]): name lines[i].strip() j i 1 # 跳过 $ 注释行和空行它们不是数据 while j len(lines) and (lines[j].startswith($) or not lines[j].strip()): j 1 if j len(lines): fields lines[j].split() if len(fields) 2: try: result.append((name, float(fields[1]))) except ValueError: pass i j i 1 return result for name, ro in scan_density(model.k): print(f{name:26s} RO {ro:.3e})脚本的逻辑很直白遇到*MAT_打头就认为进了材料卡跳过紧随其后的$#注释行取第一个数据行的第二个字段当密度。之所以跳过注释行是因为 LS-DYNA 的 k 文件里$#开头的那一行只描述字段名不参与计算如果直接读第二行会拿到中文字符串。errorsignore是为了兼容某些在注释里写了中文、编码不是 UTF-8 的旧文件。跑完之后看输出mm-ton-s 模型里的密度应该全部落在 1e-10 到 1e-8 之间如果冒出一个 7.85 或者 7850基本就能锁定问题卡片。参数换一下路径就能扫别的模型脚本不用改。2.3 用 matsum 与 d3hsp 反向验证材料有没有被正确处理改完卡别急着信自己让求解器回话。第一次运行时把求解时间设到极短比如 1e-5 秒只让它做初始化和一个增量步看输出文件里的质量、体积、初始时间步长是不是符合预期。# 先跑一个只做初始化的短算例确认材料被正确读入 lsdyna imodel.k jobidchk ncpu4 memory200m这里 jobid 决定输出文件前缀ncpu 是并行核数memory 给的是每核内存。跑完重点看三个地方d3hsp 里搜total mass把总质量和手算的体积乘密度对一下搜initial time step看初始步长是不是和你估的 L/c 一个量级matsum 文件里看每个 part 的内能和体积是不是单调合理。如果总质量差了 1000 倍不用往下算了回去查单位制。3. 时间步长因子 TSSFAC把 0.9 改成 0.7 改的到底是什么材料查完一圈没问题负体积还在下一个嫌疑就是时间步长。*CONTROL_TIMESTEP第二个字段从 0.9 改成 0.7是很多现场记录里被反复提到的一招它确实有用但得先知道它改的是稳定裕度不是物理模型。3.1 临界时间步长的算法与 TSSFAC 的位置显式求解的时间步长按单元逐个算取全局最小值dt_e TSSFAC × L_e / c。L_e 是单元的特征长度对实体单元来说约等于体积除以最大面的面积c 是材料声速。TSSFAC 就是那个安全系数理论稳定极限在 1.0 附近实际取值要打折实体单元默认 0.9。0.9 这个值在多数算例里够用但遇到几类情况就会紧张单元长宽比大、材料声速高、单元在变形过程中特征长度塌缩得很快。一旦某个单元的实际步长超过稳定极限该单元的体积更新就会越界先振荡后反号最终报负体积。把 TSSFAC 降到 0.7等于给每个增量步多留了 22% 的裕度通常能把这类“临界失稳”压回去。代价也很直接步长小了步数就多。0.9 除以 0.7 约等于 1.286同一个算例的增量步数大约增加 29%墙钟时间大致按这个比例往上走。所以这一招不是无条件好用它适合“就差一点点”的算例不适合拿来掩盖一个已经错得很离谱的模型。3.2 一次可复现的 0.9→0.7 修改与前后对比改法本身没什么技术含量关键是卡片的完整字段要摆对别把参数写到隔壁字段去*CONTROL_TIMESTEP $# dtinit tssfac isdo tslimt dt2ms lctm erode ms1st 0.0 0.7 0.0 0.0 0.0 0.0 0.0 0.0字段逐个说清楚dtinit 是初始步长填 0 表示让程序自己算tssfac 就是这次要改的安全系数isdo 控制是否对壳单元用不同的步长算法0 是默认tslimt 是最小步长下限dt2ms 是质量缩放的目标步长填 0 表示不做质量缩放lctm 用曲线控制步长一般留 0erode 是当步长低于 tslimt 时是否删除单元0 是不删ms1st 控制质量缩放是否只在第一步做。排查负体积的阶段dt2ms 和 erode 都建议保持 0先把问题暴露出来别让程序替你悄悄删单元。改完跑一遍用同样的命令只换 jobidlsdyna imodel.k jobidrun_t07 ncpu8 memory300m对比要看的是三件事d3hsp 里的初始时间步长是否按比例缩小了约 22%glstat 里步长随时间的曲线是否还在合理区间没有出现连续的台阶式下跌负体积报错是否消失或者至少推迟到了更大的变形时刻。如果报错只是推迟、没有消失说明 TSSFAC 只是在拖时间真正的问题还在材料或者网格上这时候应该回到上一章而不是继续把 TSSFAC 往下压到 0.5。3.3 质量缩放、DT2MS 与负体积之间的相互影响这一段值得单独拿出来说因为它正好能解释“改完材料以后速度变成无穷大”这个现象。质量缩放的作用是给步长过小的单元人为加密度把 dt 顶到目标值。dt2ms 填正值时所有单元都加填负值时只给步长小于 |dt2ms| 的单元加。加质量能让步长回到目标但它同时改变了惯性动量等于质量乘速度质量被放大几十倍同样的冲量下去速度响应就完全不是原来那个物理过程了。把 TSSFAC 从 0.9 改到 0.7 之后全局步长变小了如果同时开着 dt2ms为了把步长顶回目标值程序需要加更多的质量added mass 比例会明显上升。某些畸变单元被加了巨量质量之后一旦受到接触力或者局部压力速度会瞬间冲高表现出来就是节点速度 1e20 或者 NaN。所以顺序应该是先关掉质量缩放用真实的步长把负体积问题解决掉确认算例能稳定跑通再打开 dt2ms 做加速并且盯住 glstat 里的 added mass实体单元一般控制在 1% 到 2% 以内。提示排查阶段把 erode 保持为 0。开着 erode 虽然能让算例跑完但被删掉的单元带走了质量和能量最后的支反力和吸能曲线全是错的这种“跑通”比报错更有害。4. 约束不足与速度无穷大从负体积到 NaN 的连锁反应材料和时间步长都过了一遍还有一类负体积是模型自己“飘”出来的。这类算例的特征很明显报错的单元往往远离加载区位置在模型角落或者某个没有被任何接触覆盖的自由面上而且经常伴随速度无穷大一起出现。4.1 无约束自由度如何使速度增大到无穷大显式求解每一步都在解 F ma。如果一个节点既没有边界约束也没有接触力或者外力作用它的加速度理论上为零速度应该保持不变。问题出在数值上初始穿透、材料初始化应力、接触搜索的微小不平衡都会给这个节点一个极小的残余力。单步看不出来但算例跑到几万步之后这个残余力累积出来的速度就可能不可忽略而且没有任何机制把它拉回来。更麻烦的是它和负体积的耦合。一个自由飘的节点所在的单元被相邻单元拖着走形状逐渐歪掉特征长度变小局部时间步长跟着变小。步长一小同样的残余力在更多步里累积速度进一步变大单元变形更快形成一个正反馈。等你在 d3hsp 里看到速度量级不对劲的时候单元早就负体积了。第二类来源是刚体。刚体材料卡上的 CMO、CON1、CON2 三个开关如果不设成全约束刚体自己可以转、可以平移整个刚体就像一个没有支点的锤子砸到哪算哪。这种情况在包含多个刚体部件的装配体算例里特别常见。*MAT_RIGID $# mid ro e pr n m alias 2 7.85e-9 210000.0 0.3 0.0 0.0 0.0 $# cmo con1 con2 1.0 1.0 1.0CMO 等于 1 表示约束质心的平动和转动CON1 等于 1 表示约束全局 x 方向平动CON2 等于 1 表示约束 y、z 方向平动和三个方向的转动。三个都置 1刚体就只剩你显式给它的运动不会自己乱跑。如果某个刚体确实需要运动就用 *BOUNDARY_PRESCRIBED_MOTION 或者载荷曲线去驱动而不是靠留自由度让它自己动。4.2 边界条件、接触与刚体的检查顺序检查要有固定顺序否则容易在同一个地方反复打转。我一般的顺序是先看有没有完全自由的节点再看边界条件覆盖对不对然后看刚体约束最后看接触。第一步在 LS-PrePost 里用 Model 下的检查工具找 free node它会列出既没有被约束、也没有被任何接触面覆盖的节点。这一步能捞出一大半问题。第二步对着 *BOUNDARY_SPC_SET 看节点集。常见错误是节点集只选了表面的节点而内部或者边角上还有一层没被选进去表面上看起来约束住了实际上模型还能沿某个方向滑。*SET_NODE_LIST_TITLE floor_fixed $# sid da1 da2 da3 da4 1 0.0 0.0 0.0 0.0 $# nodes 1 2 3 17 18 *BOUNDARY_SPC_SET $# nsid cid dofx dofy dofz dofrx dofry dofrz 1 0 1 1 1 1 1 1dofx 到 dofrz 六个数1 是约束0 是放开。做压溃或者跌落的时候底面通常全约束侧面要看是否用无反射边界或者对称边界。有人图省事只约束了 z 方向x、y 两个方向放着不管短时间算例看不出来一旦单元开始大变形那点残余力就够把模型推出去了。第三步看接触。自接触算法里的初始穿透检查一定要打开初始穿透会在第一步就给单元一个巨大的法向力直接把单元压翻。接触刚度系数也别用默认值硬扛软材料和硬材料混搭的模型材料刚度差两个数量级以上时默认算法容易在软的那一侧压穿。*CONTROL_CONTACT $# slsfac rwpnal islchk shlthk penopt thkopt snlog 0.1 1.0 1 1 1 1 0islchk 是初始穿透检查1 表示检查并处理penopt 控制穿透的处理方式snlog 控制是否在 d3hsp 里输出接触搜索信息排查阶段建议打开能直接看到初始穿透发生在哪几个节点。第四步才是看刚体约束顺序反了容易把时间浪费在一个本来就没问题的刚体上。4.3 用 glstat、d3hsp 与单元历史变量定位第一个关键单元模型改完还是要靠输出来验证。设几个数据库输出把过程记录下来*DATABASE_GLSTAT $# dt binary lcur ioopt 1.0e-4 0 1 *DATABASE_HISTORY_SOLID $# id1 id2 id3 id4 38412 38413 38414 38415 *DATABASE_BINARY_D3PLOT $# dt lcdt beam npltc 5.0e-4 0 0 0glstat 记录全局能量和步长dt 给 1e-4 足够密DATABASE_HISTORY_SOLID 把报错单元和它周围几个邻居的历史变量写进 elout 或者 nodout之后可以精确看到哪个时刻它的体积开始掉d3plot 的输出间隔别太大否则回看时找不到出事前最后一个正常状态。拿到输出之后用脚本先筛一遍步长什么时候开始掉def find_dt_drop(path, drop_ratio0.3): glstat 前两列是时间和时间步长筛出步长突然下降的时刻 prev None with open(path, encodingutf-8, errorsignore) as f: for line in f: parts line.split() # 跳过表头、注释和不足两列的行 if len(parts) 2 or not parts[0][0].isdigit() and parts[0][0] ! -: continue try: t, dt float(parts[0]), float(parts[1]) except ValueError: continue if prev and dt prev[1] * drop_ratio: print(ft {t:.6e} dt 从 {prev[1]:.3e} 掉到 {dt:.3e}) prev (t, dt) find_dt_drop(glstat)这里假设 glstat 的第一列是时间、第二列是时间步长多数发行版都是这个顺序如果你的输出列序不同打开文件看一眼表头把下标改掉就行。drop_ratio 设为 0.3 表示相邻两次输出之间步长掉到 30% 以下就报警这个阈值调小一点能更早发现异常但噪声也多。找到时刻之后回到 d3plot 看那一步的单元形态通常能直接看出是哪个位置的单元被压穿是接触区、约束边缘还是自由面。定位到位置再回去看那个位置的边界条件和材料比盲改参数快得多。5. 负体积排查避坑五个曾经的失败案例排查顺序讲完剩下的都是具体操作中的细节问题。下面五条是同一类问题反复出现的地方每条按现象、原因、解决三步写能对上号的直接照着改就行。5.1 材料和单位制的两个坑坑一算例第一步就报负体积。现象是求解器刚初始化完还没到第一个 d3plot 输出时刻就 terminated3hsp 里的 initial time step 比预期大三个数量级。原因基本可以锁定为弹性模量单位填错比如本该填 210000MPa结果填了 210GPa 的数值搬到了 MPa 里。模量小了一千倍声速小约 31 倍临界步长大 31 倍——等等这里我算错了重新捋一遍声速正比于 sqrt(E/ρ)E 小 1000 倍则声速小约 31.6 倍临界步长 dt L/c 变大 31.6 倍。步长大三十多倍远超稳定极限第一步就把单元压翻。解决办法是把所有材料卡的 E、屈服、切线模量统一到 MPa改完用 2.2 那个脚本扫一遍密度再手算一个单元的 L/c 和 d3hsp 里的初始步长对一下量级对不上就说明还有别的卡混了单位制。坑二压力异常大、单元体积几乎不变但应力爆炸。现象是应力云图上一片通红数值远超材料强度但变形很小。原因是泊松比填了 0.5体积模量 K E/(3(1-2ν)) 分母为零程序内部要么取了一个极大值要么直接溢出了。这个坑在橡胶密封件、缓冲垫这类模型里出现频率很高。解决办法是把 ν 改到 0.49 到 0.499或者换成专门处理不可压材料的模型用应力-应变曲线描述压实不要靠泊松比硬顶。改完记得回看压力云图的量级正常橡胶件的静水压力应该是几 MPa 到几十 MPa不是几万。5.2 TSSFAC 改小与质量缩放的两个坑坑三负体积消失了速度变成无穷大。现象是把 TSSFAC 从 0.9 改到 0.7 之后负体积报错确实没了但跑到中途出现节点速度 1e20能量曲线直接垂直上升。原因是同时开着 dt2ms 质量缩放全局步长被压小之后为了让步长顶回目标值程序给畸变单元加的质量比原来更多added mass 比例可能从几个百分点涨到几十个百分点惯性被严重放大。遇到这个现象先看 glstat 里的 added mass 那一列超过 5% 就说明质量缩放参与过度。解决办法是排查阶段把 dt2ms 设为 0用真实步长跑通确认物理过程合理之后再把 dt2ms 设成一个负值只给步长过小的单元加质量并且保持 added mass 在 1% 到 2% 以内。坑四TSSFAC 一路往下压报错只推迟不消失。现象是从 0.9 压到 0.5负体积出现时刻从 3 ms 推迟到 8 ms但终究还是会来。原因是真正的问题在单元畸变或者接触穿透TSSFAC 只是把失稳往后拖。这种情况下继续压系数没有意义再往下压计算成本翻倍报错照样出现。正确做法是回到 d3plot看 8 ms 附近是哪个单元、什么位置、周围有没有接触。我一般会把这段时间的 d3plot 输出间隔改到 1e-5 秒逐帧看单元形态变化通常能发现是接触边缘的单元被反复挤压或者是某个尖角处的单元被压成了薄片。5.3 约束、接触与单元类型的一个坑坑五报错单元总在模型角落里而且换一个算例还在类似位置。现象是负体积的单元号每次都不同但位置都在远离加载区的边角速度无穷大也出现在那一带。原因是有自由节点没被约束住加上边角处单元本身长宽比偏大特征长度小、步长小正反馈最先从这里起爆。解决办法是先在 LS-PrePost 里用 free node 检查把所有完全自由的节点列出来逐个确认是不是漏选了节点集然后检查刚体卡的 CMO、CON1、CON2 是不是都没设最后看接触面有没有覆盖到那个角落自接触的搜索范围是不是被某个阈值卡掉了。顺带说一句单元类型。四面体单元在大变形下特别容易负体积因为它的四个面在压扁过程中很快退化雅可比先接近零再反号。如果你的模型里大量用的是四面体又恰好是最需要大变形的区域那负体积几乎是一定会来的。这类区域要么换成六面体要么把 *CONTROL_SOLID 的 ESORT 打开让程序在单元畸变时自动重划*CONTROL_SOLID $# esort fmatrx niptets swlocl psfail t10jtol icohprs 1 0 1 0 0 0 0esort 置 1 只对四面体有效它会在单元质量变差时按角度准则重新划分局部网格代价是计算量上升而且重划后的网格和原来的节点编号对不上做节点历史变量追踪时要留意。六面体单元不受这个参数影响所以如果你的模型以六面体为主这条不用设。6. 把负体积提前检查并处理掉一份上线前的自检表和一个取舍改到能跑通只是及格线能不能稳定跑完、结果可不可信是另一回事。我现在每次提交长算例之前会强制走一遍下面这张表前后大概十分钟能省掉一晚上的重算。检查项判断标准用什么看单位制一致性所有材料卡的密度量级相同密度扫描脚本泊松比实体材料 ≤ 0.499材料卡逐条过初始时间步长与手算 L/c 同量级d3hsp 搜 initial time step自由节点无完全无约束、无接触的节点LS-PrePost free node 检查初始穿透d3hsp 无大量初始穿透警告打开接触搜索日志added mass全程 2%glstat 对应列这张表里最容易被跳过的是初始穿透那一条。有人觉得穿透一点点没关系程序会自动推出去。实际上一旦初始穿透发生在薄壁或者尖角处第一步给的推开力就足以把单元压翻后面所有关于材料和步长的调整都是白费。单元算法这块也有取舍。六面体单元ELFORM 取 -1 或 -2抗畸变能力强大变形下不容易负体积但前处理划网格的时间可能是四面体的好几倍。四面体单元划起来快配合 ESORT 也能跑代价是每一步都要检查重划计算时间上去了而且重划区域的应力会有轻微不连续。我的习惯是主变形区一定用六面体过渡区和远离变形区的地方可以用四面体加 ESORT两边的界面用绑定接触连起来别指望共用节点能扛住大变形。还有一个后悔药*CONTROL_TERMINATION 里的 dtmin 别设得太小设成 1e-9 这种量级等于允许程序在步长崩塌之后继续磨几个小时才退出白白占着机器。设成正常步长的百分之一左右一旦步长掉下去就让它尽早停你也能早点拿到报错信息去定位。从那以后我每次提交长算例之前都会先跑一个只算初始化的短算例把 d3hsp 里的总质量、初始步长和手算值对一遍再开正式计算。这个习惯救过我很多次也让我少熬了不少夜。希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
