指数时间差分法在旋转柔性梁动力学仿真中的应用策略①

方思明 王贤明 郭书豪

(浙江工业大学机械工程学院 杭州 310014)

0 引言

在诸如机器人机械臂、直升机旋翼、涡轮机叶片等的实际工程应用中,其力学模型往往涉及旋转柔性梁[1]。当使用有限单元法等数值方法离散其动力学方程后,其离散动力学方程通常会具有数值刚性,且刚性病态程度与空间离散化的网格大小相关[2]。对刚性动力学方程进行高效、精确、稳定的数值求解,是柔性体系统动力学研究的重要内容,也是动力学仿真中的难点之一[3]。

刚性动力学方程数值求解的关键,在于选择一个有效的微分方程时间积分算法。常用的动力学时间积分算法有四阶经典Runge-Kutta 法、Gear 法、Adams 法、Newmark 法,等等[4]。时间积分算法一般分为显式算法和隐式算法[5]。一般地,这两类算法求解刚性问题时都可能存在计算效率的问题,其中显式算法单步计算效率高。但是,当求解刚性问题时,为了避免数值发散,一般的显式算法必须使用非常小的时间步长才能稳定地求解,而问题本身的解总体变化缓慢,并不需要很高时间分辨率的小步长。这是数值刚性问题的典型特征之一。隐式算法比显式算法更稳定,时间步长一般只需要满足时间分辨率的要求,没有稳定性的额外要求;但隐式算法一般每一步都需要求解代数方程组,对于大型非线性动力学系统,计算代价较高甚至无法计算。因此,对于大型刚性系统,目前主要采用显式算法进行求解[6]。于是,人们期望发展出一种显式算法,在满足时间分辨率的基本要求下,能尽可能地放大时间步长以提高整体求解效率。其中这样一类显式算法称为指数积分法[7],指数积分法是利用了微分方程Jacobian 矩阵或其近似矩阵的矩阵指数和相关函数的一类数值格式,能够有效地处理一些数值刚性问题。通过矩阵指数运算,指数积分法避免了对引起方程刚性的部分直接进行数值积分,从而可以使用较大的步长进行稳定的数值求解。

指数时间差分法(exponential time differencing,ETD)[8]是一种典型的指数积分法。ETD 方法根据常数变易法,得到微分方程的一个积分形式解析解公式,通过多项式逼近积分中包含的“非齐次项”函数,再半解析地解出积分而得到数值解。Cox 和Matthews[8]给出了任意阶ETD 格式的推导以及Runge-Katta 类型的ETD 格式。指数积分法中,指数矩阵的计算是限制该方法的因素之一。随着矩阵指数计算方法[9]的发展,ETD 格式在各个领域得到广泛应用[10-13]。ETD 格式在运用过程中,其线性算子的选取会很大程度地影响其格式的形态[7,14],然而针对线性算子选取的相关研究目前并不是十分充分。

本文针对非惯性系下的旋转柔性梁这一典型动力学系统的数值仿真,研究采用ETD 求解时其中线性算子的选取策略。旋转柔性梁存在动力刚化现象,梁的刚度会随角速度变化。一次近似理论[15-16]建立的动力学模型可以反映动力刚化现象,此时动力学方程的刚度矩阵随系统大范围摆动的情况而变化。考虑到采用ETD 求解刚性问题时,线性算子选取对格式的稳定性会有重要的影响,因此,本文在应用ETD 格式求解旋转柔性梁问题中,考虑将系统大范围摆动的信息纳入线性算子的构造,并通过数值实验验证了方法的有效性。

具体地,本文的后续内容安排如下。第1 节介绍指数时间差分法(ETD)。第2 节给出了旋转柔性梁的动力学方程。第3 节给出了相应的线性算子的选取策略。第4 节是数值算例,用来验证上述的选取策略。最后是总结与讨论。

1 指数时间差分法

本节介绍指数时间差分法的基本思想以及具体的几个格式。本文将用该格式求解非惯性系下旋转柔性梁的动力学方程。考虑如下一般形式的微分方程:

其中,y=y(t)是维度为n×1 的未知变量。当采用指数积分法求解时,将式(1)右端项改写f(t,y)成如下形式:

其中,n×n的常系数矩阵A称为指数积分法的线性算子,文中简称为线性算子。A取定之后,可以确定拟非齐次项(t,y)=f(t,y)-Ay。在给定y(tk+h)=yk下,求得下一时间步tk+h时刻的解,就是求解微分方程的逐步法。根据常数变易法,可以得到式(2)如下积分形式的解析解公式:

式(3)是精确成立的,但是,关于下一时刻解y(tk+1)本身是隐式的。对于一般问题,这并不方便直接求解,需要将式中的积分项进行显式化计算。ETD 方法一般采用只含时间t的函数,在tk≤ζ≤tk+1内对(ζ,y(ζ)) 进行逼近,通过分部积分就可以解析地求得积分。例如:若在tk≤ζ≤tk+1上采用零阶多项式的常数逼近,即(ζ,y(ζ))=(tk,yk),再解析地积出积分项,就得到了指数欧拉法,用ETD1 表示[8]:

其中,I 为n×n单位矩阵。

指数时间差分法之所以能够有效处理刚性问题,关键在于通过解析地求得了包含快速变化指数项的积分,从而避免了对方程的刚性部分进行直接的数值积分。而一般显式格式数值求解刚性问题之所以不稳定,就是由于对于快变刚性部分的直接数值积分带来了数值误差的快速放大。对于“拟非齐次项”采用不同的函数逼近方式,研究者发展出了不同的ETD 方案。一方面,产生了众多的指数时间差分线性多步法[17-18];另一方面,进一步发展出了指数时间差分Runge-Kutta(Runge-Kutta ETD,ETDRK)方法等[19]。其中,二阶和四阶指数时间差分Runge-Kutta 格式分别如下。

(1) 指数Runge-Kutta 的二阶格式,用ETDRK2表示:

其中,r1、ak、φ1均为中间变量,公式推导见文献[8],具体形式如下。

(2) 指数Runge-Kutta 的四阶格式,用ETDRK4表示:

其中,r1、r2、r3、ak、bk、ck、φ1均为中间变量,公式推导见文献[8],具体形式如下所示。

另外,为保证ETD 格式求解刚性微分方程,需要选取合理的线性算子。通常线性算子的选取更多是经验性的,一般直接选取方程线性项的系数矩阵作为线性算子[7,14]。当线性项不是引起方程刚性的主要部分,甚至方程中并不存在显式的线性项时,这样选取可能会出现问题,影响指数积分法的精度和稳定性。

2 旋转柔性梁动力学方程

本节中,针对旋转柔性梁这一典型系统,采用一次近似理论[16]建立其动力学方程,并使用有限单元法进行空间离散。

考虑图1 所示的旋转柔性梁系统。平面柔性梁绕O点做定轴转动。假设柔性梁的材料分布均匀且各向同性,变形时梁的横截面仍然保持为平面且与中轴线垂直。于是,可以采用平面细长Euler-Bernoulli 梁模型。梁的静长度为L,横截面积为S,截面惯性矩为I,材料的质量体积密度为ρ,弹性模量为E。OX0Y0为惯性坐标系,固定不动,OXY为浮动坐标系,随柔性梁转动。位于旋转运动轴点的柔性梁端部的轴线与惯性坐标系横轴OX0的夹角为θ。对于柔性梁的变形位移描述如图2 所示。P点为梁上任意一点,P0点为P点未变形前的位置。

图1 旋转柔性梁

图2 柔性梁在浮动坐标系下的变形位移描述

柔性梁上的P点的变形可以用2 个变量来描述:纵向变形量w1(x,t) 和横向弯曲挠度w2(x,t)。那么,柔性梁上P点的变形位移可表示为

其中,wc1(x,t) 是横向弯曲引起的梁纵向变形量:

P点在惯性基下的矢径为

Θ浮动坐标系相对于固定坐标系的方向余弦阵:

P0点关于浮动坐标系OXY的坐标列向量r0=[x0]T。于是,系统的动能和势能表达式为

根据Hamilton 原理得到连续动力学方程,并使用有限单元法以n个单元进行空间离散,得到离散动力学方程[20]:

其中,v=[θ,q]T,q是维度3(n+1)×1 的有限元节点广义坐标向量。总体质量矩阵、总体陀螺效应矩阵、总体刚度矩阵分别为

右端QB向量和FB向量分别为惯性力阵和广义外力阵。

为了采用ETD 格式求解,首先将式(18)写成标准的一阶常微分方程组形式:

式中,系数矩阵中M-1K和2M-1G都与系统大范围摆动的信息相关。当角速度随时间变化时,动力学方程为时变常微分方程组。指数积分法求解时,需要仔细考虑线性算子的选取。

3 线性算子的选取

当使用ETD 求解时变或非线性问题时,为了保证求解稳定,需要合理地选取线性算子。针对系统大范围运动已知的旋转柔性梁问题,考虑到其动力学特征,根据系统大范围转动的角速度来构造ETD的线性算子。

在指数积分法的实际应用中,通常经验性地取方程线性项的系数矩阵作为格式中的线性算子。当方程的刚性主要来源于其线性项时,这种选取线性算子的方式是有效的。但是,当方程的刚性并不完全来自于显式线性项时,若直接取显式线性项作为线性算子,相应的拟非齐次项仍会有较大的刚性。而这将弱化ETD 格式求解刚性问题的优势,甚至会令其退化为常规的显式格式。因此,对于刚性并不完全简单地来自于显式线性项的问题,指数积分法中线性算子的选取,需要有进一步的考虑。

文献[21]研究了求解动力学问题时,ETD 格式中线性算子的选取对格式稳定性的影响,并从稳定性的角度给出了线性算子的选取建议。对于变系数或非线性动力学问题,在无法完全动态选取Jacobian 矩阵,所选线性算子偏差始终存在的一般情形下,选取对应刚度略偏大而非偏小的线性算子是相对安全的。在这样选取下,指数积分法底层方程解的变化速率大于所求解原始方程解的变化速率;相较于对应偏差相当而刚度偏小的线性算子选取方式,此时的指数积分法展现出了更好的稳定性。这里,所谓的底层方程是指以指数积分法的线性算子作为线性项系数矩阵的齐次微分方程[22]。

对于旋转柔性梁系统,大范围的旋转运动会影响柔性梁的变形位移响应[22]。由于离心力的存在,当转动的角速度变大,柔性梁的刚度会随之变大,出现动力刚化现象[23]。这在高速旋转中尤为明显。对于大范围运动已知的旋转柔性梁系统,根据前述以一次近似理论建立的动力学方程式(18),其总体刚度矩阵与大范围摆动的信息相关:

其中,Ki、Mi、D1、GT都是常数矩阵。GT矩阵元素的数量级很小,对于角加速度不太大的一般情况,其对系统总体动力特性的影响较小。总体刚度矩阵中的第2 项,与角速度的平方相关。当角速度达到最大值时,系统的各阶固有频率一般将会达到最大值。此时,方程各个解分量的变化率是最快的。于是,在应用指数积分法求解时,需要考虑到旋转梁系统的这一动力学特征。如果仅以方程的线性常系数项系数矩阵作为线性算子,不能完整地涵盖系统的快变“刚性”以及“振荡”信息,拟非齐次项仍会含有较大的刚性,这将会导致ETD 格式求解时出现数值发散。

考虑到旋转柔性梁的动力学特征,本文将系统大范围摆动的信息纳入指数积分法线性算子的构造。总体刚度矩阵式(25)相应的Jacobian 矩阵作为线性算子:

由式(25)可知,总体刚度矩阵与角速度和角加速度相关。考虑到高速旋转但角加速度较小的多数情况,如前所述,项对系统的固有特性影响很小,而角速度对柔性梁刚度和系统动力学特性的影响较大。因此,本文在线性算子的构造中只考虑角速度的变化。如果需要动态选取线性算子,这样依据单一的角速度指标进行调整是比较方便的。而对于实际上角速度变化范围并不是特别大的大多数算例来说,可以一个固定的线性算子贯穿始终,这在应用简便性和计算经济性上是比较理想的。结合旋转柔性梁的动力学特征,以最大角速度时对应的Jacobian 矩阵作为线性算子。此时,指数积分法底层方程的刚度略偏大而非偏小,从稳定性的角度来看,一般是较为合适的。

4 数值仿真

本节通过数值实验来验证上述的线性算子选取策略。以系统大范围转动已知的旋转柔性梁这一典型问题为数值算例[24],从仿真精度、稳定性以及CPU计算时间3 个方面进行对比。梁的物理参数以及大范围转动的运动规律均根据文献[24]给定。其中,长度L为8 m,梁的横截面积S为7.2968 ×10-5m2,惯性矩I为8.2189 ×10-9m4,密度ρ为2.7667 ×103kg/m3,弹性模量E为6.8952 ×1010N/m2。系统大范围转动的运动规律为

该运动规律表示旋转柔性梁从静止开始旋转加速至稳态转速的情况。在下面的算例中,取加速过程时间T=15 s,最大转速参数Ω分别取4 rad/s 和15 rad/s。对此均匀梁,以5 个单元离散,此时已有很好的空间精度。该算例是典型的数值刚性问题。

4.1 线性算子对比

为了验证本文线性算子选取方式的合理性,用ETD 格式对上述问题进行数值求解,并比较采用不同的线性算子所获得的数值结果。

对于Ω=4 的数值算例,角速度∈[0,4]。分别取不同旋转角速度对应的线性算子A(0)、A(1)、A(2)、A(3)、A(4)、A(5)、A(6)。其中,A(0)表示以=0 来构造的线性算子,即取方程的线性项矩阵作为线性算子。A(4)表示以最大角速度=4 来构造的线性算子,以此类推。A(1)、A(2)、A(3)是以中间角速度,而A(5)和A(6)是以更大角速度分别构造的线性算子。本算例中,角加速度数值只在0～0.533 之间,构造以上线性算子时角加速度均按初始时取为0,而只考虑角速度对系统动力特性的影响来构造线性算子是合理的。

首先,采用ETD1 分别根据上述的线性算子进行数值计算,考察各参数下的ETD1 的稳定性。表1中对比表中列出的4 个时间步长下的结果。可以看到,在求解时段内,用较小角速度的线性算子A(0)～A(3),只有以步长h=10-4计算时的数值结果没有出现发散。作为对比,较大角速度的A(4)、A(5)和A(6),在所有考察的步长下都没有出现数值发散。这是由于当线性算子A(0)～A(3)没有包含系统主要的刚性、振荡信息,ETD 格式求解就容易发散。而当线性算子包含了全部或者略为偏多的刚性、振荡信息时,格式相对不容易发散。

表1 采用ETD1 求解的CPU 计算时长(仿真时长50 s)

由上可知,通过选取合适的线性算子,指数积分法就可以放大步长进行数值计算,从而提高求解效率。当然,求解的步长能够放大到多少,不仅取决于稳定性条件,也取决于解本身的变化速度。

图3 是梁末端位移的参考解。可见,相对于求解参考解所用的时间步长h=10-5来说,解本身的变化是极为缓慢的。从时间分辨率来说,并不需要如此小的步长,h=10-2大概已经足够。下面,考察以步长h=10-2求解时的数值误差情况,来比较ETD1 用以上各线性算子得到的数值结果。

图3 参考解,以四阶Runge-Kutta 法求解,时间步长h=10-5

图4(a)是以A(0)～A(3)求解时的数值误差。从图中看数值解都是发散的,发散的趋势随着线性算子呈现规律性,即线性算子中选用的角速度越接近最大角速度,数值解的发散时刻逐渐推迟。可见,线性算子的选取对ETD 格式的发散速度会产生影响,较接近最大角速度的线性算子表现相对较好。图4(b)是以A(4)～A(6)求解时的数值误差。可以看到,A(4)的误差最小,而A(6)的误差最大。单从误差上看,A(5)和A(6)的数值误差随时间变大,似乎有发散的趋势。但结合图5(a)的解曲线可知,误差增大是由于A(5)、A(6)的数值结果呈现衰减趋势。实际求解到2000 s,两者仍未出现数值发散。在使用更大步长h=10-1后,数值衰减现象会更明显,如图5(b)所示。这是相应ETD 格式的数值阻尼效应。可见,当采用大于最大角速度的A(5)和A(6)时,在整个求解时间内线性算子都包含着过强的刚性、振荡信息,虽然在有限时间内可能不会出现数值发散,但长时间的误差累积会影响数值解的精度。

图4 末端位移误差图

图5 旋转柔性梁的末端位移

综上可知,采用ETD 格式数值求解非惯性系下的旋转柔性梁系统,线性算子需要根据角速度来构造。从稳定性的角度考虑,以最大角速度时的Jacobian 矩阵作为线性算子是相对合适的。这样选取的线性算子可以兼顾稳定性以及精度。若无法精确取到相应线性算子,与偏小角速度相比,偏大角速度的线性算子稳定性相对来说会更好一些,虽然在精度上会由于格式的数值阻尼而损失。

4.2 不同阶ETD 格式的比较

本节中,对比不同ETD1、ETDRK2 以及ETDRK4 的数值结果。数值算例取Ω=15,此时,∈[0,15]。线性算子分别取A(0)、A(15) 2 种来进行计算。仿真时长为50 s。其中,A(0)表示以=0 来构造的线性算子,即取方程的线性项矩阵作为线性算子;A(15)是以最大角度构造的线性算子。

表2 是各阶ETD 格式在时间步长h=10-2下计算所花的CPU 时间,其中,出现发散的情况用*表示。由表2 可知,就CPU 计算时间而言,ETD1 的计算效率几乎是ETDRK2 的2 倍,是ETDRK4 的4 倍。另外,在以A(0)作为线性算子时,各阶ETD 格式均出现了数值发散;以A(15)作为线性算子时,各阶ETD 格式在表中所列的计算时间内都没有出现数值发散。这一结果,与本文的线性算子选取建议和前一算例的结论相符。

表2 仿真时间50 s,各格式的数值结果,步长h=10 -2

图6 是旋转柔性梁末端位移的数值误差随时间的变化情况。图6(a)是以A(0)作为线性算子时的数值误差,可以清楚地看到ETD1、ETDRK2 以及ETDRK4 先后出现数值发散。图6(b)是A(15)作为线性算子时的数值误差,可以看到各阶ETD 格式均没有出现数值发散。从误差的数值上来看,ETD1格式的误差最大,而ETDRK2 和ETDRK4 的误差相差不大。

图6 末端位移绝对误差图,步长h=10 -2

图6 中出现一个“反常”的现象。ETDRK4 的精度阶次比ETDRK2 高,但在图6(b)中,后者的误差反而比前者的小。当0≤t≤15 时,旋转柔性梁处于加速旋转状态,此时,固定的线性算子始终与包含实时刚性信息的Jacobian 矩阵会存在偏差。作为高阶格式,ETDRK4 相对于ETDRK2,对线性算子上的偏差应该更为敏感。因此,在加速阶段0≤t≤15,有可能在一定的时间段内,ETDRK4 格式是弱不稳定的。由于时间段不长,不会引起数值发散,但导致数值解在此阶段累积了较多误差。于是,出现了二阶格式ETDRK2 比四阶格式ETDRK4 的误差更小的现象。随着时间步长的缩小,ETDRK4 在加速阶段的弱不稳定情况逐渐消失,逐渐呈现收敛性。例如,图7 中的ETDRK4 的精度远高于ETDRK2。对于本问题,采用ETDRK2 格式的CPU 计算时间,仅是采用ETDRK4 时的一半。综合考虑,对实用性仿真来说,ETDRK2 是更适合用来求解该问题的ETD 格式。

图7 末端位移绝对误差图,步长h=10 -3

5 结论

本文针对非惯性系下旋转柔性梁的数值仿真,给出了指数时间差分法求解该问题时的应用策略。结合旋转柔性梁的动力学特征,将其大范围转动的角速度纳入指数积分法线性算子的构造。由于动力刚化现象,系统的刚度会随角速度变化,线性算子需要根据大范围转动的角速度来选取。对于本文中的数值算例,以上述的线性算子选取方式,指数积分法的稳定性得到了改善,使得各阶ETD 格式能够采用更大的时间步长来求解问题。这能够提高旋转柔性梁数值仿真的求解效率。其次,本文对比了不同阶的ETD 格式,从精度、计算效率等方面综合考虑,指出对于一般问题采用ETDRK2 进行数值求解可能更适合。

值得指出的是,本文中所考察的旋转柔性梁系统,其大范围运动相对简单。如果采用本文的线性算子构造策略,单一固定的线性算子就能贯穿整个计算过程而获得较好的数值结果。但是,对于系统大范围运动未知、角速度变化范围较大等更复杂的情况,更好的做法是阶段性地根据角速度来动态选取线性算子。而本文的相关结论可以为进一步的研究提供基础。