时间谱元法在动态响应优化中的应用

本资料包含pdf文件1个，下载需要1积分

机械结构几乎都在动载荷作用下运行，机器的各种动态性能是与时间相关的函数。欲使机器的性能最优，动态优化是很有必要的。故而必须处理由于动态载荷作用机器的各种动态响应，包括 目标函数与约束函数。在时间域上，较精确地求出相关响应，并且要求满足时间函数约束，而且不能漏掉任何响应，特别是响应极值；而应用很小步长计算动态响应，-定程度上可以减小时间函数约束的可能。在优化过程中，每-次迭代都要重新计算时间函数目标函数、时间函数约束，因此，要求求解动态响应必须高精度而且高效。

有限差分法是比较流行的求解动态响应微分方程的方法。应用很小时间步长计算二阶微分方程，自适应调节时间步长，直到获得要求的精度。时间步长直接影响计算的稳定性和准确性，这极大地影响了计算效率。特别是当结构受脉冲载荷时，差分法很难获得满意的结果。文献[1]应用谱元法的 h型精细方案对受脉冲的动态系统分析求解，在没有增加单元数量的基础上获得了谱收敛精度。

文献[2]中，提出将设计变量和位移响应、设计变量和位移响应与速度响应、设计变量和位移响应与速度响应加速度响应分别作为优化变量 3种方案，应用有限差分法近似和时间相关的约束，运动微分方程被作为等式约束处理。对于动态响应最大值最小问题，文献E33提出了直接处理 目标函数的方法，解决了迭代过程中最大响应震荡造成收敛困难的问题。文献[4]应用时间有限元法将动态响应微分方程转化为代数方程，然后进行了二阶灵敏度分析，为动态响应优化奠定基矗文献[5]中，应用了时间谱元法离散时间微分方程，GLL点法和关键点法处理与时间相关的约束，然而 GLL点法不稳定，对单元数和插值次数比较敏感，关键点法由于每-次迭代的最大值的位置发生变化，因此收敛速度非常慢甚至是发散 。

20世纪 70年开始研究承受瞬态载荷的结构优化。2006年，Kang等首次把动态响应优化作为优化的-个分支 ]。在其优化中，设计变量每迭代-次，时间响应需要计算-次，同时约束必须在整个时间区间内被满足。处理时问约束有几种方法。-种是只在响应的全局最大值处满足约束；另-种更加稳定的方法是使在更小时间步长上满足约束，保证中间点处的约束失效不 可能发生。在 准静态方法中，应用了这种方法使多个等静态载荷满足约束而不是动态载荷7]。在这些方法中，约束数大大增加。

因为在优化迭代过程中，要计算这些约束的灵敏度，收稿 日期 ：2012-03-26；修 订 日期 ：2012-11-14基金项目：国家自然科学基金资助项 目(51275489，51175198)；中北大学校基金资助项 目(2012)；博士科研启动基金资助项 目(2011)g 撮-学E程工V动m振振 动 工 程 学 报 第 26卷这样优化耗费也在增加。再-种更加有效的方法是只在响应的局部极值点处处理约束。这样可以减少约束的数量和灵敏度的次数 ]。

本文提出关键点及其相邻的 GLL点法处理与时间相关的约束，在最大值附近的-些点能够被包括在约束集，克服当前迭代的局部最大值可能在下- 次迭代不是局部最大值的问题。

1 时间谱元法1984年 Patera提出了谱元法g]，在流体动力学中得到应用，其融合了有限元法的处理边界、结构的灵活性与谱方法的快算收敛性。在要求相同精度的情况下，谱元法能够采用较少的单元减小了计算开销。在-个单元内，将时间离散为与 GLL多项式零点相对应的网格点，在这些点上进行 Lagrange插值。

从理论上分析，在-定点数上插值，当这些点是对应的正交多项式的零点时获得的插值精度最高l 。如图1所示，均匀点分布和 GLL点分布分别近似 Runge函数式。

1厂(z)- ，z∈ [-1，1] (1)1.0 -0.5 0.0 0.5 1.O图 1 GLL点和平均点 Lagrange插值 Runge函数Fig.1 Lagrange interpolation in GLL points and averagepoints for Runge function在有限元法中，随着插值次数增加，近似的数值误差极大或近似失败；而 GLL插值有明显的优势。

时间谱元法的主要步骤：(工)首先转化二阶微分方程组为-阶微分方程组，应用 Bubnov-Galerkin原理转化为等效的积分形式，代入单元积分表达式，得到单元谱元方程。(I)将时间段划分成若干单元，每个单元由若干个节点组成。(1I)用 LegendrePIOGLL点图 2 GLL点分布(p表示单元内插值节点数，在[-1，11区间上)Fig.2 GLL point distribution(The represents the numberof interpolation nodes.Intervels are[-1，～1])正交多项式作为基函数，将每-个单元的近似解表达为基函数的线性组合。(IV)将所有单元的谱元方程按-定规则合成其总体谱元方程。(V)解时间的总体谱元方程，得到其全局近似解。

1.1 动态响应微分方程及其形式转化动态响应微分方程可表达为4-d 4-Kx-F (2)式中 M 为质量矩阵，c为阻尼矩阵，K为刚度矩阵，F为时间函数载荷向量，X为位移向量，呈为加速度向量， 为速度向量。初始条件 (O)-b。，x(0)-1，。；其中，M，C，K不随时间变化；X，王， 为时间的函数，F为时间的任意函数，t∈Et。，tn]。

在时间谱元法中，为了用数值方法解动态响应方程 ，通过设 - ，X - ，将式(2)转化为-阶线性常微分方程组式M0 。K0] X- 0 F0)， . 十l l l 1 l l zJ LIL-I Jlz JLIJIJ1(O)- ，X2(O)：X (3)1.2 单元划分如图3所示，将仿真时间t∈Et。，t ]分割为 T/个互不相交的单元，即[ ，t ]，Et ，t ]，[t。，t。]，，Et ，t ]，对每单元配置若干个点。在配置点时，可以根据动态载荷的变化情况，在载荷突变的地方配置多-些点，在载荷平缓的地方可以配置少-些点。这样能保证在动态载荷的所有地方都能有足够的精度。

tL-------J-----J ， -、、- - ， - - -J-------J图 3 时间谱元法离散时间Fig.3 Divided the time domain temporal SEMO 5 0 5 O 5 2 l 1 O O 0 第 3期 毛虎平，等：时间谱元法在动态响应优化中的应用1.3 单元分析区间分段以后 ，将每-个时间段划分单元 ，为提高近似解的精度，可以在单元区间上再增加插值点来增加近似解的阶次 ，这在有限元里被称为高阶元 ，即 hp方法。

在时间谱元法中，为了达到谱方法的收敛速度，并且增加数值计算的精度，单元内插值点应放置在特殊的配置点上，它们是 由 N1个 GLL点构成的，而基函数由有限项正交多项式的和来表示，构成单元上任意配置点的形函数，如此可在有限的插值点上获得指数级的收敛速度。

在每-个单元上，将定义在插值点上的形函数通过 Lagrange插值得到近似函数。换句话说，任意- 个定义在参考单元上的函数 ( )，都可以用下式近似， ( )- ：z， ( )P：J)( ) (4)式 中 P ( )为 J单元 的忌次 Lagrange多项式 ，为定义在[-1，1]上的 GLL点，z ( )为单元 J上未知节点在 GLL点的值。

有两种类型正交多项式可用于时间谱元法中：分别为 Chebyshev和 Legendre正交多项式。这里只讨论 Legendre展开，勒让德多项式是 Sturm-Liouville Form方程的解[(1- )P ( )] n(n1)P( )-0 (5)它 的权函数 叫-1。忌阶 Legendre多项式可 以定义如下P ( )- 广dE(g- ) ] (6)由于 Legendre多项式的零点不包括区问端点，因此引入 Lobbato多项式。Lobbato多项式是通过Legendre多项式的微分定义的口。

L。 ( )- P 1( ) (7)将式(3)展开，其中-个方程可以表达为A sz 厂(纠- (8) o-d 0，z0 0，t∈Et ，ti1式中 A 为耦合矩阵。

将每-个状态变量在给定时间段内离散化，近似为m次 Lagrange多项式(4)，在每-个单元上应用 Bubnov-Galerkin法使得插值误差最小1 ，得出霎 [ -厂 ( ， ))ld -0 (9)将每-个单元用矩阵形式表示为L ( )- F。(t) (10)式中 L ： -A I ，F ( )- -J f ， ([ ( )z( ) .27( )] ) ，L 为 的-般函数，∞ 。 为Gauss-Lobatto积分在 q点的权值口 。

1.4 集成总体时间谱元方程并求解总体合成的任务是将所有单元谱元方程按序叠加，得到总体谱元方程，即线性方程组。对于 N 个状态变量，通过耦合矩阵A (N ×N 方阵)的张量叉乘得到全部状态变量的全局组装式(r B )x - (A： B )x -(r B )F (x ) (11)式中 B ，B 为第 i时间单元的全局微分矩阵与全局权矩阵，F (x )为第 i个时间单元的激励力的全局形式 ，I 为 N ×N 单位矩阵 ，x 为第 i个时间单元的所有状态变量在时间节点的集合。

化简式(11)得到G X - - B F(X ) (12)其中，G -B mn i。未知向量 x 包含了第i个时间单元的所有单元中所有 GLL配置点上的位移和速度离散解。 为第 i个时间单元的全局线性矩阵。式(12)所示的线性代数方程组可以直接求解。

2 动态响应优化机械结构动态响应优化设计问题的目标函数为选择设计变量向量 x(z ，z ，，z ) ，以使系统在受到瞬变载荷作用下，在给定时问区间[0，T]内满足瞬时动态性能(振幅或相对位移极限、动应力、动应变破坏极限等)及设计变量允许变化范围等约束条件，并使系统某些关心位置或坐标的最大动态响应(位移、速度、加速度)在某种意义或准则下达到最优。其数学模型可表示为X - (z1，.2C2， ，z )z( )-Ez(t)，乏，兰]之- x，t， ( )，F( )]Z(to)z ZO， (。)- 。

(13)J(X ， ， )h Ex，t，z( )]-0 i1，2，，mg x，t， ( )]≤ 0 J-1，2，，V t∈Eo，刀式中 x为系统的设计变量向量，由系统的几何参数和物理参数组成；z( )为系统的状态变量向量，由系统运动状态的广义坐标如位移、速度、加速度等组成，其要满足机械系统运动规律的运动方程，即动态特性的数学描述，之(f)-pEx，t，z，F]；J为目振 动 工 程 学 报 第 26卷标函数；h ，g，分别表示等式约束与不等式约束，它们在所有时间点(V t∈E0，T])上都要满足。

处理与时间相关的约束有两种方法(见文献[12])：-种是约束在所有 GLL点得到满足；另-种是约束在每-个单元的绝对值极值点得到满足。第- 种方法要求 GLL点之间的距离尽量的小，这样可以减少漏掉 GLL之间的绝对值极值点得不到满足的可能性，因此约束数量比较庞大，同时谱元法的求解开销也较大；第二种方法可以用满足精度要求的较少谱单元求解运动微分方程，在每-个单元内，对其高次 Lagrange函数进行-维搜索找到单元绝对值极值点，当然单元绝对值极值点是震荡的，因此每迭代-步，单元绝对值极值点要重新计算。文献[5]中说明了 GLL点法处理与时间相关的约束，非常不稳定，但单元数和插值点数达到-定大时，才能找到最优值，否则找不到最优值；而关键点法，当前迭代的局部最大值可能在下-次迭代不是局部最大值，这样收敛速度变慢甚至不收敛。本文提出将所有单元响应绝对最大值以及其相邻的 2个 GLL点作为约束，在最大值附近的-些点能够被包括在约束集，如图4所示。划分 3个单元，即[t。，t ]，[t ，t ]，It。，t。]，关键点分别为 P ，P。，P。，其相邻的 GLL点分别为P1，P12， 21，P 2，P3l，P 32。

g lI0 /f0 、/ /f2 i / l图 4 单元绝对极值点及其相邻的两个 GLL点作为约束Fig.4 The constraints are the absolute extreme valueand the two adjacent GLL points这样，关键点集 P 可以采用下式表示Pc-t i / ，tGLLo，tGLI 1 (14)式中 t -argminL，tGIJL。为 t 。 /ma 的左端附近 GLI 点 ，tG ，. 为 t /ma 的右端 附近 GLL点 。从式(14)中可以很清楚地看出本文与作者前期研究的不同，关键点法仅仅考虑 train/max，而没有考虑左右端的点。

3 算例分析在各种各样 的机械系统 (如 ，减振器、武器后座机构、飞机起落架等)中，系统的主质量在该系统达到稳态之前可能作大幅度振动，当激励力的频率接近于主系统的固有频率时，系统可能发生破坏，对激励的瞬时动态响应必须约束。对于这种问题，可以用抽象模型的优化来解决。

3.1 线性两自由度减震器最优设计已知两 自由度系统，参数如图 5所示，m -4.534 kg，m20.453 4 kg， 1-1 749.03 N/m，激励频率是主质量固有频率 n 的 1.2倍，即 -23.57。求阻尼器的刚度系数 k 和阻尼系数 c使得主质量的最大位移响应最小，同时满足振动空间约束条件减振器位移响应和主质量位移响应之差不能超过 3倍的主质量最大响应，以及稳态约束条件。

图 5 减振器Fig.5 Shock absorberfm1-z1( )k1Lzl( )k2[z1( )-322]1 c[z1( )-z2( )]-Fsincotm2星2k2[ 2～321]c[主2( )- (15)f . 1Xl)]-. . 2121(0)-z2(O)-o，321(O)-322(0)-0引入符号 ：z -/F/忌 是由力 F产生的主质量的静位移，n -/愚 m 为主系统的非耦合固有频率， -/走 /m为阻尼器系统的非耦合固有频率， m。/m 为减振器质量与主质量之比，厂- /n 为减振器与主质量的非耦合固有频率之比， -叫/ 为激励频率与主质量的非耦合固有频率之比，C -2m。 为临界阻尼，c为阻尼系数， -c/C 为阻尼比，z -z ( ，厂， )为主质量的位移。令 -32 /32 ，i-1，2，式(15)转化为fz1Q： 12 (z1-z2) (z -z2)- ：sin(J2 ) (16)2w (乏 -乏 )l厂 (z。-z )-0再令z1-乏1，z2-zl，323-乏2，z4- 2，那么式(16)转第 3期 毛虎平 ，等：时间谱元法在动态响应优化中的应用 399化为 A X - F那么，优化模型表示为f min m-a x I 2( ，f， ) l t∈[0，刀(17)S.t. A X - Fl ( ，f， )-lzz( ，f， )l≤3max f zz( ，厂，) (18)X4 ( ，f， )≤ 口z ( ，f，)-z z( ，f，)I≤ 3口mi ≤ ≤ mi ≤ 厂≤ -厂根据文献E12]，引入人工设计变量 z ，可以把式(18)转化为o(X)- z5S.t.X A X - F-z2 l≤ z5l z ( ，厂 -x2( ，厂， )l≤ 3x5(19)z4 ( ，f，)≤ 口J z ( ，f， )-.27 2( ，f，t)l≤ 3≤ ≤ 色t≤ f≤ 式中 n-0.030 48 m；z ，z 分别表示 主质量和阻尼器的稳态位移响应，表达式为仨 进2：± 二 z6兰堡 !√ z6 (2O)式中 Xb-(2 )。( -1) [ 厂 -( -1)( -厂。)] ，阻尼比 ∈[io～，0.167 85]，减振器与主质量的非耦合固有频率之比.厂∈[1O-，2.O]，减振器质量与主质量之比 -0.1不变。应用黄金分割法进行-维搜索寻找关键点，将 目标函数、时间设计变量的收敛上限设置为 1.0×10 ；为了不漏掉-切可能的关键点，在每-个单元进行两次-维搜索，那么要求找到的所有单元关键点满足约束条件。以其中-组数据为例，优化迭代过程如图 6所示。由文献[13]中知，梯度投影法计算的已知条件- 0.1，a- 1.2，10 ≤ f≤ 2.0，lO ≤ ≤0.167 85， -1.2，初始条件f-1.6， -0.02，d-0.081，最优结果 .厂-1.338， -0.021 21，d0.060 1。

表 1和2为与文献[13]优化结果比较，可以看出本文提出的关键点及其相邻的 GLL点法处理与时间相关的约束非常稳定，对不同的单元和插值数都找到了最优值，而且迭代次数 比较少，最大的 13次；而在作者的前期研究(文献[5])中关键点法处理与时间相关的约束不稳定，比如当 Nel4， -5以及当Nel10， -6时，没有找到最优值，并且迭代次数是关键点及其相邻的 GLL点法的两倍多。两自由度问题是最简单的多 自由度问题，这样就能说明本文研究结果更具有实际意义 。

其实，从耗时来分析，还可以给 SQP提供解析梯度函数或者简单的数值梯度函数，这样函数的评估次数能明显减少口 ]；SQP算法采用差分法计算梯度，每次迭代需要多次进行仿真函数评估(评估次表 1 关键点及 其相邻的 GLL点 法数值 实验结果Tab.1 Numerical experimental results of the critical points and adjacent GLL points203O40501.328 896 27O.O22 616 0O0.060 040 141.328 857 69O.O22 91O O7O.O59 882 831.328 857 470.O22 911 73O.O59 881 971.328 857 46O.O22 911 79O.O59 881 911.328 857 46O.O22 911 820.059 881 89881.328 857 46O.O22 911 82O.O59 881 921.328 857 46O.O22 911 82O.O59 881 891.328 857 46O.O22 91l 82O.O59 881 891.328 857 46O.O22 911 82O.O59 881 891.328 857 46O.O22 9l1 82O.O59 881 891171.328 857 460.022 911 82O.O59 881 891.328 857 46O.O22 911 820.O59 881 891.328 857 460.022 911 82O.O59 881 891.328 857 46O.O22 911 82O.O59 881 891.328 857 46O.O22 91l 820.059 881 899400 振 动 工 程 学 报 第 26卷表 2 关键点法数值实验结果Tab.2 Numerical experimental results of the critical points0.09O.08目§0.07撂0.060.050.O40.03O 5 10 15 2O 25 30迭代次数n/次图 6 主质量动态响应Fig.6 The dynamic response of main mass数等于设计变量维数的两倍乘以实例数1)，优化求解很费时口 ；对于基于谱元法的动态响应优化，每评估-次函数值就要运用谱元法更新-次，很明显，单元数越多及其插值次数越高，耗时越多，减少函数评估次数，自然耗时就能下降。从文献[12]可知，作者研究的单 自由度的缓冲器优化结果波动的原 因是是否有网格点在 响应的绝对值最大点处，文献E]23第 1O页图 4说明了此问题。因此这里只比较了关键点法与本文提出的方法。表 3是