企业官网建设流程全解析

大多数人第一次要在ABAQUS里做复合材料冲击损伤仿真，第一反应都是直接用内置的Hashin模型。但真把低速冲击模型跑起来就会发现，内置模型在很多场景下不够灵活：损伤变量不能自定义，没法输出自定义的失效状态，也没法把应变率、温度这些因素耦合进来。所以VUMAT子程序几乎成了绕不开的一步。

这篇文章围绕铺层复合材料的冲击损伤仿真，把基于三维Hashin准则的VUMAT子程序开发从理论到落地的完整路径拆开讲。从Hashin公式里每个分量到底在算什么，到Fortran代码怎么组织，再到冲击模型里材料方向、接触、网格这些前置条件，最后是调试手段和结果验证。适合正在做低速冲击、落锤冲击、弹道冲击仿真，或者准备用子程序做自定义损伤模型的研究生和工程师参考。内容尽量说人话，尽量给可以直接抄作业的东西。

1. Hashin准则的三维修订：先搞清楚每种失效模式到底在算什么

1.1 四种失效模式和对应的判定公式

Hashin准则把单向复合材料层内的失效分为纤维拉伸、纤维压缩、基体拉伸、基体压缩四种模式。冲击工况下应力状态通常是三维的，所以必须用三维形式，不能用平面应力简化。这里先给出最常用的三维Hashin判定公式，后面所有代码都按这套公式写。

纤维拉伸失效（σ11 ≥ 0时触发）：

F_ft = (σ11 / X_T)^2 + α * (τ12 / S_L)^2 + (τ13 / S_L)^2 ≥ 1

纤维压缩失效（σ11 < 0时触发）：

F_fc = (σ11 / X_C)^2 ≥ 1

基体拉伸失效（σ22 + σ33 ≥ 0时触发）：

F_mt = ((σ22 + σ33) / Y_T)^2 + (τ23^2 - σ22 * σ33) / S_T^2 + (τ12 / S_L)^2 + (τ13 / S_L)^2 ≥ 1

基体压缩失效（σ22 + σ33 < 0时触发）：

F_mc = ((σ22 + σ33) / (2 * S_T))^2 + ((Y_C / (2 * S_T))^2 - 1) * (σ22 + σ33) / Y_C + (τ23^2 - σ22 * σ33) / S_T^2 + (τ12 / S_L)^2 + (τ13 / S_L)^2 ≥ 1

各符号的含义如下：X_T是纵向拉伸强度，X_C是纵向压缩强度，Y_T是横向拉伸强度，Y_C是横向压缩强度，S_L是纵向剪切强度，S_T是横向剪切强度。α是剪切应力贡献系数，常见取1.0，也有人按文献取0到1之间的值，目的是让纤维拉伸的预测偏保守还是偏激进，后者只是工程折中，没有严格的物理推导。

有两点容易被忽略。第一，公式里的σ11、σ22、σ33是材料局部坐标系下的应力分量。单向复合材料1方向是纤维方向，2方向和3方向是面内横向和厚度方向。如果你的模型材料方向没指对，Hashin公式计算出来的全是错的。第二，τ23^2 - σ22 * σ33这一项很有意思——如果σ22和σ33同号且数值很大，这一项会变小甚至为负，物理含义是横向压应力抑制了基体剪切开裂的倾向。所以千万不要把这个括号拆掉，也不要误写成τ23^2 + σ22 * σ33。

1.2 三维和二维的区别，以及为什么冲击必须用三维

二维Hashin是基于平面应力假设的，即σ33 = τ23 = τ13 = 0。这个假设在薄板面内受载时基本够用，但低速冲击过程中，接触区域下方有显著的厚度方向压应力，冲击背面则会出现明显的弯曲正应力，平面应力假设与真实状态出入很大。三维Hashin保留了σ33和两个横向剪切分量，对局部应力状态的描述更接近真实物理。

从代码角度说，二维和三维的差异不仅是多两个应力分量的问题，还牵涉弹性矩阵的维度。三维正交各向异性本构有9个独立工程常数：E11、E22、E33、ν12、ν13、ν23、G12、G13、G23。初学VUMAT的人最容易在这块出问题——只给了六个常数，刚度矩阵里G13、G23不知道填什么。有些材料厂商给的数据表里没有G13和G23，可以用近似关系G12 ≈ G13，或者用半经验公式估算，但一定要在文档里写清楚假设，否则后面审稿或验收时说不清。

1.3 应力更新中的符号约定和坐标系统一

在VUMAT里，ABAQUS传入的是应变增量strainInc，每个增量步的应力更新逻辑是：

σ_new = σ_old + C_d : Δε

这里C_d是当前损伤状态下折减后的刚度矩阵。VUMAT中应力张量和应变张量的存储顺序是固定的：直接分量在前（11、22、33），剪应力分量在后（12、13、23）。在显式分析中，这个顺序和坐标系的旋转变换由ABAQUS内部完成，子程序里拿到的应力分量已经是材料局部坐标系（由ORIENTATION决定）下的值。

这意味着在子程序里不需要自己做坐标变换，除非你的材料定义在全局坐标系下而铺层是曲面的。记住这一点，可以省掉大量调bug的时间。

2. VUMAT子程序的骨架搭建：从接口变量到损伤状态变量的组织

2.1 VUMAT接口中真正需要关心的变量

VUMAT的官方接口参数很多，一眼看过去容易吓住。实际编程时，读入侧重点关注这几个：nblock是向量化块大小，通常理解为一批材料点；nprops是材料参数个数；props是材料常数组；stateOld和stateNew分别是旧状态变量和新状态变量；strainInc是应变增量；stressOld是要被更新成stressNew的旧应力。写出行侧主要是stressNew和stateNew。

代码结构的大致框架如下：

subroutine vumat( & nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal, & stepTime, totalTime, dt, cmname, coordMp, charLength, & matDes, stateOld, stateNew, fieldOld, fieldNew, & stressNew, stateNew, fieldNew, intEneNew, sseNew, & spdNew, svNew, stressOld, stateOld, strainInc, relSpinInc, & tempOld, tempNew, dtOld, dtNew, snew, cnew, rhoOld, & rhoNew, engOld, engNew, plasticStrain, plasticStrainOld, & tangent, strainOld, sseOld, spdOld, svOld ) C include 'vaba_param.inc' C character*80 cmname dimension props(nprops), stateOld(nstatev, nblock), & stateNew(nstatev, nblock), stressOld(nblock, ndir+nshr), & stressNew(nblock, ndir+nshr), strainInc(nblock, ndir+nshr) C do k = 1, nblock ! 主循环体 end do return end

主循环体里的逻辑一定要按增量步组织。显式分析本身的增量步很小，所以直接用线性增量格式更新应力即可。不要在VUMAT里尝试做Newton-Raphson迭代，那是有意给自己找麻烦。显式程序的稳定时间步通常都在纳秒到微秒量级，用增量形式误差可控。

2.2 状态变量规划：后处理和调试都靠它

状态变量的设计在动手写代码前就要定下来，不然后面后处理的时候一个个试SDV编号，心态会崩。我推荐一个比较通用的规划：

• SDV1：纤维拉伸失效标志（0或1）
• SDV2：纤维压缩失效标志（0或1）
• SDV3：基体拉伸失效标志（0或1）
• SDV4：基体压缩失效标志（0或1）
• SDV5：纤维损伤变量df（0到1连续值，取历史最大值）
• SDV6：基体损伤变量dm（0到1连续值，取历史最大值）
• SDV7：当前的损伤模式编号（辅助调试用）

这里有一个非常关键但新手容易犯的错：失效标志和损伤变量是两回事。失效标志是Hashin准则是否被触发的布尔量，损伤变量是刚度退化程度的连续量。一个单元在某个增量步触发了Hashin准则，不代表它立刻失去全部刚度。如果直接把失效标志当损伤变量，损伤为1时刚度直接清零，极容易导致冲击局部产生严重的应力振荡，甚至单元畸变到计算崩掉。

损伤演化应该让刚度逐渐退化。最简单的做法是写成指数退化形式：

df_new = max(df_old, 1.0 - exp(-β * (F_ft - 1.0)))

β是退化速率参数，控制损伤从0增长到1的快慢。这个形式的好处是产出的应力-位移曲线比较平滑，不容易出现数值跳变。退化函数的物理意义不深究，但工程上非常实用。如果不引入断裂能，至少用这种方法把曲线抹平。

2.3 三维正交各向异性刚度的折减逻辑

无损伤时，三维正交各向异性材料在局部坐标系下的柔度矩阵可以写成工程常数的形式，再求逆得到刚度矩阵。损伤发生后，对刚度矩阵做折减。常用的折减策略是：纤维损伤df主要折减C11以及涉及纤维方向的剪切项（C12、C13、C55、C66），基体损伤dm折减其余方向（C22、C33、C23、C44）。

具体实现时，可以先组装无损伤刚度矩阵C0，然后构建折减因子矩阵。

d11 = 1.0 - df d22 = 1.0 - dm d12 = sqrt(d11 * d22) d13 = sqrt(d11 * d22) d23 = d22 d44 = d22 d55 = sqrt(d11 * d22) d66 = sqrt(d11 * d22)

这种方法比统一乘一个(1-d)要合理得多，因为基体损伤对纵向模量的影响和对横向模量的影响本来就不一样。如果四种失效模式对应四个独立损伤变量（df_t、df_c、dm_t、dm_c），折减逻辑会更复杂，但思路完全一致：区分哪些刚度分量受哪种损伤影响。从工程验证角度，两变量的方案已经能覆盖绝大多数低速冲击工况，先跑通再说。

3. 冲击模型里真正的隐形杀手：材料方向、接触设置和网格尺寸

3.1 材料方向指错，子程序再对也是白写

这是我见过最常见的“VUMAT算出来全错”的原因，而且这个错误极其隐蔽。ABAQUS里铺层复合材料通常用ORIENTATION定义局部坐标系，然后用SOLID SECTION把某个铺层区域关联到对应的坐标系。典型操作是每个铺层角度的单元都单独建一个Solid Section，并在其中指定ORIENTATION名字。

这里必须做到两件事。第一，所有铺层必须使用三维坐标系。不要用默认的全局坐标系硬充，除非你的铺层真的全是0度，否则结果必错。第二，检查坐标系的箭头方向。后处理里可以在Visualization模块下打开Material Orientation视图，直接看每个单元的1方向箭头是否沿着纤维方向。我在项目里遇到过0度层和90度层的单元局部坐标看起来对了，但45度层和-45度层的1方向箭头完全反了的情况，原因是定义ORIENTATION时旋转轴的顺序写错。检查这一步花两分钟，省下来的是几个晚上的排查时间。

另外，如果你的几何体是曲面或带角度的变厚度结构，建议用离散坐标系（discrete）而不是解析坐标系。离散坐标系可以精确贴合曲面形态，尤其是在曲率大的接触区域，解析坐标系很容易让纤维方向跟着法向跑偏。

3.2 接触、沙漏和单元删除的配合

低速冲击仿真里接触设置看似简单，实则影响结果稳定性。推荐直接用通用接触（General Contact），接触域包括所有外表面，冲击物和层合板之间用罚函数法，摩擦系数按试验工况取0到0.5之间。如果只关心冲击力峰值，摩擦影响有限；如果关心试件的整体变形模式和分层面积，摩擦系数值得标定。

层合板单元大多数情况下用C3D8R。这个单元是减缩积分，沙漏问题躲不掉。在冲击接触区域，网格细化之后沙漏模式容易被激发，必须在Section Controls里开启增强沙漏控制（Enhanced hourglass control），同时在后处理里监测ALLAE（伪应变能）占总内能的比例。稳妥的经验值是沙漏能控制在5%以内，超过10%说明网格或控制参数有严重问题。

单元删除也要提前想好。如果VUMAT里把某方向损伤变量推到接近1，下一步刚度几乎为零，单元可能产生极端畸变。打开ELEMENT DELETION = YES可以缓解，但这里的坑是：一旦删除单元，接触面上的压力分布会突变，冲击力曲线可能出现锯齿状振荡。不删除单元，畸变单元又会拉低稳定时间步。实际做法是在单元主要承载方向的所有损伤变量都超过阈值时才删除，不要损伤一触发就删。

3.3 网格尺寸、时间步和质量缩放

冲击仿真里网格尺寸和时间步是互相牵制的。显式分析的稳定时间步大致正比于最小单元边长除以材料波速。复合材料中碳纤维增强树脂的面内波速很高，如果接触区网格细到0.5mm，稳定时间步可能只有纳秒级别，计算量会非常可观。

网格方案建议这样：冲击接触区域（以接触点为中心，半径大约是冲击头直径的1.5到2倍区域）加密到0.5到1mm，远离区域逐步过渡到2到4mm。厚度方向每个铺层至少切一个单元，如果铺层较厚或者想捕捉弯曲应力梯度，建议一个铺层切两个单元。单元长宽比不能太夸张，冲击区域附近的单元长宽比在3以内比较保险。

质量缩放几乎必用，但必须控制住。显式分析里常用的方式是固定时间增量步，让ABAQUS自动调整密度。判断标准是看模型总动能占内能的比例——低速冲击的动能变化主要是冲击物的动能转移，如果因为质量缩放导致整个试件的动能比重明显异常，结果就不可信了。我在工程判据上习惯把质量增加控制在5%以内，也就是在Step里设置Time Scaling Factor时，增量和无质量缩放时真实稳定时间步的比值不要拉得太大。

4. 单单元测试与代码调试：把排查链路的每一步走通

4.1 为什么一定要做单单元测试

很多人写完VUMAT第一件事就是直接跑整板冲击模型，然后对着满屏的报错发呆。这就像没测过发动机就装整车——完全没有必要。正确的路径是先建立一个包含一个C3D8R单元的模型，分别做几个基本工况：纵向单轴拉伸（验证纤维拉伸）、横向单轴拉伸（验证基体拉伸）、横向压缩（验证基体压缩）、纯剪切（验证剪应力项）。

单单元测试的操作要点是：一个单元，一个ORIENTATION与全局方向一致，底面固定，顶面给位移边界条件，用光滑的幅值曲线加载，避免瞬时冲击产生应力波。输出该单元的正应力分量和SDV值。比如纵向拉伸工况，理论上S11应该随着施加的拉伸应变线性增长，增长到某个临界点时触发SDV1，然后S11逐渐下降或维持在一个平台附近，具体形态取决于损伤演化方式。

这里有一个新手很容易踩的坑：直接把单元某个方向的应变加到很大，结果S11已经远超强度阈值，但SDV1始终是0。原因通常在材料局部坐标系下该方向并不是1方向，或者你在代码里误把σ11写成了应力张量的第几个分量。建议在子程序里加debug输出，把每个KD的应力分量打印出来，对照ABAQUS自带的单元输出S11、S22、S12，很快就能定位。

4.2 常见报错和异常行为的排查链路

先列一个常见异常的排查顺序表：

关于NaN问题多说一句。在损伤变量接近1时，1-df趋近于0，刚度矩阵里对应行几乎为零。这时候如果应力更新计算里还有除以该项的表达式，很容易得到无穷大。解决办法是给损伤变量设一个上限，比如0.999，而不是让它精确等于1。同时在步进更新后判断一次，如果某个应力分量绝对值超过物理上限（比如超过了材料模量乘以特征应变几个量级），就强制修正并输出一条警告信息。

实际调试时，用WRITE打印非常有效。很多人不敢在VUMAT里打印，担心输出太多影响性能。在排查阶段完全无所谓，去掉或加开关都行。建议用INT格式加状态变量一起打印：

if (k <= 3) then write(6,*) 'STEP=', totalTime, 'K=', k, & 'S11=', stressNew(k,1), 'SDV1=', stateNew(1,k) endif

这样能实时观察最开始几个材料点的状态变化，不用等整个作业跑完再分析。

4.3 用内置Hashin模型当基准结果

ABAQUS自带了一个基于Hashin起始准则和损伤演化的复合材料损伤模型。这个内置模型在层内损伤的预测上可以作为VUMAT的对照基准。做法是建两个一模一样的冲击模型，一个用内置模型，一个用你自己的VUMAT，材料参数保持一致。

对比时重点关注三条曲线：冲击力-时间曲线、冲击力-位移曲线、总吸收能量。VUMAT结果和内置模型有差异是正常的，因为损伤演化方式不同，但只要峰值冲击力和能量吸收差别在合理范围（比如20%以内），说明你的主逻辑方向是对的。如果差了一个量级，回去查本构参数和折减逻辑。这个方法在毕业设计和项目验收阶段都非常好用，能快速证明你的子程序“整体合理”。

5. 结果验证与后处理：如何用损伤云图和能量曲线反推代码正确性

5.1 SDV输出和后处理设置

做完计算之后，判断VUMAT是否正确的重要依据就是SDV云图。在Step输出设置里记得选State Variables，否则后处理看不到SDV。对应上面的SDV规划，SDV1到SDV4分别显示四种失效标志的分布区域。典型的低速冲击损伤形态是：接触面下方以基体压缩/剪切损伤为主，冲击背面以基体拉伸损伤为主，纤维断裂主要出现在背面高弯曲拉应力区域。如果你的SDV云图形态完全不符合这些规律，比如纤维损伤在接触面正下方大面积出现，多半是材料方向或铺层方向设置有问题。

更细致的验证可以看损伤变量SDV5和SDV6的分布。失效标志只能看出哪些区域失效，损伤变量能看出失效的严重程度梯度。冲击区域中心和外周的损伤变量应该有明显递减梯度，而瞬时退化模式下这个梯度往往非常生硬，过渡区特别窄。这种生硬梯度虽然反馈了损伤的局部性，但和实际断口形貌的连续分布有差距，这也是我坚持用指数软化的原因之一。

5.2 冲击力曲线和能量曲线的诊断价值

冲击力-时间曲线是冲击仿真的核心输出。低速落锤冲击下，冲击力曲线一般有几个特征：初期有一个快速上升段，随后出现振荡平台，回程阶段冲击力下降。振荡的物理来源主要是层合板局部厚度方向的振动模态，因此有一定高频分量是正常的。如果振荡幅度异常大，优先检查沙漏能和接触刚度。

能量曲线里，内能（ALLIE）和动能（ALLKE）的转换是关注重点。冲击物初动能大部分转化成层合板的内能，一部分变成动能和摩擦耗散。如果总能量不守恒，检查接触是否漏了耗能项，或者沙漏能在持续增长。沙漏能ALLAE占比超过5%就要警惕，超过10%基本可以判定网格或沙漏控制参数有问题。在做参数研究的时候，这些曲线也能帮你快速筛掉那些明显不合理的工况组合。

5.3 与文献和实验数据的对比思路

子程序开发完不能只停留在“曲线看着合理”。有条件的话，拿文献里的冲击力峰值、分层面积、凹坑深度做定量对比。文献数据通常在低速落锤冲击的论文里很常见，注意匹配试件尺寸、铺层顺序、冲击能量和边界条件。

VUMAT相比内置模型的价值在这里就体现出来了：你可以把SDV输出映射成用户自定义的损伤指标，比如把基体拉伸损伤区域的投影面积和实验的C扫描分层面积对比，把纤维断层的投影区域和实验的X射线照片对比。只要材料参数和边界条件对得上，这类对比的说服力比单纯展示应力云图强得多。

6. 从Hashin到完整的冲击仿真模型：cohesive界面、Voronoi和效率优化

6.1 层内和层间失效必须分开处理

Hashin准则解决的是层内损伤——铺层内部的基体裂纹和纤维断裂。但复合材料冲击仿真里还有一个占比极重、能量耗散贡献很大的失效模式：层与层之间的分层。分层靠VUMAT是模拟不出来的，必须在层间建立cohesive单元或者定义cohesive contact。

最常见的做法是层间插入零厚度的cohesive单元，用双线性牵引-分离本构描述界面。界面参数主要是初始刚度、界面强度（法向和两个剪切方向）和断裂能。这里有个工程经验：初始刚度不是越大越好。初始刚度太大会导致cohesive单元自身的稳定时间步急剧下降，严重拖慢显式分析。通常取界面刚度为主材料模量的10到100倍即可，兼顾精度和时间步。

cohesive层对网格质量非常敏感，尤其厚度方向相邻单元的尺寸如果差异大，cohesive层会提前触发损伤或出现穿透。建议把cohesive层两侧的体单元尺寸控制得尽可能一致。

6.2 Voronoi建模的应用场景

Voronoi在ABAQUS复合材料仿真里的用途主要分两类：一种是在RVE（代表性体积单元）层面把材料细观结构离散成多边形晶粒或纤维束的胞元，用来做细观损伤分析；另一种是在宏观模型中沿层间界面或面内把cohesive层切分成Voronoi多边形碎片，模拟裂纹沿薄弱界面随机扩展的路径。

如果你的研究内容涉及层间裂纹的随机扩展路径，Voronoi界面是比连续cohesive层更贴近物理实际的选择。Voronoi在ABAQUS里通常靠Python脚本生成——先随机撒点，再计算Voronoi图，然后把多边形区域赋予不同的cohesive材料属性或切分成不同的单元集。这块工作量不小，而且网格质量是难点，尤其是多边形顶角处的过渡网格。我的建议是：如果不是研究课题明确需要，初期版本先别上Voronoi，把连续cohesive层跑通跑稳定再说。

6.3 GPU加速、并行和瑞利阻尼的取舍

显式冲击仿真计算量巨大，优化手段分成两类：硬件层面和模型层面。硬件层面，ABAQUS从2020版开始更完整地支持GPU并行计算。在Editor Environment中设置GPU设备后，显式分析在多GPU并行下提升明显，但效果和模型规模强相关。小模型GPU加速可能还不如CPU多核，二三十万单元以上的模型收益才比较显著。模型层面，域分解并行（Domain Decomposition）是Explicit分析默认的并行方式，你只需要在Job处理器里把并行核数调高。

关于瑞利阻尼，低速冲击仿真里很多人纠结要不要加。我的经验是：如果研究重点是冲击峰值力和损伤形貌，可以不加，加了反而容易掩盖真实的动力学响应。如果仿真和实验的冲击力曲线在卸载阶段总是对不上，或者低频残余振动明显，可以引入非常小的质量比例阻尼α来消耗低频振荡，但要对α的取值做敏感性分析。刚度比例阻尼β会显著影响稳定时间步，建议慎用，除非你有明确的频率标定数据。

6.4 断裂能引入才是从“能跑”到“可信”的关键

前面提到的指数退化是一种工程技巧，但它不是物理意义上的损伤演化方式。真正严谨的渐进损伤模型应该引入断裂能：Hashin准则只负责判断损伤起始，损伤起始后按照断裂能G和特征单元长度L，用线性或指数软化曲线计算损伤演化。这样网格尺寸敏感性会大幅降低，同一套材料参数在不同网格密度下能得到相对一致的冲击响应。

这算是VUMAT开发路线上一个比较长远的目标。起步阶段可以用简单的刚度折减把整个过程跑通，但当模型较准、做网格敏感性研究时，断裂能演化几乎是必须的。代码结构上只需要把损伤更新部分换成基于能量准则的演化解算，Hashin判定入口保持不变。整个框架不用推倒重来，这也正是自己从零写一份代码最大的优势——每一行都熟悉，改起来完全可控。

最后说一个老经验：VUMAT的调试过程一定要有耐心，出了问题先怀疑自己的代码和建模，再怀疑材料数据。大多数人写的Hashin代码第一版都有一些小问题，比如正负号、方向、状态变量序号错位。建议把每一版代码和对应的结果文件都留一个副本，哪怕方向错了或者参数错了的版本也别删。回头对照的时候，你会发现那些“错的版本”恰恰是帮你定位问题的最好线索。