考虑边界稳定的结构振动显式多时间步长计算方法

    马志强 楼云锋 金先龙

    

    

    

    摘要：多时间步长方法常用于求解有限元动力学空间多尺度问题，可分为单元重叠与非重叠两种形式。单元非重叠形式的多时间步长方法涉及边界数据的插值过程，造成计算不稳定、精度不高。为了提高多时间步长方法计算稳定性与求解精度，基于节点分割与重叠边界单元，提出一种改进的多时间步长计算方法。结合显式积分格式数据传递规律与重叠单元，分区边界数据由动力学方程显式求出，算法实施过程简单。采用能量方法分析了算法的稳定性能，给出了算法一般性稳定条件。稳定性分析结果表明改进方法的稳定性只由分区内部单元特性决定，重叠单元不影响算法的稳定性能，这扩展了改进方法的适用范围。数值算例结果表明该算法具有较高的计算精度，能量校验进一步证实了算法的稳定性能。

    关键词：多体系统动力学;显式积分;多时间步长;稳定性分析;能量校验

    中图分类号：0313.7;0242.2

    文献标志码：A

    文章编号：1004-4523 （2019）06-1041-09

    DOI：10. 16 385/j. cnki. issn. 1004-4523. 2019. 06. 013

    引言

    多体系统动力学以及车辆碰撞有限元仿真等问题中常涉及不同时间与空间尺度的耦合计算。这类问题需要根据结构设计特点、载荷分布以及分析区域划分不同属性的网格。采用中心差分法分析时整个模型采用单一时间步长。受限于稳定性要求，临界时间步长受模型最小单元属性的限制，部分尺寸较小的网格会使得整个仿真模型采用较小的时间步长，这会导致计算效率大幅下降。

    多时间步长方法或者子循环方法是处理这类问题的常用方法。所谓多时间步长方法是指根据需要将模型划分成若干区域，不同区域根据分区单元特性选择时间步长。一个多时间步长计算过程包含一个系统时间步与多个子循环时间步。

    Liu和Belytschko首先采用网格分割将多时间步长方法应用于瞬态动力学问题求解[1]。而后众多学者进一步研究改进方法，并将其运用到特定工程中。根据模型单元分割特点可以分为重叠与非重叠单元形式，不同多时间步长方法的区别在于分区边界数据的处理方法。

    非重叠单元形式的多时间步长方法研究中，以Smolinski和Daniel为代表的学者，在系统时间步内，以位移、速度或者加速度满足特定插值规律来处理各子区域边界节点间的相互作用[2-4]。Gra-vouil和Combescure以FETI（the finite elementtearing and interconnecting method）方法为基础，约束边界节点速度在子循环时间步内连续，提出了显隐混合多时间步长的GC方法[5]。Prakash和Hj elmstad改进了GC方法，通过拉格朗日乘子约束边界节点速度在系统时间步连续，提高了求解效率[6]。Lindsay等使用多时间步长方法处理诸如裂纹之类的结构非连续问题[7]。Ruparel等结合域分解方法与多时间步长方法处理非协调边界网格有限元模型，通过能量分析给出算法稳定的分区边界数据连续性条件[8]。

    国内学者中，陈丽华等使用区域分解方法，构造了冲击动力问题的混合时间步长显式积分并行算法[9]。高晖等比较了几种常用多时间步长方法的稳定性与精度，提出了一种应用于汽车碰撞过程仿真的基于常速度的阻尼子循环法[10]。缪建成等基于中心差分方法提出一种适用于柔性多体系统仿真的显式子循环方法[11]。

    由上可知，以非重叠边界为基础的多时间步长方法，处理不同步长分区边界数据传递，涉及边界变量插值，这会恶化多时间步长方法的稳定性能，限制分区步长比的选择[12]。现有的隐式或者显隐子循环算法都无法保证无条件稳定，边界数据插值会降低算法的稳定性，需要严格限制分区边界数据的连续性要求，对于显式子循环方法而言，还要限制各自分区内部临界时间步长。

    同时，这类方法存在着变量不一致性误差，进一步降低计算精度[13]。需要指出的是采用重叠单元形式的多时间步长方法的研究中，Arlequin方法是典型的方法[14]。Ghanem等将Arlequin方法用于结构瞬态动力学显隐混合多尺度分析[15]，而后Fernier等研究显式积分格式的Arlequin模型[16]。Arlequin方法使用拉格朗日乘子约束边界连续性，拉格朗日乘子可以理解为不同分区的边界耦合力。处理时间多尺度问题时Arlequin方法引入特定插值形式的拉格朗日乘子，这同非重叠形式的方法一样需要严格限制边界数据连续性条件。Rama等将刚度、质量与阻尼矩阵拆分为块对角矩阵与耦合矩阵，耦合的分区内力通过重叠单元获得[17]，但该方法只适用于分区同时间步长的计算。

    为了改善现有多时间步长方法对边界数据连续性的要求，提高不同分区步长比的选择范围，本文提出一种基于边界单元重叠形式的显式多时间步长计算方法。有限元仿真模型采用节点分割办法被划分为若干分区，分区边界重叠多层单元。各分区采用预测校正形式的显式Newmark積分格式。重叠边界单元方法符合显式动力学求解过程数据传递规律，重叠单元其实质为相邻分区预测波形的传递区域，因而处理分区多时间步长时不涉及边界数据插值过程，其实质是在子循环过程中，以重叠边界的变化代替边界数据的插值过程。

    由于不同分区采用不同的时间步长，一般形式的多时间步长方法的稳定性研究依旧是个热点课题。本文基于能量法分析算法的稳定性，给出了一般形式的稳定性条件。最后以典型数值算例验证方法的有效性。

    1 显式多时间步长计算方法

    条件稳定的显式求解格式，其临界时间步长受CFL条件（Courant-Friedrichs-Lewy condition）限制。为了保证显式计算的数值稳定性，采用的时间步长要小于临时时间步长。所有分区最小时间步长为△tr，由分区单元属性决定

    特别地，当γE一γL一1/2时，满足式（18）的稳定性条件即含有阻尼形式的中心差分法的稳定性条件。中心差分格式也适用于本文的多重叠网格形式的多时间步长计算方法。

    采用能量增量形式计算的重叠网格多时间步方法。重叠分区矩阵BEB同样满足式（17）的分析结论，除一般显式算法临界步长限制外，无需额外约束重叠分区内节点变量的连续性要求，就可以满足稳定性条件。

    3.2 能量校验方法

    多时间步长方法中，数值计算过程中的不稳定现象可以直观的由能量校核方法体现。对线性弹性结构有限元分析模型中，模型的内能W int、动能Wkin以及外力做的功wext可以表示为式中 K为刚度矩阵;M为质量矩阵;u，v为节点位移与速度矢量;fext为等效节点外力矢量。若计算过程数值稳定，则在系统时间步n内，能量校验满足式中 ε为一很小的误差容许极限（10-2）。通过检测任意时间步有限元模型能量的变化，可进一步验证稳定性分析的结论。

    4 数值算例

    采用经典的实体梁承受集中载荷模型为数值算例，梁两端简支。实体梁结构模型长L为20 m，高2m，厚1m。采用六面体网格离散，为了分析不同时间尺度对计算结果的影响，模型劃分了3种不同尺寸的有限元网格，离散后的实体梁模型有13720个单元，17370个节点，52110个自由度。实体梁模型有限元网格截面如图2所示。

    实体梁赋予各向同性材料，密度p=2200 kg/m3，弹性模量E=3×1010 Pa，泊松比v—0.22。梁沿着厚度方向划分5层单元，点A所在梁截面上所有节点施加等效集中阶跃载荷fext，载荷总大小为4×l06 N。仿真时间为0.4 s，采用显式多时间步长方法计算，不同网格尺寸施加不同的时间积分步长，计算梁中间位置点B的横向位移。模型的临界时间步长受分区3单元尺寸的约束，临界步长为1. 37×10-5 s。

    不同分区采用的时间步长如表1所示。其中m为分区1与3步长比

    采用步长比条件2即m=5计算梁中间点横向位移，参考方法采用经典的多时间步长方法[2]，采用模态叠加法计算梁中点响应的理论解。参考方法将显式与隐式积分格式统一，在子循环过程中，小步长所需边界节点位移实质为相邻大步长节点位移的线性插值。参考方法中求解格式采用预测校正New-mark格式，对大步长时间步边界节点位移校正值做线性插值传递到小步长分区边界节点。

    图3与4为本文所述方法与参考方法比较的结果，图4为计算结果的局部放大图。图3与4可以看出经典的多时间步长方法与本文方法在一定条件下都可以计算出梁横向位移。从图4中可以直观地看出，在相同步长比的条件下，多重叠形式的网格允许动力学预测波形在整个分区内传递，不涉及插值与波形截断，参考方法中子循环过程边界节点需要大步长分区节点的线性插值。相比于非重叠形式的参考解法，本文方法具有更高的计算精度。

    分别以2倍步长比、5倍步长比与1 0倍步长比计算梁中间节点的位移，梁中间位置横向位移计算结果局部放大图如图5所示。选取时间段0. 295-0.3 s，3条放大横向位移曲线表明，随着步长比的增加，梁中间位置横向位移曲线基本保持一致，计算结果表明本文方法具有较高的精度一致性。

    仿真时间为0. 297 s时，梁横向位移计算云图如图6所示。

    为了定量地分析不同算法位移计算结果，定义最大误差MaxEr与计算误差和SumEr两个计算指标。其中最大误差MaxEr=max‖Xm -Xref‖i，误差和SumEr=∑ Xm - Xref‖i。其中n为时

    间步数，Xref为某一变量参考值（理论值）。

    由于分区采用的时间步长不同，在相同仿真时间下，时间步长越小，计算时间步越多，这样就不具有相同的度量标准，因此全部采用步长比为m =2的条件下，大步长分区所需的时间步数作为标准n。统计梁中间点B横向位移的不同步长比下的方法精度如表2所示。

    从表2的计算结果可以看出，随着步长比的增加，本文方法最大误差与误差和均逐步增加。也就说明给定网格划分，采用步长比越大，计算精度相对降低。但是与经典的非重叠多时间步长相比，本文所述的方法在相同计算步长的前提下具有较高的计算精度，同时步长比越大，本文方法所述计算精度保持更高。

    采用本文3.2节能量计算方法计算模型在给定时间部内的外力功、内能、动能与检验能量，检验能量反映的是因为部分重叠网格产生的附加能量，这可以从侧面检验稳定性结论。图7为10倍步长比条件下，实体梁模型能量分析结果。

    在模型计算时间步内，采用显式多时间步长计算方法，当各个分区步长均满足CFL条件时，分区产生检验能量相对于模型内能、动能而言，数量级极小（10-4），未出现数值不稳定现象。

    为验证步长比对计算效率的影响，统计实体梁在不同步长比下所需计算时间，计算资源采用普通微机CPU i7-4790K。构造统一的计算时间统计标准，此时分区1的时间步长为6. 25×10-5 s，分区3为1. 25×10-6 s，分区2步长变动。分区2与3的步长比为m。仿真时间为0.4 s。

    由表3统计的模型计算时间可知：随着步长比m的增加，分区2的时间步长逐渐由1. 25×10-6 s增加到2.5×10-5 s，分区2的计算时间步数相应减少，计算时间逐步减少，也应看到步长比增加到5以后，计算时间没有明显继续下降，出现了负载不均衡现象。因而涉及多时间步长计算，负载均衡是要重点考虑的问题。

    为了验证时间步长效应对计算精度的影响，将图2中分区2与分区3统一划分较大尺寸网格，如图8所示。梁的几何尺寸与单元所采用材料属性均保持不变，梁沿着厚度方向划分3层网格。悬臂梁自由端所有节点施加正弦等效集中载荷F，载荷总大小为4×l06 N。载荷周期为0. 05 s，仿真时间为0.2 s。

    整个模型含有1020个单元，1588个节点，显式积分临界时间步长受分区2网格尺寸限制，为4. 54×10-5 s，选取分区2的时间步长为4×10-5s，受限于分区2的临界时间步长，此时步长比m最大可选为2。梁右端点上部节点C位移分别采用模态叠加法、本文所述方法以及显式MGMT方法计算[8]。同步长与2倍步长条件下位移计算结果如图9所示，参考的显式MGMT方法采用2倍步长比，小分区步长4×10-5 s。局部放大如图1 0所示。

    由图9中可知悬臂梁自由端受正弦载荷时，采用不同方法得到的自由端节点C位移得到大致相同的计算规律。2倍步长比MGMT求解不同分区边界通过Mortar方法耦合节点外力，同FETI方法相似，涉及插值过程。由图1 0可以看出，本文所述重叠单元的方法相对与固定边界的多时间步长方法，求解精度有所提高，同时提高步长比求解精度会降低。

    5 结 论

    多时间步长方法是处理局部网格多尺度问题的经典方法，采用分区网格边界重叠的方法，本文提出一种改进的显式多时间步长求解方法。通过能量法分析了算法的稳定性条件，给出了该方法一般性的稳定性条件。该方法的优点如下：

    （1）在子循环时间步内，预测波形得以完整的传递，不涉及插值、截断等，提高计算精度;可以进一步扩大分区步长比使用范围，缩短有限元空间多尺度问题的求解时间。

    （2）重叠多层网格不影响分区稳定性条件。这也意味着显式子循环过程只要各自分区满足临界步长稳定性要求，无需为重叠部分单元额外施加更加严格的稳定性条件。这提高了多时间步长方法的适用范围，同时方法易于编程实现。

    （3）从边界数据的传递规律来看，预测校正方法作为显式积分方法的特例，子循环过程是以重叠边界的变化代替边界数据的插值过程。本文方法可以推广到任意显式时间积分格式（如中心差分法）。

    由于显式积分方法的解耦特性，该方法容易并行化求解大规模以及超大规模结构动力学问题。这将在提高精度的基础上，进一步提升计算效率，这也是下一步研究的重点。

    参考文献：

    [1] Liu W K， Belytschko T. Mixed-time implicit-explicitfinite elements for transient analysis[J]. Computers&Structures， 1982，15 （4）：445-450.

    [2] Smolinski P，Belytschko T，Neal M. Multi-time-stepintegration using nodal partitioning[J]. InternationalJournal for Numerical Methods in Engineering， 1988，26（2）：349-359.

    [3] Daniel W J T. Analysis and implementation of a newconstant acceleration subcycling algorithm[J]. Inter-national Journal for Numerical Methods in Engineer-ing， 1997，40（15）：2841-2855.

    [4]Daniel W J T. A partial velocity approach to subcy-cling structural dynamics[J]. Computer Methods inApplied Mechanics and Engineering， 2003， 192 （3）：375-394.

    [5] Gravouil A， Combescure A. Multi-time-step and two-scale domain decomposition method for nonlinearstructural dynamics[J]. International Journal for Nu-merical Methods in Engineering， 2003，58 （10）：1545-1569.

    [6] Prakash A， Hjelmstad K D. A FETI-based multi-time-step coupling method for Newmark schemes instructural dynamics[J]. International Journal for Nu-merical Methods in Engineering， 2004， 61 （13）： 2183-2204.

    [7] Lindsay P，Parks M L，Prakash A. Enabling fast，stable and accurate peridynamic computations usingmulti-time-step integration[J]. Computer Methods inApplied Mechanics and Engineering， 2016， 306： 382-405.

    [8] Ruparel T， Eskandarian A， Lee J D. Concurrentmulti-domain simulations in structural dynamics usingmultiple grid and multiple time-scale （MGMT） meth-od[J]. International Journal of Computational Meth-ods， 2018，15（4）：1850021.

    [9] 陳丽华，程建钢，姚振汉，冲击动力问题的混合积分并行算法及应用[J].工程力学，2003， 20（2）：15-20.

    Chen Lihua， Cheng Jiangang， Yao Zhenhan. Mixedtime integration parallel algorithm and its applicaitionto dynamic impact problem[J]. Engineering Mechan-ics， 2003，20（2）： 15-20.

    [10]高 暉，李光耀，钟志华，等，汽车碰撞计算机仿真中的子循环法分析[J].机械工程学报，2005， 41（11）：9 8-101.

    Gao Hui， Li Guangyao，Zhong Zhihua， et al. Analysisof subcycling algorithms for computer simulation ofcrashworthiness[J]. Chinese Journal of MechanicalEngineering，2005， 41（11）：98-101.

    [11]缪建成，朱 平，陈关龙，等，多柔体系统响应计算的子循环计算方法研究[J].力学学报，2008， 40（4）：511- 519.

    Miao Jiancheng， Zhu Ping， Chen Guanlong， et al.Study on sub-cycling algorithm for flexible multi-bodysystem[J]. Chinese Journal of Theoretical and AppliedMechanics， 2008， 40（4）：511-519.

    [12] Klisinski M. Inconsistency errors of constant velocitymulti-time step integration algorithms[J]. ComputerAssisted Mechanics and Engineering Sciences， 2001，8（1）：121-139.

    [13] Karimi S，Nakshatrala K B.On multi-time-step mon-olithic coupling algorithms for elastodynamics[J].Journal of Computational Physics， 2014， 273： 671-705.

    [14] Dhia H B， Rateau G.The Arlequin method as a flexi-ble engineering design tool[J]. International Journalfor Numerical Methods in Engineering， 2005， 62（11）：1442-1462.

    [15] Ghanem A， Torkhani M， Mahjoubi N， et al.Arlequinframework for multi-model， multi-time scale and het-erogeneous time integrators for structural transient dy-namics[J]. Computer Methods in Applied Mechaniesand Engineering， 2013， 254：292-308.

    [16] Fernier A， Faucher V， Jamond 0. Multi-model Arle-uin method for transient structural dynamies with ex-plicit time integration[J]. International Journal for Numerical Methods in Engineering， 2017， 112 （9）：119 4-1215.

    [17] Rao A R M， Rao T V S R A， Dattaguru B. A newparallel overlapped domain decomposition method fornonlinear dynamic finite element analysis[J]. Comput-ers& Structures， 2003，81（26-27）：2441-2454.

    [18] Hughes T J R. The Finite Element Method： LinearStatic and Dynamic Finite Element Analysis [M].Courier Corporation， 2012.