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

来源:用户上传      作者:

　　摘要：多时间步长方法常用于求解有限元动力学空间多尺度问题，可分为单元重叠与非重叠两种形式。单元非重叠形式的多时间步长方法涉及边界数据的插值过程，造成计算不稳定、精度不高。为了提高多时间步长方法计算稳定性与求解精度，基于节点分割与重叠边界单元，提出一种改进的多时间步长计算方法。结合显式积分格式数据传递规律与重叠单元，分区边界数据由动力学方程显式求出，算法实施过程简单。采用能量方法分析了算法的稳定性能，给出了算法一般性稳定条件。稳定性分析结果表明改进方法的稳定性只由分区内部单元特性决定，重叠单元不影响算法的稳定性能，这扩展了改进方法的适用范围。数值算例结果表明该算法具有较高的计算精度，能量校验进一步证实了算法的稳定性能。
　　关键词：多体系统动力学;显式积分;多时间步长;稳定性分析;能量校验
　　中图分类号：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，由分区单元属性决定 　　△tr=a min<△tker） （k=l，2，…，nelem，） （1）式中△tkr为分区第k个单元属性决定的时间步长;a为比例系数，常用0. 9-0. 95;nelem，为有限元模型单元数目。
　　定义分区步长比mi为第i个分区采用的时间步长△ti与所有分区最小时间步长的比值
　　为了表示子循环时间步内的变量，以Xn+j/m代表某一变量，其中n为系统时间步，j为子循环时间步
　　采用基于预测校正形式的显式Newmark格式，求解步骤分为预测步与校正步。预测步中，任意分区内部节点的速度Vil，n+j/mi与位移uiI，n+j/mi的预测公式可以写为[17]
　　式中 β，γ为Newmark积分参数;下标i为分区编号，I代表分区内部节点。u，v与a分别表示节点位移、速度与加速度。时间离散的分区结构动力学方程可以改写为式中 质量矩阵m采用集中质量矩阵，为对角矩阵;k为刚度矩阵;c为阻尼矩阵，采用Rayleigh阻尼形式;f为分区节点外力向量;E为分区外部节点。分区有限元方程自由度按照先内部节点后外部节点的顺序编号。
　　分区内部节点的加速度求解公式可以表示为
　　边界节点加速度求解因未與分区通信，包含求解信息不完整，一个系统时间步后无法正确求出。基于重叠单元求解方法，边界节点为相邻分区的内部节点。边界节点加速度信息由相邻分区内部节点传递而来。
　　节点的位移与速度通过Newmark求解格式的校正步计算得到
　　公式（4）一（7）为多时间步长的显式求解格式，预测校正形式的显式多时间步计算流程如图1所示。图1中边界节点与外部节点构成重叠网格区域。B，与E2，E1i与B2i分别具有相同节点编号。
　　2 分区边界数据的更新顺序
　　非重叠形式的多时间步长方法，边界节点位移、速度与加速度均假定为时间的线性插值函数。事实上这存在数据的内在不一致性，因为假定加速度为线性变化，速度则为二次函数，位移为时间的三次函数。
　　而重叠网格形式的多时间步长方法其数据传递满足一致性要求。求解过程中，节点加速度是位移、速度校正过程的核心数据。为满足重叠部分位移、速度与加速度连续性要求。边界数据更新顺序如下列所示：
　　（1）给定分区初始节点位移、速度信息，计算节点初始加速度，给定分区采用的时间步长△ti与积分参数β，γ;
　　（2）小步长分区所有节点采用子循环时间步长△tr进行预测步计算，由式（6）计算节点加速度进而校正子循环步节点速度与位移预测值;
　　（3）大步长分区采用系统时间步长值，计算节点加速度，进而校正速度与位移预测值;
　　（4）判断小步长分区子循环步数是否达到分区步长比;若否小步长分区继续（2）;
　　（5）将各自分区的边界节点加速度、速度以及位移组成广义数据结构，传递到相邻分区外部节点，完成分区外部节点的信息更新。
　　3 显式多时间步长方法的稳定性分析
　　多时间步长法的稳定性校验多采用能量法与算子谱半径法两种。算子分析法多针对线性非重叠网格分析。本节采用能量分析方法分析网格重叠多时间步长算法的稳定性，给出有关该算法稳定条件的一般性结论。
　　3.1 基于能量法的稳定性分析
　　能量法分析稳定性中，数值积分算法的稳定性条件为：对于积分过程的任意时间步n，存在关于计算结果的范数Sn，满足条件：Sn≤kSo，其中k为正常数。能量分析方法多用增量能量形式，即时间步n+1的模型能量相对于时间步n时的系统能量增量。重叠网格的多时间步长方法，分区边界采用重合单元连接。在一个时间步内，模型的总能量包含大步长分区、小步长分区以及边界重叠单元能量三个部分。
　　为了便于分析，首先定义变量前向差分算子与平均算子式中 x代表某一变量，下标为时间步。采用能量增量形式的模型总能量式中 下标L表示采用系统步长△t分区，inter为重叠分区，E为采用子循环时间步长△t，分区。小步长分区采用预测校正格式显式Newmark格式，其在时间步j时能量增量为[18]△tr（7E一1/2）CE。a与y为节点加速度、速度矢量。ME，CE与KE分别为小步长分区质量、阻尼与刚度矩阵。γE为Newmark时间积分参数。则小步长分区在一个系统时间步内能量增量可以表示为
　　重叠分区采用子循环时间步求解，其能量增量形式可以采用小步长分区的表述式中与KEB分别为网格重叠区域质量、阻尼与刚度矩阵。
　　对于采用系统时间步长△t的分区，一个系统时间步内的能量增量可以表示为式中分别为大步长分区质量、阻尼与刚度矩阵。
　　在一个系统时间步内，模型中总能量增量的表示形式可以写成
　　将分区能量计算公式（10），（11）与（12）带人系统时间步（n+1）时刻的系统总能量计算公式（9），并简化可以得到
　　若数值方法保持稳定，那么需要计算过程中不产生附加伪能量，即对有限元模型需要满足
　　根据模型能量增量表达式（15）与稳定性要求（16），可进一步分析网格重叠多时间步长方法的稳定性。
　　结构动力学计算过程中阻尼矩阵常采用Ray-leigh阻尼形式，则在表达式（15）中涉及阻尼矩阵的后三项可以统一表示为（n）T C（n）≤O（不含阻尼形式包含在内）满足形式（16）。由式（15）知该多时间步长算法为条件稳定的算法，若算法稳定需要满足条件：γE≥1/2，γL≥1/2且矩阵B9B与BEB为非负定矩阵。因为B，B与BEB具有相同的表达形式B=M一1/2c一1/2γL△t2K，统一分析该矩阵非负定条件，不加上标区别。对于结构任意阶模态均需要满足式中 ξ为结构第n阶模态阻尼比，△t为时间步长，wn”为第n阶结构模态频率。对于有阻尼系统，本文所述方法的稳定性条件可以改写为式中 Tn为系统的最小固有周期。对于采用子循环时间步的分区步长应该满足mi Δt，≤△ter。 　　特别地，当γ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.
转载注明来源:https://www.xzbu.com/8/view-15126099.htm