结构有限元分析中的时域瞬态法

瞬态结构动力学分析（又称时间历程分析）是用于确定结构随时间变化的动力学响应的一种方法。通过瞬态动力学分析，我们可以确定结构在稳态载荷、瞬态载荷、和简谐载荷的随意组合作用下的随时间变化的位移、应变、应力及力。瞬态分析最大的优点就是考虑了惯性力和阻尼作用。

瞬态问题常见的分析方法有：

• 时域瞬态分析

• 特征值提取(自然频率和模态)

• 稳态响应(频域的谐响应分析)

• 响应谱分析(冲击的峰值响应计算)

• 随机响应分析(由随机激励引起的振动)

由于时域瞬态方法最为直观，适应各种线性、非线性的工况，使用率高。本文只介绍时域瞬态法。其他方法将会在以后的文章中介绍。

控制方程

瞬态动力学的基本运动方程是：

其中：[M]是质量矩阵。[C]是阻尼矩阵。[K]是刚度矩阵。{u_tt}是节点加速度。{u_t}是节点速度。{u}是节点位移。{F}是载荷。可以看出结构瞬态方程是一个含有二阶时间导数的方程。求解二阶时间导数的方法很多，结构有限元中使用最为广泛的是Newmark隐式时间积分法。WELSIM默认的结构时间求解器也是Newmark法。

常见求解方法

瞬态动力学有常见三种有限元求解方法：完全法、缩减法及模态叠加法。 

完全法采用完整的系统矩阵计算瞬态响应。它是功能最强的，支持求解各类非线性特性（塑性、大变形、大应变等）。优点是：容易使用，不必关心选择主自由度或振型。允许各种类型的非线性特性。采用完整矩阵，不涉及质量矩阵近似。一次分析就能得到所有时间历程下的位移和应力。允许所有类型的边界条件。完全法的主要缺点是开销大，费时，计算所得数据大。

缩减法通过采用主自由度及缩减矩阵压缩数据规模。在主自由度处的位移被计算出来后，将解扩展到完整自由度集上。优点是：比完全法快且开销小。缺点是：初始解只计算主自由度的位移。只能施加节点边界条件。所有载荷必须加在用户定义的主自由度上。时间步长必须恒定，不支持自动时间步长。不支持非线性（点对点接触除外）。

模态叠加法通过对模态分析得到的振型（特征值）乘上因子并求和来计算结构的响应。优点是：比缩减法或完全法更快、开销更小。允许考虑模态阻尼（阻尼比作为振型号的函数）。缺点是：时间步长必须恒定，不支持自动时间步长。不支持非线性（点对点接触除外）。不能施加非零位移边界条件。

由于完全法的优越性，及其在非线性问题的广泛应用。本文只介绍完全法。

阻尼（Damping）

考虑阻尼的作用是瞬态分析优点之一。阻尼可以看做是一种能量的耗散，它来自于诸多因素，而有限元分析中，时常将阻尼看做是一种综合作用，而从数值角度施加阻尼。 

有限元中常见的阻尼设置有三种：直接阻尼（direct damping），Rayleigh阻尼（Rayleigh damping）和复合阻尼（composite damping）。直接阻尼可以定义与每阶模态相关的临界阻尼比，其典型的取值范围是在临界阻尼的1%到10%之间。Rayleigh阻尼假设阻尼矩阵是质量和刚度矩阵的线性组合。尽管阻尼是正比于质量和刚度矩阵的假设没有严格的物理基础，实际上我们对于阻尼的分布知之甚少，而实践证明Rayleigh阻尼在有限元法中是有效的，并被广泛应用。复合阻尼可以根据每种材料定义一个临界阻尼比，这样就得到了对应于整体结构的复合阻尼值。当结构中有多种不同的材料时，复合阻尼更为有效。

在大多数线性动力学问题中，准确地定义阻尼对于结果十分重要。但是，阻尼算法和参数只是近似地模拟了结构吸收能量的特性，并非从原理上模拟引起这种效果的物理机制。因此，有限元分析中确定阻尼数据是很困难的。有时我们可以从试验中获得这些数据，有时必须要通过查阅参考资料或者经验来确定阻尼参数。

时间求解器

瞬态结构问题，由于二阶时间导数的引入，我们需要能求解时间导数的时间求解器。我们一般将此求解器分为两大类：显示求解器，和隐式求解器。

显式求解器用上一步的结果和当前步的结果计算下一步的计算结果。显示求解器有不稳定区域，需要很小的时间步长。好处是可以不需要非线性求解器（牛顿迭代求解器）来求解非线性特征，也不需要组装刚度矩阵。同时编码和算法相对简单，不需要额外的内存来存储额外中间数据。显示求解器的稳定步长可以估算得出，但是由于实际计算的复杂性，会设置的比理论步长还要小一些。显式算法要求质量矩阵为对角矩阵，而且只有在单元级计算尽可能少时速度优势才能发挥, 因而往往采用减缩积分方法，容易激发沙漏模式，影响应力和应变的计算精度。显示求解器的代表算法有中心差分法，Euler向前差分法，Runge-Kutta，线加速度法等。

隐式求解器用当前步结果和下一步未知结果反复迭代下一步结果，必须通过迭代得到。隐式求解需要组建刚度矩阵，需要牛顿迭代求解非线性问题。常用方法有Newmark法和Wilson-Theta法，其变种算法可以通过修改算法中的alpha、beta和theta参数来实现，如HHT算法在一些问题上的计算效率和精确度有显著提高。隐式求解法的最大优点是它具有无条件稳定性，即时间步长可以任意大。但是实际运算中上要受到迭代次数及非线性程度的限制，需要取一个合理值。同时如Newmark等算法具有二阶精度。

显示和隐式算法的比较如图所示

为了确保得到结构的全部反应并保证解的稳定性和收敛，选择正确的时间步长是很重要的。一般说来，时间步越短，解越精确。然而时间步越小，求解步数就越多，运行时间将显著增加。因此，求解所需要的时间将限制步长不可能太小。同时，时间步不能太大，否则计算容易漏掉结构的很多高阶频率，从而导致所得到的解不真实。

自动时间步长是一种能够优化瞬态求解效率的自动化算法，用户只需要设置初始步长和步长的最大最小值，在求解过程中，求解器将根据需要自动减小时间步，以解析解中的任何快速变化。在结果变化较小的过程中，求解器能增大时间步，从而提高计算效率。

边界条件和初始条件

和静力分析不一样的是，由于瞬态和惯性的特征。瞬态分析支持速度和加速度边界条件。这些边界条件，尤其是加速度可以用来模拟激振的效果。

对于初始条件，瞬态分析一般会支持一下三种：

• 线性静力结果。它指定了结构在某种静荷载作用下的初始位移。

• 一个瞬态分析结果，它确定了结构在某时刻的瞬态反应。求解可以从任何指定的时间步开始。

• 结构中所有自由节点的初始速度和初始加速度。

初始速度条件是比较常见的初始条件，它常以体条件施加在某个结构体的所有节点上。

以上是结构时域瞬态分析中的核心知识。下面我们来看如何在有限元软件中如何来实现计算和分析。

结构瞬态分析步骤

1. 建立或导入模型

2. 网格划分

3. 载荷步与时间步设置

Number Of Steps：载荷步，表示分几步施加载荷。

Current Step：当前载荷步。

Current End Time：当前载荷步结束时间。

Auto Time Stepping：自动时间步长是否打开。WELSIM目前对结构分析只支持固定步长，所以关闭此选项。

Define By：定义载荷子步的方式。可以通过时间和载荷子步数来定义，用时间定义的意思是每一载荷子步经历的时间是多少，用载荷子步定义的意思是一个载荷步有多少个载荷子步。

Time Step：载荷子步的时间。

4. 设置边界条件，初始条件，接触等。

5. 计算并验证结果

下图显示Y方向变形和von-mises应力，已经在整个时间历程下的最大最小值。

可以看出本分析中不含有阻尼，所以往复振动无衰减。关于阻尼的作用与设置，我们会在以后的文章中介绍。