Dyna动力分析实战指南从基础入门到高级应用解决工程难题与优化策略

引言：LS-DYNA在工程动力分析中的核心地位

LS-DYNA作为全球领先的显式动力学分析软件，在汽车碰撞安全、航空航天冲击防护、电子设备跌落测试以及金属成型制造等领域发挥着不可替代的作用。它能够模拟结构在极短时间内的大变形、材料失效、接触摩擦等复杂非线性行为。本指南旨在为工程师提供从基础理论到高级应用的完整学习路径，重点解决实际工程中的难题并探讨优化策略。

第一部分：基础入门——掌握LS-DYNA的核心概念

1.1 显式动力学与隐式动力学的区别

在进行动力分析之前，必须理解LS-DYNA所采用的显式积分算法（中心差分法）与传统的隐式积分法（如Newmark法）的根本区别。

显式动力学（Explicit）：
- 特点：无需迭代求解，利用当前时间步的已知量直接计算下一时间步的未知量。
- 优势：非常适合求解短时间、高动态、包含复杂接触和材料失效的问题。计算成本低。
- 劣势：时间步长受制于Courant-Friedrichs-Lewy (CFL) 条件，必须非常小以保证数值稳定性。
隐式动力学（Implicit）：
- 特点：需要在每个时间步进行迭代求解，形成并求解刚度矩阵。
- 优势：时间步长可以很大，适合求解低频、静态或准静态问题。
- 劣势：对于高度非线性问题，迭代可能不收敛，且单步计算成本高。

1.2 关键输入卡片详解：关键字文件（Keyword File）

LS-DYNA的输入文件是基于关键字的文本文件（通常以.k或.key为后缀）。理解核心关键字是建模的基础。

*CONTROL 系列：控制求解参数。
- *CONTROL_CONTACT：接触控制。
- *CONTROL_SOLUTION：求解器设置。
- *CONTROL_TIMESTEP：时间步长控制。
*PART：定义部件（Part），将材料、单元属性、壳/体厚度关联起来。
*SECTION：定义单元类型（如壳单元SHELL、体单元SOLID）及其属性（如积分点数、壳单元的算法）。
*MAT：定义材料模型。LS-DYNA拥有庞大的材料库，如*MAT_RIGID（刚体）、*MAT_PLASTIC_KINEMATIC（弹塑性）等。
*EOS：状态方程，主要用于流体、气体及高应变率下的金属材料。
*PART：定义接触界面。

1.3 第一个案例：简单的钢球撞击刚性墙

为了直观理解，我们构建一个简单的物理场景：一个钢球以5m/s的速度撞击刚性墙。

1.3.1 几何建模与网格划分（前处理）

通常使用HyperMesh、ANSA或LS-PrePost进行前处理。

创建一个半径为10mm的球体。
将其离散化为四面体或六面体单元（推荐使用四面体单元SOLID164）。
创建一个平面作为刚性墙。

1.3.2 关键字文件代码示例

以下是一个简化的K文件片段，展示了核心定义：

*KEYWORD
*TITLE
Ball Impact Example
*CONTROL_CONTACT
$      sidc       sid       mco        oic       ocs       mst        msf     sldthk    sldvo
         0         0         0         0         0         0         0         0         0
*CONTROL_SOLUTION
$      soln       nlq        ilc        its       cfl        dt0       etol      rtol
         1         0         1         0         0         0         0         0
*CONTROL_TIMESTEP
$      dtinit      tssfac      isdo       slso       dtmin      dtmax      tscel
      0.0000      0.9000         0         0      0.0000      0.0000         0
*CONTROL_ENERGY
$      hgen        rwen       slnten       ryco
         2         2         2         2
*DATABASE_BINARY_D3PLOT
$       dt        lcdt       beam       nplt       pset
      0.0100         0         0         0         0
*DATABASE_RCFORC
$       dt        lcdt       beam       nplt       pset
      0.0010         0         0         0         0
*PART
$      pid       secid       mid       eosid      hgid      grav      adpopt      tmid
         1         1         1         0         0         0         0         0
*SECTION_SOLID
$      secid       elform       aet
         1         1         1
*MAT_PLASTIC_KINEMATIC
$      mid        ro         e         pr         sy         fs         vp
         1   7.83e-9   2.06e+5      0.30   4.00e+2      0.00      0.0
$       c           p
         0         0
*PART
$      pid       secid       mid       eosid      hgid      grav      adpopt      tmid
         2         2         2         0         0         0         0         0
*SECTION_SHELL
$      secid       elform       shrf       nip       propt       qr/irid       icomp
         2         2      1.0000         0         0         0         0
$       t1        t2        t3        t4
      1.0000      1.0000      1.0000      1.0000
*MAT_RIGID
$      mid        ro         e         pr         n        couple         m
         2   7.83e-9   2.06e+5      0.30         0         0         0
$       c           v1        v2        v3
         0         0         0         0
*SET_PART_LIST
$      sid       da1       da2       da3       da4
         1         2         0         0         0
*CONTACT_AUTOMATIC_SINGLE_SURFACE
$      ssid        msid        sst        mst        sbox        mbox       spr       mpr
         1         0         0         0         0         0         0         0
*NODE
$      nid         x         y         z
         1   0.00000   0.00000   0.00000
...
(此处省略大量节点和单元数据)
*ELEMENT_SOLID
$      eid       n1       n2       n3       n4       n5       n6       n7       n8
         1         1         2         3         4         5         6         7         8
...
*BOUNDARY_PRESCRIBED_MOTION_SET
$      sid        dof        vad       lcid        sf        vid       death       birth
         1         1         0         1      1.000         0      1.0000      0.0000
*DEFINE_CURVE
$      lcid        sfa        sfo        offa        offo        dtp        tmin        tmax
         1      1.000      1.000      0.000      0.000         0      0.000      1.000
$         a               b
          0.0000000000000000E+00   0.0000000000000000E+00
          0.0010000000000000E+00   5.0000000000000000E+00
          1.0000000000000000E+00   5.0000000000000000E+00
*END

代码解析：

*MAT_PLASTIC_KINEMATIC：定义了钢材的弹塑性属性，密度7.83e-9 tonne/mm³，弹性模量206000 MPa，泊松比0.3，屈服应力400 MPa。
*CONTACT_AUTOMATIC_SINGLE_SURFACE：定义了单面接触，用于防止钢球自身穿透或穿透刚性墙。
*BOUNDARY_PRESCRIBED_MOTION_SET：给刚性墙（Part ID 1）施加约束，这里实际上应该给钢球施加初始速度。在LS-DYNA中，通常使用*INITIAL_VELOCITY_GENERATION来赋予初始速度，但为了演示约束，本例展示了如何通过曲线控制运动。更常见的做法是：
```
*INITIAL_VELOCITY_GENERATION
$      sid        dof        vad       lcid        sf        vid       birth       death
       1         3         0         0      1.000         0      0.0000      1.0000
$       vx        vy        vz        vx1       vy1       vz1
       0         0      5000.0         0         0         0
```
这里将Part 1（钢球）的Z方向速度设为5000 mm/s (即5 m/s)。

第二部分：工程难题解决——接触与材料非线性

2.1 接触算法的深度解析

接触问题是LS-DYNA中最容易出错也是最核心的部分。

2.1.1 接触类型的选择

*CONTACT_SINGLE_SURFACE：自动搜索接触，不区分主面（Master）和从面（Slave）。适合复杂装配体，但计算量大。
*CONTACT_AUTOMATIC_SURFACE_TO_SURFACE：需要定义接触对，算法更稳健，适合已知接触面的情况。
*CONTACT_ERODING：当单元失效被删除后，暴露出来的内部面可以继续参与接触。这在模拟侵彻、穿甲时至关重要。

2.1.2 接触刚度与穿透问题

工程难题：模拟中经常出现明显的穿透现象，或者由于接触刚度过大导致沙漏模式（Hourglassing）激增。

解决策略：

调整接触参数：
- SOFT：控制接触刚度的计算方式。设为1或2通常比默认的0更柔和，减少震荡。
- SLSFAC：从面刚度比例因子。如果出现穿透，适当减小该值（如0.1或0.01）。
网格质量：主面和从面的网格密度差异不应过大（通常建议比例小于4:1）。
壳单元厚度考虑：在*CONTROL_CONTACT中开启SHLTHK，确保壳单元的厚度被计入接触距离。

2.2 材料失效与损伤模型

在解决汽车碰撞或跌落测试时，材料往往会发生断裂。

2.2.1 常用失效模型

塑性应变失效：在*MAT卡片中定义FS（失效塑性应变）。当等效塑性应变达到该值时，单元被删除。
基于损伤的失效：如*MAT_ADD_EROSION，可以定义更复杂的失效准则，包括主应力失效、剪切失效、时间步长失效等。

2.2.2 案例：手机跌落屏幕破裂模拟

场景：手机以1.5m/s跌落，屏幕（玻璃材质）需要模拟破裂。

材料定义策略：玻璃通常表现为脆性材料。我们可以使用*MAT_BRITTLE_DAMAGE（材料号96）或*MAT_JOHNSON_HOLMQUIST_CONCRETE（材料号111，常用于陶瓷和玻璃）。

代码示例（玻璃材料片段）：

*MAT_JOHNSON_HOLMQUIST_CONCRETE
$      mid        ro         g        a        b        c        f
         3   2.50e-9   30000.0   0.009   0.007   0.006   0.008
$       n        fc        epso        es        mu        tl        s
       0.07   0.0000000000000000E+00   1.0000000000000000E-06   0.0000000000000000E+00   0.0000000000000000E+00   0.0000000000000000E+00   0.0000000000000000E+00
$       d1        d2        lcf       epsf       sflg       mflg       lcfm
       0.00       0.00         0       0.00         0         0         0
$       afit      bfit      cfit      dfit      efit      ffit      gfit
       0.00       0.00       0.00       0.00       0.00       0.00       0.00

注：实际应用中，玻璃材料参数需根据实验标定。

优化策略：

网格细化：在屏幕区域使用精细网格以捕捉裂纹路径。
沙漏控制：脆性材料删除单元后容易产生沙漏，需在*CONTROL_HOURGLASS中选择合适的算法（如Type 4或Type 5）并调整系数。

第三部分：高级应用——热-力耦合与优化

3.1 热-力耦合分析（Thermal-Mechanical Coupling）

在制动盘制动、焊接或高温成型中，温度场和应力场是相互影响的。

3.1.1 分类

单向耦合：先计算温度场，将温度作为预加载施加到结构上。
双向耦合：结构变形产生热量（塑性功、摩擦热），热量改变材料属性（软化），同时影响结构响应。

3.1.2 实现方法

在LS-DYNA中，需要使用*BOUNDARY_PRESCRIBED_TEMPERATURE定义边界温度，并在材料模型中定义随温度变化的属性（*MAT_..._T系列）。

摩擦生热代码示例：在接触定义中，开启摩擦生热开关：

*CONTACT_AUTOMATIC_SURFACE_TO_SURFACE
$      ssid        msid        sst        mst        sbox        mbox       spr       mpr
         1         2         0         0         0         0         0         0
$       fs        fd        dc        vc       vdc    penchk        bt        dt
      0.15      0.15       0.0       0.0       0.0         0       0.0       0.0
$       sfs       sfm       fsf       vsf
       1.0       1.0       1.0       1.0

同时需要定义摩擦生热参数（通常在*CONTROL_CONTACT中设置SLSFAC和RWPNAL等，或者在接触卡片中指定摩擦系数随温度变化的曲线）。

3.2 DOE（实验设计）与优化策略

面对复杂的工程问题，如何找到最优解？LS-DYNA通常与优化软件（如LS-OPT）结合使用。

3.2.1 响应面法（Response Surface Methodology, RSM）

通过在设计空间内选取有限的样本点，拟合出一个近似的数学模型（响应面），然后在该模型上寻找最优解。

流程：

确定设计变量（Design Variables, DV）：如梁的厚度、材料的屈服强度、泡沫的密度。
确定目标函数（Objective Function, OBJ）：如最小化加速度峰值、最大化吸能。
确定约束（Constraints, CON）：如质量不超过某值、侵入量不超过某值。

3.2.2 优化案例：汽车保险杠吸能盒优化

目标：在满足碰撞法规（如低速碰撞）的前提下，最小化重量。

步骤：

定义DOE实验：选取吸能盒的厚度（T1, T2）和诱导槽深度（D）作为变量。

生成K文件模板：使用LS-PrePost或脚本将变量参数化（使用*PARAMETER）。


*PARAMETER
$ name      value
T_thick     1.5
D_slot      2.0
*PART
$      pid       secid       mid       eosid      hgid      grav      adpopt      tmid
     1         1         1         0         0         0         0         0
*SECTION_SHELL
$      secid       elform       shrf       nip       propt       qr/irid       icomp
     1         2      1.0000         5         0         0         0
$       t1        t2        t3        t4
{T_thick} {T_thick} {T_thick} {T_thick}

运行LS-OPT：配置采样方法（如拉丁超立方），提交计算。
结果分析：查看响应面，找到满足约束下的最小重量组合。

第四部分：常见错误排查与性能优化

4.1 沙漏模式（Hourglassing）的控制

现象：网格出现零能量变形（像手风琴一样折叠），导致计算结果不准确甚至发散。

原因：单点积分单元（为了计算速度）缺乏刚度，无法抵抗这种变形模式。

解决方案：

在*CONTROL_HOURGLASS中选择合适的Hourglass类型（Type 4或Type 5最常用）。
设置Hourglass系数（通常在0.05到0.1之间，不要超过0.1，否则会过度硬化）。
避免在关键区域使用单点积分的壳单元，改用全积分单元（ELFORM=1或2），但这会显著增加计算时间。

4.2 负体积（Negative Volume）与时间步长过小

现象：计算在几步后崩溃，提示Negative Volume或时间步长为零。

原因：

单元畸变严重（通常由大变形或接触刚度不匹配引起）。
材料模型参数错误导致刚度为负。

解决方案：

检查接触定义，确保接触刚度设置合理。
使用*CONTROL_SHELL中的WARP参数来增强壳单元的抗翘曲能力。
如果是体单元，考虑使用*CONTROL_SOLID中的相关设置，或更换材料模型算法。

4.3 计算加速技巧

质量缩放（Mass Scaling）：为了增加最小单元的临界时间步长，可以人为增加单元的质量（密度）。
- 原则：增加的总质量应小于模型总质量的5%，且不影响整体惯性。
- 实现：在*CONTROL_TIMESTEP中设置DT2MS（目标时间步长）和MS1ST（仅对第一步生效）。
并行计算：利用LS-DYNA的MPP（Massively Parallel Processing）版本，在多核服务器上运行。通常核心数在16-32核时效率较高，超过此数需考虑通信开销。

结语

LS-DYNA是一个强大的工具，但其威力取决于使用者对物理问题的理解和对软件参数的把控。从简单的静态分析到复杂的流固耦合，核心在于“建模-求解-验证”的循环。希望本指南提供的代码实例和策略能帮助您在工程实践中攻克动力分析的难题，实现产品性能的优化。记住，仿真不仅仅是画网格，更是对物理世界的数学重构。