AUXILIARY HARMONIC EXCITATION GENERALIZED METHOD FOR RANDOM VIBRATION ANALYSIS OF LINEAR STRUCTURES UNDER STATIONARY GAUSSIAN EXCITATION
-
摘要: 相比于时域法, 频域法是更为高效、易行的随机振动分析方法, 但对于平稳激励下的随机振动分析, 现有频域方法常需振型截断或功率谱矩阵分解, 将会影响计算精度和效率. 为此, 本文在频域法的框架下, 针对平稳高斯激励下线性结构的随机振动分析提出了一种精确且高效的辅助简谐激励广义法. 首先, 引入广义脉冲响应函数和广义频响函数的概念, 推导了与响应功率谱计算的完全二次项组合法等价的广义分析方法. 其次, 通过辅助简谐激励的响应乘积代替广义频响函数的乘积, 在广义分析方法的基础上进一步提出了更易于实现的辅助简谐激励广义法. 再次, 根据辅助简谐激励下结构响应求解方式的不同, 提出了具有不同适用性的两种辅助简谐激励广义法实现方案, 即基于振型叠加的辅助简谐激励广义法和基于时程分析的辅助简谐激励广义法; 同时, 给出了上述两种实现方案的计算性能及其与已有方法的对比分析. 最后, 通过两个算例验证本文所提方法的计算精度和效率. 由算例结果可知, 本文提出的辅助简谐激励广义法在计算响应功率谱时与完全二次项组合法和虚拟激励法的计算精度保持一致, 而计算效率相对完全二次项组合法和虚拟激励法具有显著的优势.Abstract: Compared with the time-domain method, the frequency-domain method is a more efficient and easy-to- implement method for random vibration analysis. However, the existing frequency-domain methods often involve truncation for degree of mode or decomposition of the power spectrum in multi-correlation conditions, which may have impact on the computational accuracy and efficiency of the methods. To this end, an accurate and efficient auxiliary harmonic excitation generalized method is proposed for the analysis of random vibration of linear structures under stationary Gaussian excitation in the framework of the frequency domain method. First, the concepts of generalized impulse response function and generalized frequency response function are introduced, and a generalized analysis method, which is equivalent to the complete quadratic combination method of response power spectrum calculation, is derived. Secondly, replacing the product of generalized frequency response function by the product of response of auxiliary harmonic excitation, a more easily implemented auxiliary harmonic excitation generalized method is further proposed based on generalized analysis method. Third, according to the different calculation methods of response for structure under the auxiliary harmonic excitation, two generalized methods of auxiliary harmonic excitation generalized method with different applicability are proposed, namely, the auxiliary harmonic excitation generalized method based on the mode superposition and the auxiliary harmonic excitation generalized method based on the time analysis. Meanwhile, the computational performance of the above two methods and their comparative analysis with the existing methods are introduced. Finally, the computational accuracy and efficiency of the proposed method are verified by two examples. The results of the examples show that the auxiliary harmonic excitation generalized method has significant advantages of the calculation efficiency over the complete quadratic combination method and the pseudo-excitation method with the same calculation accuracy.
-
引 言
线性结构在平稳高斯随机激励下的响应分析是随机振动研究的基础[1-2]. 线性随机振动理论分析方法可以分为时域分析法和频域分析法[3-6]. 时域分析法是以激励的自相关函数为基础的分析方法, 比如矩方程法[2,7]、递推矩法[8]等. 频域分析法是以激励的谱密度为基础的分析方法. 相对于时域分析法, 频域分析法数学形式简单, 因此成为了随机振动领域的主流方法. 频域分析方法根据输出结果形式的不同主要分为蒙特卡洛法[9-10]、均方根法[11-15]和功率谱法[16-19]. 由于功率谱法是频域法中最基本的分析方法, 所以本文主要针对功率谱法开展研究.
功率谱法是直接建立激励功率谱和响应功率谱之间关系的分析方法. 根据两者关系表达式的不同, 功率谱法可以分为振型法和虚拟激励法. 振型法顾名思义是以振型为基础计算响应功率谱的分析方法. SRSS法[16]作为振型法最为原始的分析方法主要是通过振型叠加法给出响应功率谱的计算表达式, 该方法由于不考虑振型之间的耦合项而有较高的计算效率. 虽然SRSS组合简单, 但当模态频率接近时, 其计算结果将存在较大误差. 随后, CQC法[4-5]通过考虑振型之间的耦合项给出响应功率谱的精确表达式. 然而, 针对一些需要考虑较多振型的结构, CQC法将面临计算量大的问题[17]. 为提高计算效率, 引入振型截断是最为常见的处理方式[18-20]. 但振型截断往往对计算精度造成不利影响[21], 且振型的截断阶数并没有一个合适的确定标准. 对于大型复杂结构, 即使引入振型截断, 仍然会涉及大型矩阵的运算进而影响计算效率[22].
虚拟激励法[23-28]首先构造与激励相关的虚拟激励, 然后通过虚拟激励下的响应获得响应功率谱. 对于复杂结构且频率离散点数目较大时, 需要计算每个频率离散点处的虚拟激励响应, 计算量有时也比较大. 值得指出的是, 针对多点激励问题虚拟激励法需要引入激励功率谱矩阵的Cholesky分解将多点相关激励转换为独立随机激励. 一方面, 当频率离散值数目和激励数目较多时, Cholesky分解计算量较大; 另一方面, 当矩阵不正定性, Cholesky分解将难以实施[29]. 换言之, 虚拟激励法尽管可以避免CQC法中振型计算及其截断导致的误差, 却衍生了功率谱矩阵分解影响计算效率的新问题.
不难发现, 现有的功率谱法常常涉及振型截断或激励功率谱矩阵分解带来的精度和效率问题. 基于此, 本文在引入广义频响函数的基础上, 推导了响应功率谱的计算表达式, 结合辅助简谐激励策略提出了一种新型随机振动分析方法, 即辅助简谐激励广义法. 同时根据辅助简谐激励下响应求解方式的不同, 给出了具有不同适用范围的两种求解方案.
1. 平稳激励下线性结构响应的广义分析方法
1.1 理论基础
对于s个自由度的线性结构, 在平稳随机激励下结构系统的动力微分方程可写为[30]
$$ {\boldsymbol{M}}\ddot {\boldsymbol{Y}}\left( t \right) + {\boldsymbol{C}}\dot {\boldsymbol{Y}}\left( t \right) + {\boldsymbol{KY}}\left( t \right) = {\boldsymbol{QF}}\left( t \right) $$ (1) 式中, M, C和K分别为结构的s × s维质量矩阵、正交阻尼矩阵和刚度矩阵, Q为s × m维激励定位矩阵, F(t)=[F1(t), F2(t), ···, Fm(t)]T为平稳随机激励向量, 上标T为转置, m为激励的数量.
在零初始条件下, 结构第k自由度的稳态位移响应Yk(t)为
$$ {Y_k}\left( t \right) = \sum\limits_{j = 1}^s {\int_{ - \infty }^\infty {{\varphi _{k,j}}{h_j}\left( \tau \right){\boldsymbol{\varphi}} _j^{\rm{T}}} {\boldsymbol{QF}}\left( {t - \tau } \right){\rm{d}}\tau } $$ (2) 式中, φj=[φ1,j,φ2,j,···, φs,j]T表示第j阶振型, φk,j为φj中第k个元素, hj(t)为脉冲响应函数, 可表示为
$$ {h_j}\left( t \right) = \frac{{\text{1}}}{{{\omega _j}\sqrt {1 - \xi _j^2} }}\sin \left( {t{\omega _j}\sqrt {1 - \xi _j^2} } \right){{\rm{e}}^{ - {\xi _j}{\omega _j}t}} $$ (3) 其中, ωj和ξj分别为结构第j阶频率和振型阻尼比.
1.2 基于广义脉冲响应函数的响应计算
不难发现式(2)涉及所有振型的求和, 如果引入广义脉冲响应函数, 则响应求解将转化为关于荷载数目的求和.
将F(t)中第l个元素用脉冲函数δ(t)代替, 其余分量都取为0, 即Fl(t)=[0, ··· , 0, δ(t), 0, ···, 0]T. 在Fl(t)作用下的运动微分方程可以写为
$$ {\boldsymbol{M}}\ddot {\boldsymbol{Y}}\left( t \right) + {\boldsymbol{C}}\dot {\boldsymbol{Y}}\left( t \right) + {\boldsymbol{KY}}\left( t \right) = {\boldsymbol{Q}}{\underline {\boldsymbol{F}} _l}\left( t \right) = {{\boldsymbol{Q}}_l}\delta \left( t \right) $$ (4) 其中, Ql为Q中的第l列向量. 结合式(2), 可得式(4)中第k自由度的稳态位移响应Yk,l(t)为
$$ \begin{split} {Y_{k,l}}\left( t \right) =& \sum\limits_{j = 1}^s {\int_{ - \infty }^\infty {{\varphi _{k,j}}{h_j}\left( \tau \right){\boldsymbol{\varphi}} _j^{\rm{T}}} {{\boldsymbol{Q}}_l}\delta \left( {t - \tau } \right){\rm{d}}\tau } = \hfill \\ {\text{ }} & \sum\limits_{j = 1}^s {{\varphi _{k,j}}{h_j}\left( t \right){\boldsymbol{\varphi}} _j^{\rm{T}}{{\boldsymbol{Q}}_l}} \buildrel \Delta \over = {\underline h _{k,l}}\left( t \right) \hfill \end{split} $$ (5) 显然, hk,l(t)为结构在脉冲激励下的响应, 亦可称之为脉冲响应函数, 但hk,l(t)所对应的脉冲激励是作用在与激励分量Fl(t)相应位置上的, 与常规的脉冲响应函数不完全相同, 因此, 本文将hk,l(t)称为广义脉冲响应函数. 由式(5)可知, 广义脉冲响应函数hk,l(t)包含了所有脉冲响应函数hj(t)的贡献, 且与F(t)中第l个分量的位置向量Ql有关.
于是, 结构在荷载Fl(t)=[0, ··· , 0, Fl(t), 0, ··· , 0]T作用下, 第k自由度的响应为
$$ {Y_{k,l}}\left( t \right) = \int_{ - \infty }^\infty {{{\underline h }_{k,l}}\left( \tau \right){F_l}\left( {t - \tau } \right){\rm{d}}\tau } $$ (6) 相应地, 结构在激励F(t)下Yk(t)可表示为
$$ {Y_k}\left( t \right) = \sum\limits_{l = 1}^m {\int_{ - \infty }^\infty {{{\underline h }_{k,l}}\left( \tau \right){F_l}\left( {t - \tau } \right){\rm{d}}\tau } } $$ (7) 比较式(2)和式(7)可知, 引入广义脉冲响应函数后, 结构响应表达式的求和次数由s次转变为m次. 工程实际中, 激励数量m通常小于结构自由度s. 因此, 若可以方便地给出hk,l(t), 则式(7)比式(2)更易于实用.
1.3 基于广义频响函数的响应功率谱计算
结合式(7), Yk1(t)和Yk2(t)的互相关函数
${R_{{Y_{{k1}}}{Y_{{k2}}}}}\left( \tau \right) $ 可表示为$$ \begin{split} {R_{{Y_{{k1}}}{Y_{{k2}}}}}\left( \tau \right) =& E\left[ {{Y_{{k1}}}\left( t \right){Y_{{k2}}}\left( {t + \tau } \right)} \right]= \hfill \\ & \sum\limits_{p = 1}^m {\sum\limits_{q = 1}^m {\int_{ - \infty }^\infty {\int_{ - \infty }^\infty {{{\underline h }_{{k1,p}}}\left( {{\tau _p}} \right){{\underline h }_{{k2,q}}}\left( {{\tau _q}} \right)} } } } . \hfill \\ & {\text{ }} E\left[ {{F_p}\left( {t - {\tau _p}} \right){F_q}\left( {t + \tau - {\tau _q}} \right)} \right]{\rm{d}}{\tau _p}{\rm{d}}{\tau _q} = \hfill \\ & \sum\limits_{p = 1}^m {\sum\limits_{q = 1}^m {\int_{ - \infty }^\infty {\int_{ - \infty }^\infty {{{\underline h }_{{k1,p}}}} } } } \left( {{\tau _p}} \right){\underline h _{{k2,q}}}\left( {{\tau _q}} \right) \cdot \hfill \\ & {\text{ }} {R_{{F_p}{F_q}}}\left( {{\tau _p} + \tau - {\tau _q}} \right){\rm{d}}{\tau _p}{\rm{d}}{\tau _q} \hfill \end{split} $$ (8) 式中, E(·)为期望运算, RFpFq(·)为激励分量Fp(t)和Fq(t)的互相关函数.
根据维纳−辛钦关系, Yk1(t)和Yk2(t)的互功率谱函数
${S_{{Y_{{k1}}}{Y_{{k2}}}}}\left( \omega \right) $ 可由$ {R_{{Y_{{k1}}}{Y_{{k2}}}}}\left( \tau \right) $ 的Fourier变换得到$$ \begin{split} & {S_{{Y_{{k1}}}{Y_{{k2}}}}}\left( \omega \right) = \int_{ - \infty }^\infty {\left[ {\sum\limits_{p = 1}^m {\sum\limits_{q = 1}^m {\int_{ - \infty }^\infty {\int_{ - \infty }^\infty {{{\underline h }_{{k1,p}}}} } } } \left( {{\tau _p}} \right){{\underline h }_{{k2,q}}}\left( {{\tau _q}} \right)} \right.} \cdot \hfill \\ &\qquad {\text{ }} {{R_{{F_p}{F_q}}}\left( {{\tau _p} + \tau - {\tau _q}} \right){\rm{d}}{\tau _p}{\rm{d}}{\tau _q}} \Bigggr]{16}{{\text{e}}^{ - {\rm{i}}\omega \tau }}{\text{d}}\tau = \hfill \\ & \qquad \sum\limits_{p = 1}^m {\sum\limits_{q = 1}^m {\left( {\int_{ - \infty }^\infty {{{\underline h }_{{k1,p}}}\left( {{\tau _p}} \right){{\text{e}}^{{\rm{i}}\omega {\tau _p}}}{\rm{d}}{\tau _p}} } \right) \cdot \left( {\int_{ - \infty }^\infty {{{\underline h }_{{k2,q}}}\left( {{\tau _q}} \right){{\text{e}}^{ - {\rm{i}}\omega {\tau _q}}}{\rm{d}}{\tau _q}} } \right)} } \cdot \hfill \\ & \qquad \left[ {\int_{ - \infty }^\infty {{R_{{F_p}{F_q}}}\left( {{\tau _p} + \tau - {\tau _q}} \right){{\text{e}}^{ - {\rm{i}}\omega \left( {{\tau _p} + \tau - {\tau _q}} \right)}}{\rm{d}}\left( {{\tau _p} + \tau - {\tau _q}} \right)} } \right] = \hfill \\ & \qquad \sum\limits_{p = 1}^m {\sum\limits_{q = 1}^m {\left\{ {\left[ {\underline H _{{k1,p}}^*\left( \omega \right)\underline H _{{k2,q}}^{}\left( \omega \right)} \right] \cdot {S_{{F_p}{F_q}}}\left( \omega \right)} \right\}} } \hfill \end{split} $$ (9) 其中, i2=−1, 上标*表示复共轭, SFpFq(ω)为Fp(t)和Fq(t)的互功率谱密度函数, Hk,l(ω)的表达式为
$$ \begin{split} {\underline H _{k,l}}\left( \omega \right) =& \int_{ - \infty }^\infty {{{\underline h }_{k,l}}\left( t \right)\exp \left( { - {\rm{i}}\omega t} \right){\rm{d}}t} = \hfill \\ {\text{ }} & \sum\limits_{j = 1}^s {{\varphi _{k,j}}{\boldsymbol{\varphi}} _j^{\rm{T}}{{\boldsymbol{Q}}_l}} \left[ {\int_{ - \infty }^\infty {{h_j}\left( t \right)\exp \left( { - {\rm{i}}\omega t} \right){\rm{d}}t} } \right]= \hfill \\ {\text{ }} & \sum\limits_{j = 1}^s {{\varphi _{k,j}}{H_j}\left( \omega \right)} {\boldsymbol{\varphi}} _j^{\rm{T}}{{\boldsymbol{Q}}_l} \hfill \end{split} $$ (10) 式中, Hj(ω)为hj(t)的傅里叶变换. 显然, Hk,l(ω)为广义脉冲响应函数hk,l(t)的傅里叶变换, 文中称其为广义频响函数. 不难发现, Hk,l(ω)包含了所有频响函数Hj(ω)的贡献且式(9)避免了关于s的双重求和.
同理, Yk的自功率谱函数
${S_{{Y_k}{Y_k}}}\left( \omega \right) $ 为$$ {S_{{Y_k}{Y_k}}}\left( \omega \right) = \sum\limits_{p = 1}^m {\sum\limits_{q = 1}^m {\underline H _{k,p}^*\left( \omega \right)} } \underline H _{k,{{q}}}^{}\left( \omega \right){S_{{F_p}{F_q}}}\left( \omega \right) $$ (11) 将式(9)与式(11)相结合, 响应的互功率谱矩阵SYY (ω)与激励互谱矩阵SFF (ω)的关系为
$$ {{\boldsymbol{S}}_{{{YY}}}}\left( \omega \right) = \underline {\boldsymbol{G}} {\left( \omega \right)^*}{{\boldsymbol{S}}_{{{FF}}}}\left( \omega \right)\underline {\boldsymbol{G}} {\left( \omega \right)^{\rm{T}}} $$ (12) 式中, G(ω)=[H1(ω), H2(ω), ···, Hs(ω)]T为广义频响函数矩阵, Hk(ω)=[Hk,1(ω), Hk,2(ω), ···, Hk,m(ω)] (k = 1, 2, ···, s).
将式(10)代入式(12)右侧, 可得
$$ {{\boldsymbol{S}}_{{{YY}}}}\left( \omega \right) = {\boldsymbol{\varPhi}} {{\boldsymbol{H}}^*}{{\boldsymbol{\varPhi}} ^{\rm{T}}}{\boldsymbol{Q}}{{\boldsymbol{S}}_{{{FF}}}}\left( \omega \right){{\boldsymbol{Q}}^{\rm{T}}}{\boldsymbol{\varPhi}} {\boldsymbol{H}}{{\boldsymbol{\varPhi}} ^{\rm{T}}} $$ (13) 式中, Φ=[φ1, φ2, ···, φs], H=diag[H1(ω), H1(ω), ···, Hs(ω)]. 显然, 式(13)为计算响应功率谱矩阵的CQC法表达式[4-5]. 换言之, 式(12)与式(13)是等价的.
若仅考察结构部分自由度数目为J的响应, 则式(12)和式(13)将分别改写为
$$ {{\boldsymbol{S}}_{{{{Y}}_r}{{{Y}}_r}}}\left( \omega \right) = {\underline {\boldsymbol{G}} _J}{\left( \omega \right)^*}{{\boldsymbol{S}}_{{{FF}}}}\left( \omega \right){\underline {\boldsymbol{G}} _J}^{\rm{T}} $$ (14) $$ {{\boldsymbol{S}}_{{{{Y}}_r}{{{Y}}_r}}}\left( \omega \right) = {{\boldsymbol{\varPhi}} _J}{{\boldsymbol{H}}^*}{{\boldsymbol{\varPhi}} ^{\rm{T}}}{\boldsymbol{Q}}{{\boldsymbol{S}}_{{{FF}}}}\left( \omega \right){{\boldsymbol{Q}}^{\rm{T}}}{\boldsymbol{\varPhi}} {\boldsymbol{H}}{\boldsymbol{\varPhi}} _J^{\rm{T}} $$ (15) 式中, GJ(ω)为考察自由度对应的Hk(ω)组成的J × m维矩阵, ΦJ为由φk,j (k=1, 2, ···, J; j=1, 2, ···, s)组成的J × s维矩阵.
比较式(14)和式(15)可知, 引入广义频响函数后, 广义频响函数矩阵的维数为J × m维, 远小于CQC法中频响函数矩阵的维数s × s维. 因此, 若可以方便地给出Hk,l(ω), 则式(14)比式(15)更易于实用.
由于式(9), 式(11), 式(12)和式(14)中涉及广义频响函数, 因此本文将基于上述表达式的随机振动分析方法称为广义分析方法. 特别地, 当m=1时式(12)则退化为一致激励问题或单点激励问题.
2. 平稳激励下线性结构响应的辅助简谐激励广义法
广义频响函数的快速、准确确定是广义分析方法的关键. 显然, 对于大型复杂结构, 由式(10)确定广义频响函数进而求解响应功率谱是不方便的. 为此本节给出一种较为便捷的求解方法.
2.1 辅助简谐激励广义法
将F(t)中第l个元素以辅助简谐激励exp(iωt) 代替, 其余量为零, 即构造辅助简谐激励向量Fe,l(ω, t)=[0, ···, 0, exp(iωt) , 0, ···, 0]T. 由式(6)可知, 结构在Fe,l(ω, t)下第k自由度的位移响应Yk,l(ω,t)为
$$ \begin{split} {\underline Y _{k,l}}\left( {\omega ,t} \right) =& \int_{ - \infty }^\infty {{{\underline h }_{k,l}}\left( \tau \right)} \exp \left[ {{\rm{i}}\omega \left( {t - \tau } \right)} \right]{\rm{d}}\tau =\hfill \\ & \left[ {\int_{ - \infty }^\infty {{{\underline h }_{k,l}}\left( \tau \right)} \exp \left( { - {\rm{i}}\omega \tau } \right){\rm{d}}\tau } \right]\exp \left( {{\rm{i}}\omega t} \right) =\hfill \\ & \underline H _{k,l}^{}\left( \omega \right)\exp \left( {{\rm{i}}\omega t} \right) \hfill \end{split} $$ (16) 于是, 式(9)中两个广义频响函数的乘积可表示为
$$ \begin{split} & \underline H _{{k1,p}}^*\left( \omega \right)\underline H _{{k2,q}}^{}\left( \omega \right) =\hfill \\ & \qquad {\left( {\underline H _{{k1,p}}^{}\left( \omega \right) \cdot {{\text{e}}^{{\text{i}}\omega {\text{t}}}}} \right)^*} \cdot \left( {\underline H _{{k2,q}}^{}\left( \omega \right) \cdot {{\text{e}}^{{\text{i}}\omega {\text{t}}}}} \right)= \hfill \\ & \qquad { {{{\underline Y }_{{k1,p}}}\left( {\omega ,t} \right)} ^*} \cdot {{{\underline Y }_{{k2,q}}}\left( {\omega ,t} \right)} \buildrel \Delta \over = {H_{A,{k1},{k2,p,q}}}\left( \omega \right) \hfill \end{split} $$ (17) 相应地, 式(9)可改写为
$$ \begin{split} & {S_{{Y_{{k1}}}{Y_{{k2}}}}}\left( \omega \right) = \sum\limits_{p = 1}^m {\sum\limits_{q = 1}^m {\left[ \left({ {{{ {{{\underline Y }_{{k1,p}}}\left( {\omega ,t} \right)} }^*} \cdot {{{\underline Y }_{{k2,q}}}\left( {\omega ,t} \right)} }}\right) \right.} } \cdot \hfill \\ & \qquad{\text{}}\left. { {S_{{F_p}{F_q}}}\left( \omega \right)} \right] \hfill \end{split} $$ (18) 类似地, 式(11)可改写为
$$ \begin{split} & {S_{{Y_k}{Y_k}}}\left( \omega \right) = \sum\limits_{p = 1}^m {\sum\limits_{q = 1}^m {\left[ {\left( {{{ {{{\underline Y }_{k,p}}\left( {\omega ,t} \right)} }^*} \cdot {{{\underline Y }_{k,q}}\left( {\omega ,t} \right)} } \right)} \right.} } \cdot \hfill \\ & \qquad\left. { {S_{{F_p}{F_q}}}\left( \omega \right)} \right] \hfill \end{split} $$ (19) 式(18)和式(19)亦可表示为如下矩阵形式, 即
$$ {{\boldsymbol{S}}_{{{YY}}}}\left( \omega \right) = \underline {\boldsymbol{Y}} {\left( {\omega ,t} \right)^*}{{\boldsymbol{S}}_{{{FF}}}}\left( \omega \right)\underline {\boldsymbol{Y}} {\left( {\omega ,t} \right)^{\rm{T}}} $$ (20) 式中, Y(ω)=[Y1(ω), Y2(ω), ···, Ys(ω)]T为辅助简谐激励下的响应矩阵, 其中Yk(ω)=[Yk,1(ω), Yk,2(ω), ···, Yk,m(ω)]. 不难发现, 式(20)中SFF(ω)可直接应用, 无需分解.
由于式(18), 式(19)或式(20)综合了广义分析法的思路和辅助简谐激励的策略, 文中将其称为辅助简谐激励广义法.
值得指出的是, 由式(10)可知, Hk,l(ω)除与结构自身特性有关外, 仅取决于激励定位矩阵Q. 因此, 对于Q相同的随机振动重分析, 辅助简谐激励广义法具有天然的优势.
2.2 辅助简谐激励广义法的高效求解方法
线性结构在辅助简谐激励下的响应分析属于谐振响应分析, 常用分析方法包括振型叠加法与时程分析法. 显然, 当只需少数低阶振型即可获得较为精确响应的前提下, 振型叠加法较时程分析法更为方便; 但对于包含较多振型贡献的大型复杂结构, 时程分析方法可以有效避免振型叠加法中振型截断导致的误差. 因此, 文中根据不同的情形给出两种可供选择的辅助简谐激励响应计算方案.
2.2.1 基于振型叠加的辅助简谐激励响应计算
当结构自由度较少或者可以预先确定准确计算辅助简谐激励响应所需的振型最高阶数且此阶数较小时, 可采用振型叠加法计算之, 即
$$ {\underline Y _{k,l}}\left( {\omega ,t} \right) = \sum\limits_{j = 1}^r {{\varphi _{k,j}}{H_j}\left( \omega \right)} {\boldsymbol{\varphi}} _j^{\rm{T}}{{\boldsymbol{Q}}_l}\exp \left( {{\rm{i}}\omega t} \right) $$ (21) 式中, r为自由度数目(即r=s)或需要考虑的振型阶数(即r<s). 需要指出的是, 由于式(17)和式(18)中Yk1,p(ω,t)与Yk2,q(ω,t)共轭相乘导致exp(iωt)相互抵消, 所以式(21)在实际计算中可省略exp(iωt).
2.2.2 基于时程分析的辅助简谐激励响应计算
对于复杂高维结构, 为避免振型叠加法中振型截断的误差, 可采用时程分析法计算辅助简谐激励响应. 由于涉及到不同频率下辅助简谐激励响应的大量重分析, 文中引入时域显式法以提高计算效率.
通常, Yk,l(ω, t)在计算中需要将频率离散化, 记频率离散值的集合为
$$ {\underline \omega _p} \in \left\{ {0,\Delta \omega , 2\Delta \omega ,\cdots ,N\Delta \omega } \right\} \buildrel \Delta \over = {\boldsymbol{W}} $$ (22) 式中, Δω为频率间隔, 即
$$ \Delta \omega = {{\left( {{\omega _{k,U}} - {\omega _{k,L}}} \right)} \mathord{\left/ {\vphantom {{\left( {{\omega _{k,U}} - {\omega _{k,L}}} \right)} N}} \right. } N} $$ (23) 其中, ωk, L为下限截止频率, ωk, U为上限截止频率, (N + 1)为离散频率值数目.
对于W中任一频率离散值对应的简谐激励exp(iωpt) (ωp∈W)均需求解如下的运动微分方程, 即
$$ {\boldsymbol{M}}\underline {\ddot {\boldsymbol{Y}}} \left( {{{\underline \omega }_p},t} \right) + {\boldsymbol{C}}\underline {\dot {\boldsymbol{Y}}} \left( {{{\underline \omega }_p},t} \right) + {\boldsymbol{K}}\underline {\boldsymbol{Y}} \left( {{{\underline \omega }_p},t} \right) = {{\boldsymbol{Q}}_l}{\underline {\boldsymbol{F}} _{e,l}}\left( {{{\underline \omega }_p},t} \right) $$ (24) 其中, Fe,l(ωp, t)=[0, ···, 0, exp(iωpt) , 0, ···, 0]T. 由于ωp∈W, 因此式(24)涉及(N + 1)次结构动力分析, 但对于不同的ωp, 式(24)中M, C, K和Ql并未改变, 因此可采用时域显式法[31-32]进行不同ωp时的高效计算.
根据时域显式法[31-32], 可以给出式(1)中F(t)为任意激励时响应Y(tτ)的显式表达式为
$$ {\boldsymbol{Y}}\left( {{t_\tau }} \right) = {{\boldsymbol{A}}_{\tau ,0}}{\boldsymbol{F}}({t_0}){\text{ }} + {{\boldsymbol{A}}_{\tau ,1}}{\boldsymbol{F}}({t_1}) + \cdots + {{\boldsymbol{A}}_{\tau ,\tau }}{\boldsymbol{F}}({t_\tau }) $$ (25) 式中, [Aτ,0, Aτ,1, ···, Aτ,τ]可通过2m次单位脉冲激励确定, 具体计算方法可参考文献[31-32]. 于是, 式(24)的响应Yk,l(ωp,tτ)可按下式计算, 即
$$ \begin{split} {\underline Y _{k,l}}\left( {{{\underline \omega }_p},{t_\tau }} \right) =& {a_{k,0}} \exp \left( {{\rm{i}}{{\underline \omega }_p}{t_0}} \right) +\\ & {a_{k,1}} \exp \left( {{\rm{i}}{{\underline \omega }_p}{t_1}} \right) + \cdots + {a_{k,\tau }} \exp \left( {{\rm{i}}{{\underline \omega }_p}{t_\tau }} \right) \hfill \end{split} $$ (26) 其中ak,n (n=0, 1, ···,τ)为Aτ, n 中的第k个元素.
显然, 基于时程分析的辅助简谐激励响应计算亦适用于结构自由度较少的情况.
将由式(21)或式(26)得到的辅助简谐激励响应代入到式(20)即可得到响应功率谱, 文中分别称之为基于振型叠加的辅助简谐激励广义法与基于时程分析的辅助简谐激励广义法. 其中, 对于基于振型叠加的辅助简谐激励广义法, 式(20)可改写为
$$ {{\boldsymbol{S}}_{{{YY}}}}\left( \omega \right) = \underline{\underline {\boldsymbol{Y}}} {\left( {\omega ,t} \right)^*}{{\boldsymbol{S}}_{{{FF}}}}\left( \omega \right)\underline{\underline {\boldsymbol{Y}}} {\left( {\omega ,t} \right)^{\rm{T}}} $$ (27) $$ \underline{\underline {\boldsymbol{Y}}} \left( \omega \right) = \underline{\underline {\boldsymbol{\varPhi}} } \;\underline{\underline {\boldsymbol{H}}} \;{\underline{\underline {\boldsymbol{\varPhi}} } ^{\rm{T}}}{\boldsymbol{Q}} $$ (28) 式中,
${\underline{\underline {\boldsymbol{\varPhi}} } } $ 为前r阶振型矩阵,$ {\underline{\underline {\boldsymbol{H}} } } $ = diag[H1(ω), H2(ω), ···, Hr(ω)].2.3 计算步骤
采用辅助简谐激励广义法计算响应功率谱的计算步骤如下:
(1) 通过式(21)或式(26)计算响应Yk,l(ω, t);
(2) 结合式(18)、式(19)和Yk,l(ω, t), 分别计算
${S_{{Y_k}{Y_k}}}\left( \omega \right) $ 和${S_{{Y_{k1}}{Y_{k2}}}}\left( \omega \right) $ .3. 计算性能分析与比较
辅助简谐激励广义法中响应的求解既可采用振型叠加法也可采用时程分析法, 但前者以振型已知为前提, 而已有文献[20, 23]对基于振型叠加的响应功率谱分析方法的性能分析时均不考虑振型计算的工作量, 因此不宜将基于振型叠加的方法与基于时程分析的方法统一对比. 为此, 文中将基于振型叠加的辅助简谐激励广义法与基于振型叠加的其他方法对比, 将基于时程分析的辅助简谐激励广义法与基于时程分析的其他方法对比. 为简化, 仅分析一个频率值处各算法的乘法次数, 且不区分实数和复数.
3.1 基于振型叠加的辅助简谐激励广义法的性能
根据式(27)可知基于振型叠加的辅助简谐激励广义法所需的乘法次数NPPM-1为
$$ {N_{{\rm{PPM}} - 1}}\left( {s,m,r} \right) = 2rsm + {r^2}m + {s^2}m + s{m^2} $$ (29) 其中, 2rsm + r2m为求解广义频响函数矩阵所涉及的乘法次数, s2m + sm2为通过广义频响函数矩阵和激励功率谱矩阵相乘所涉及的乘法次数.
为比较, 本节亦给出了传统CQC法和基于振型叠加的虚拟激励法的计算性能.
传统CQC法的表达式为
$$ {{\boldsymbol{S}}_{{{YY}}}}\left( \omega \right) = \sum\limits_{j = 1}^r {\sum\limits_{k = 1}^r {H_j^*{H_k}{{\boldsymbol{\varphi}} _j}{\boldsymbol{\varphi}} _j^{\rm{T}}} } {\boldsymbol{Q}}{{\boldsymbol{S}}_{{{FF}}}}\left( \omega \right){{\boldsymbol{Q}}^{\rm{T}}}{{\boldsymbol{\varphi}} _k}{\boldsymbol{\varphi}} _k^{\rm{T}} $$ (30) 其所需的总乘法次数NCQC为[20]
$$ {N_{{\rm{CQC}}}}\left( {s,m,r} \right) = {r^2}\left( {2sm + m + {s^2} + {m^2}} \right) $$ (31) 基于振型叠加的虚拟激励法的基本公式为[23]
$$ {{\boldsymbol{S}}_{{{FF}}}}\left( \omega \right) = \sum\limits_{j = 1}^c {{{\boldsymbol{l}}_j}} {\boldsymbol{l}}_j^{\rm{T}} $$ (32) $$ {{\boldsymbol{y}}_j} = \underline{\underline {\boldsymbol{\varPhi}} }\; \underline{\underline {\boldsymbol{H}}}\; {\underline{\underline {\boldsymbol{\varPhi}} } ^{\rm{T}}}{\boldsymbol{Q}}{{\boldsymbol{l}}_j} $$ (33) $$ {{\boldsymbol{S}}_{{{YY}}}}\left( \omega \right) = \sum\limits_{j = 1}^c {{\boldsymbol{y}}_j^*{\boldsymbol{y}}_j^{\rm{T}}} $$ (34) 式中, lj为SFF(ω)分解的m × 1维的虚拟激励向量, yj为lj对应的虚拟响应, c为SFF(ω)的秩, 此处取c = m. 基于振型叠加的虚拟激励法所涉及的乘法次数NPEM-1为
$$ {N_{{\rm{PEM}} - 1}}\left( {s,m,r} \right) = m\left( {2rs + {r^2} + sm + {s^2}} \right) + {{{m^3}}/6} $$ (35) 其中, 式(32)所需要的乘法次数为m3/6, 式(33)按从右到左的顺序进行矩阵相乘所涉及的乘法次数为2rs + r2 + sm + s2, 式(34)需要进行m次响应求解.
对比式(29)、式(31)和式(35), 当s很大时, 传统CQC法、基于振型叠加的虚拟激励法和基于振型叠加的辅助简谐激励广义法的计算量大致为r2s2, ms2和ms2. 显然, 传统CQC的计算量远大于后两者. 结合式(29)和式(35), 基于振型叠加的虚拟激励法和基于振型叠加的辅助简谐激励广义法的计算量差值为激励功率谱矩阵分解所需乘法次数m3/6>0. 显然, 当m≥2时, 基于振型叠加的辅助简谐激励广义法的计算效率高于基于振型叠加的虚拟激励法; 当m=1时, 由于激励功率谱矩阵无需分解, 两者计算量相等. 需指出的是, 若需考虑多个SFF(ω), 式(28)和式(35)可将不包含SFF(ω)的计算结果存储, 而式(35)则需重新分析, 因此式(29)对于此类情形优势将更为突出.
3.2 基于时程分析的辅助简谐激励广义法的性能
根据式(20)与式(26)可知基于时程分析的辅助简谐激励广义法所需的乘法次数NPPM-2为
$$ \begin{split} & {N_{{\rm{PPM}} - 2}}\left( {s,m,p} \right) = m\bigg[ {184{s^3} + {s^2}m\left( {2p - 3} \right)} + \hfill \\ & \qquad {{\text{ }} {{p\left( {p + 3} \right)m\left( {m + 1} \right)s}/4}} \bigg] + {{s^2}m + s{m^2}} \hfill \end{split} $$ (36) 式中, p为时域积分步数.
为比较, 本节亦给出了基于时程分析的虚拟激励法的计算性能. 与基于振型的虚拟激励法相比, 只需将式(33)代替以时程分析结果. 若采用精细积分法进行时程分析, 则虚拟激励法所需乘法次数NPEM-2为[33]
$$ \begin{split} & {N_{PEM - 2}}\left( {s,m,p} \right)= \hfill \\ & \qquad m\left( {184{s^3} + 16{s^2}pm + {s^2}} \right) + {{{m^3}}/6} \hfill \end{split} $$ (37) 图1给出了NPEM-2与NPPM-2的比值和自由度数s的关系图, 其中p=1000. 易知, 当s较小时, 辅助简谐激励广义法的计算效率低于虚拟激励法, 当s增大时, 辅助简谐激励广义法的计算效率将高于虚拟激励法, 且激励数目越多优势越明显.
![]()
图 1 NPEM-2/NPPM-2和自由度数目s的关系Figure 1. The relationship between NPEM-2/NPPM-2 and the number of degrees of freedom s4. 算例分析
算例采用文中建议的辅助简谐激励广义法计算响应功率谱, 且将其与虚拟激励法和传统的CQC法进行对比验证建议方法的精度和效率. 为便于区分, 本节将基于振型叠加的辅助简谐激励广义法记为MAHEGM, 基于时程分析的辅助简谐激励广义法记为TAHEGM, 基于振型叠加的虚拟激励法记为MPEM, 基于时程分析的虚拟激励法记为TPEM, 传统CQC法记为CQC.
4.1 算例1
考察如图2所示的4层剪切框架在随机激励下的随机响应. 每层质量为2.5×104 kg, 刚度为 8.0×106 N/m , 采用 Rayleigh 阻尼模型且阻尼矩阵C =aM + bK, 其中 a = 0.3689, b = 0.003 3. 框架的每一层均作用有随机激励. 根据激励种类的不同分为两种工况, 即工况1 (记为case 1) 为一致激励和工况2 (记为case 2) 为非一致激励.
4.1.1 工况1
该工况下, 向量Q为[1,1,1,1]T. 平稳随机激励F(t)=[F(t)]的功率谱密度函数S(ω)为
$$ \begin{split} & S(\omega ) = \frac{{\omega _g^4 + 4\xi _g^2\omega _g^2{\omega ^2}}}{{ {{{\left( {{\omega ^2} - \omega _g^2} \right)}^2} + 4\xi _g^2\omega _g^2{\omega ^2}} }} \cdot\hfill \\ &\qquad {\text{ }} \frac{{{\omega ^4}}}{{ {{{\left( {{\omega ^2} - \omega _f^2} \right)}^2} + 4\xi _f^2\omega _f^2\omega _{}^2} }} {S_0} \times 6.25 \times {10^8} \hfill \end{split} $$ (38) 其中, ωg=20 rad/s, ωf=2 rad/s, ξg=0.85, ξf=0.85, S0=0.0107m2/s3. 此外, 计算中取ωL=0, Δω =0.1047 rad/s, N =1500.
由于结构自由度数目为4, 本文提出的MAHEGM和TAHEGM都可用于求解响应功率谱. 其中, MAHEGM中取r=4; TAHEGM中, 时域积分区间为[0, 40] s, 积分步数为4000步. 图3和图4分别给出部分楼层位移响应的自功率谱和互功率谱计算结果. 为比较, 图中亦给出了CQC, MPEM[23]和MAHEGM的结果.
通过图3和图4可以看出MAHEGM, TAHEGM, MPEM, TPEM和CQC得到位移响应功率谱曲线在整个频域区间吻合较好, 并根据各方法的位移响应功率谱得到的楼层位移方差也是相等的, 其中1层、2层、3层和4层的位移方差依次为5.46×10−5, 1.89×10−3, 3.41×10−3, 4.44×10−3.
此外, 在TAHEGM中, 根据式(17)可通过Yk,l(ω,t)确定HA,k1,k2,p,q(ω). 图5给出了部分HA,k1,k2,p,q(ω)的计算结果. 与此同时, 图5中亦给出了
$ \underline H _{{k_1},p}^*\left( \omega \right)\underline H _{{k_2},q}^{}\left( \omega \right) $ 的精确解(记为HE,k1,k2,p,q(ω)). 通过对比不难发现, 由TAHEGM计算出的HA,k1,k2,p,q(ω)的幅值与HE,k1,k2,p,q(ω)的幅值在整个频率域上吻合的较好, 验证了本文提出由HA,k1,k2,p,q(ω)代替HE,k1,k2,p,q(ω)的有效性.4.1.2 工况2
该工况下, 向量Q的数值为
$$ {\boldsymbol{Q}} = {\left[ {\begin{array}{*{20}{c}} 1&1&1&0 \\ 1&0&1&1 \end{array}} \right]^{\rm{T}}} $$ (39) 平稳随机激励向量F(t)=[F1(t), F2(t)]T的互功率谱密度函数
${S_{{F_i},{F_j}}}\left( \omega \right) $ 为$$ \begin{split} & {S_{{F_i},{F_j}}}\left( \omega \right) = \exp ( - \omega \left| {i - j} \right|/5) {S_{F,F}}\left( \omega \right) \hfill \\ &\qquad {\text{ }}\left( {i = 1,2;{\text{ }}j = 1,2} \right) \hfill \end{split} $$ (40) 其中, SF,F (ω)为自功率谱密度函数且表达式为
$$ {S_{F,F}}\left( \omega \right) = 6.25 \times {10^8} \times \omega /{\left( {10 + {\omega ^2}} \right)^{4/3}}$$ (41) 此外, 计算中取ωL = 0, Δω =0.1047 rad/s, N =1500. 该工况下仍分别采用MAHEGM和TAHEGM进行功率谱分析, 其中r值、时域积分区间和时间步数与工况1一致.
图6和图7分别给出MAHEGM, TAHEGM, MPEM, TPEM和CQC得到的部分楼层位移的自功率谱和互功率谱. 从图中曲线和表中数值可以看出, MAHEGM和TAHEGM能够得到与MPEM, TPEM和CQC相同精度的结果. 根据各方法的位移响应功率谱曲线得到的楼层位移方差也是相等的, 其中1层、2层、3层和4层的位移方差依次为2.45×10−4, 8.36×10−4, 1.50×10−3, 1.93×10−3.
表1列出工况1和工况2中各方法得到的楼层位移方差时所消耗的时间. 本文所有程序均在处理器为i5-9400 F、主频为2.90 GHz和内存为48 GB的计算机上运行. 通过表1可知, 对于低维结构, 基于振型分解的方法相对时程分析的方法具有明显优势, 且MAHEGM计算效率优势更为显著. 对于一致激励和非一致激励情况下, 基于振型分解的方法中MAHEGM计算效率最高; 基于时程分析的方法中TAHEGM的计算效率最高. 在振型阶数、激励的积分点数和频率离散值数目相同的前提下, 随着激励数目的增多, MAHEGM的计算效率提高越明显. 因此, 对于低自由度系统, 本文建议采用MAHEGM进行随机响应分析.
表 1 不同方法的计算时间Table 1. The computation time of different methodsCase MPEM CQC MAHEGM TPEM TAHEGM 1 0.034 s 1.17 s 0.033 s 80.49 s 12.09 s 2 0.95 s 2.06 s 0.057 s 293.86 s 55.71 s 4.2 算例2
考察如图8(a)所示的三维框架结构在随机激励下的随机响应. 柱截面和梁截面绘制于图8(c), 其中柱截面h=0.35 m, 梁截面h=0.45 m. 材料密度为7.85 × 103 kg/m3, 质量比例阻尼参数和刚度比例阻尼参数分别为0.0778和0.020 3. 该结构在Y方向承受荷载的随机过程, 且荷载作用点的位置如图8(a)和图8(b)中① ~ ㉔. 对于任意i点处与j点处的荷载互功率谱密度函数定义为
$$ {S_{{F_i}{F_j}}}\left( \omega \right) = \sqrt {{S_{{F_i}{F_i}}}\left( \omega \right){S_{{F_j}{F_j}}}\left( \omega \right)} \cdot {\gamma _{ij}}\left( \omega \right) $$ (42) 其中,
${S_{{F_i}{F_i}}}\left( \omega \right)$ 为双边自功率谱密度函数, γij(ω)为第i点和第j点的空间相干函数$${\gamma _{ij}}\left( \omega \right) = \exp \left\{ { - \left| \omega \right|{{\left[ {256{{\left( {{x_i} - {x_j}} \right)}^2} + 100{{\left( {{z_i} - {z_j}} \right)}^2}} \right]}^{0.5}}} \right\} $$ (43) 式中, xi和zi为第i点的X方向和Z方向坐标.
以荷载方差为2.27×105为控制条件, 考虑双边自功率谱密度函数
${S_{{F_i}{F_i}}}\left( \omega \right)$ 的类型对响应的影响, 设计3种工况, 且功率谱形式来源于Davenport谱[34]、Kaimal谱[35]和Karman谱[36], 分别为(1) 工况1(记为case 1)
$$ {S_{{F_i}{F_i}}}\left( \omega \right) = \left( {8.53 \times {{10}^7}} \right)\omega /{\left( {1 + 144{\omega ^2}} \right)^{4/3}} $$ (44) (2) 工况2(记为case 2)
$$ {S_{{F_i}{F_i}}}\left( \omega \right) = \left( {7.967 \times {{10}^7}} \right)/{\left( {1 + 100\omega } \right)^{5{\rm{/}}3}} $$ (45) (3) 工况3(记为case 3)
$${S_{{F_i}{F_i}}}\left( \omega \right) = \left( {5.283 \times {{10}^7}} \right)/{\left( {1 + 700\omega } \right)^{5/6}}$$ (46) 取ωL = 0.01 rad/s, Δω =0.01 rad/s, 同时 N = 100. 上述3种激励的功率谱如图9所示, 由图9曲线对比可知, 尽管3种工况下激励的方差一致, 但3种功率谱曲线在不同频率点处的幅值信息却大不相同.
在工况1中, 采用CQC, TAHEGM和虚拟激励法, 其中, CQC振型阶数取为1, 5和30, TAHEGM的时域积分区间为[0, 20] s, 时域积分步数为2000步, 虚拟激励法采用谐响应分析方法[28](记为SPEM) 并将结果作为标准解. 通过上述方法得到部分楼层位移方差, 如表2所示. 由表可知, CQC在取振型阶数1和5的结果与标准解的误差很大, 这说明主观确定振型阶数对CQC的结果有很大影响. 同时CQC在取振型阶数30的结果与SPEM的结果接近, 说明30阶振型是能够精确计算响应所需阶数, 为此本文将该算例下30阶振型认为是可以预先确定准确计算辅助简谐激励响应所需的振型最高阶数, 并将其作为已知条件, 采用MAHEGM计算位移功率谱并将方差列于表2所示, 显然, MAHEGM-30的结果也能与SPEM的结果接近.
在工况2和工况3中, 采用SPEM, CQC-30, TAHEGM和MAHEGM-30计算位移功率谱并将位移方差列于表2. 表中, CQC-i是CQC通过前i阶振型计算响应功率谱; MAHEGM-30是MAHEGM通过前30阶振型计算响应功率谱; ε为相对标准差. 由表可知, 相对SPEM, CQC-30, MAHEGM-30和TAHEGM在这两种工况中都可得到比较精确的结果.
表 2 不同楼层位移方差Table 2. The variances of displacements at different floorsCase Method Floor 3 5 7 Value/10−9 ε/% Value/10−9 ε/% Value/10−9 ε/% 1 SPEM 3.61 − 7.28 − 8.41 − CQC-1 2.17 39.89 5.01 31.18 6.72 20.01 CQC-5 2.97 17.73 6.82 6.32 9.10 8.20 CQC-30 3.61 0.00 7.29 0.14 8.49 0.95 MAHEGM-30 3.61 0.00 7.29 0.14 8.49 0.95 TAHEGM 3.61 0.00 7.28 0.00 8.41 0.00 2 SPEM 6.10 − 11.90 − 13.50 − CQC-30 6.11 0.16 11.90 0.00 13.60 0.74 MAHEGM -30 6.11 0.16 11.90 0.00 13.60 0.74 TAHEGM 6.10 0.00 11.90 0.00 13.50 0.00 3 SPEM 4.44 − 8.81 − 10.10 − CQC-30 4.44 0.00 8.82 0.11 10.20 0.99 MAHEGM -30 4.44 0.00 8.82 0.11 10.20 0.99 TAHEGM 4.44 0.00 8.81 0.00 10.10 0.00 图10给出工况1中各方法得到的第6层自功率谱. 由图中曲线可知, CQC-1明显偏离SPEM对应的曲线, CQC-5的曲线与SPEM的曲线比较接近, TAHEGM, MAHEGM-30, 和CQC-30和SPEM的曲线比较吻合.
![]()
图 10 在工况1中第6层位移自功率谱SY6Y6 (ω)对数图Figure 10. The auto-power spectrum SY6Y6 (ω) of displacement at 6th floor in case 1 with logarithmic scale图11和图12给出由SPEM, TAHEGM, MAHEGM-30和CQC-30在不同工况下的楼层位移自功率谱和互功率谱的对数坐标图. 由图中曲线可知, MAHEGM-30和TAHEGM与SPEM在整个频域区间吻合较好.
![]()
图 11 第4层位移自功率谱SY4Y4(ω)对数图Figure 11. The auto-power spectrum SY4Y4 (ω) of displacement at 4th floor with logarithmic scale![]()
图 12 楼层位移的互功率谱SY3Y5 (ω)对数图Figure 12. The cross-power spectrum SY3Y5 (ω) of displacement with logarithmic scale除位移响应外, 文中亦对内力响应进行了分析. 由于直接采用振型叠加类方法计算内力功率谱时需传递矩阵[23], 计算过程复杂, 不便于应用, 因此, 文中仅给出了TAHEGM和SPEM的结果, 如图13和图14所示. 从图中可以明显地看出TAHEGM所得曲线与SPEM的对应曲线在整个频域区间保持一致.
![]()
图 13 点1的剪力功率谱对数图Figure 13. The power spectrum of shear force at point 1 with logarithmic scale![]()
图 14 点2的弯矩功率谱对数图Figure 14. The power spectrum of bending moment at point 2 with logarithmic scale表3给出了不同方法在不同工况下计算位移方差的耗时和总耗时. 由表可知, MAHEGM-30的计算效率远高于CQC-30. 虽然MAHEGM-30比TAHEGM有较高的计算效率, 但MAHEGM需要以振型已知为前提, 对大型复杂结构, MAHEGM适用范围将会受限. TAHEGM的耗时明显低于SPEM, 值得指出的是TAHEGM在工况1中计算位移方差所消耗的时间为7579.91 s, 其中在TAHEGM中确定系数所消耗的时间为7543.38 s, 计算简谐激励对应的响应时间为32.52 s, 与激励的功率谱对应相乘的时间仅为0.017 s. 由于TAHEGM在工况1中已经存储HA, k1, k2, p, q(ω), 在工况2或工况3作用中求解时只需将HA, k1, k2, p, q(ω)与激励功率谱矩阵对应相乘相加即可, 不需要对激励功率谱矩阵进行重分析. 因此, 在工况2和工况3中TAHEGM计算每层位移方差所消耗的时间为0.017 s. 因此, TAHEGM在针对高维复杂结构遭受多工况的响应功率谱计算具有显著优势.
比较表2中同一楼层在不同工况下的位移方差数值可以看出在激励的方差相同而激励形式不同时, 位移方差却有明显差异, 相对差异最大达到69%左右, 说明在结构设计中对激励形式的选择需要慎重考虑, 进一步说明有必要对同一结构进行激励参数分析.
通过图10 ~ 图14也可以发现工况2和工况3在低频处对响应功率谱幅值的影响比工况1大, 说明不同的激励形式对响应功率谱的幅值影响也不尽相同.
表 3 不同方法的计算时间Table 3. The computation time of different methodsMethod Time/s Total time/s Case 1 Case 2 Case 3 SPEM 34993.39 34993.39 34993.39 104980.17 CQC-1 2.31 − − 2.31 CQC-5 13.47 − − 13.47 CQC-30 335.22 128.22 128.22 591.66 MAHEGM -30 17.17 0.22 0.22 17.61 TAHEGM 7579.91 0.017 0.017 7579.944 注: − 表示根据计算结果不需要计算. 5. 结论
对于平稳高斯激励下线性结构随机振动问题, 针对现有方法必须依赖振型截断或激励功率谱矩阵分解的不足, 本文首先推导了与CQC法等价的广义响应分析方法, 然后引入辅助的简谐激励实现了其直接计算, 从而形成了一类新的随机振动分析方法—辅助简谐激励广义法. 在此基础上, 结合振型分解法或时域显式法对辅助简谐激励响应的高效求解, 给出了辅助简谐激励广义法的两种高效算法. 算例结果表明, 辅助简谐激励广义法在保证与虚拟激励法计算精度相当的同时, 其计算效率相对虚拟激励法明显提高, 特别是针对同一结构的多工况问题, 辅助简谐激励广义法的整体计算效率显著提高.
值得指出的是, 尽管本文建议方法是以平稳激励下线性结构响应为研究对象, 但类似思路可以拓展至非平稳激励[37]和非线性结构[38-39]的情形, 亦是后续的研究重点.
-
![]()
图 1 NPEM-2/NPPM-2和自由度数目s的关系
Figure 1. The relationship between NPEM-2/NPPM-2 and the number of degrees of freedom s
![]()
图 10 在工况1中第6层位移自功率谱SY6Y6 (ω)对数图
Figure 10. The auto-power spectrum SY6Y6 (ω) of displacement at 6th floor in case 1 with logarithmic scale
![]()
图 11 第4层位移自功率谱SY4Y4(ω)对数图
Figure 11. The auto-power spectrum SY4Y4 (ω) of displacement at 4th floor with logarithmic scale
![]()
图 12 楼层位移的互功率谱SY3Y5 (ω)对数图
Figure 12. The cross-power spectrum SY3Y5 (ω) of displacement with logarithmic scale
![]()
图 13 点1的剪力功率谱对数图
Figure 13. The power spectrum of shear force at point 1 with logarithmic scale
![]()
图 14 点2的弯矩功率谱对数图
Figure 14. The power spectrum of bending moment at point 2 with logarithmic scale
表 1 不同方法的计算时间
Table 1 The computation time of different methods
Case MPEM CQC MAHEGM TPEM TAHEGM 1 0.034 s 1.17 s 0.033 s 80.49 s 12.09 s 2 0.95 s 2.06 s 0.057 s 293.86 s 55.71 s 
下载: 导出CSV
表 2 不同楼层位移方差
Table 2 The variances of displacements at different floors
Case Method Floor 3 5 7 Value/10−9 ε/% Value/10−9 ε/% Value/10−9 ε/% 1 SPEM 3.61 − 7.28 − 8.41 − CQC-1 2.17 39.89 5.01 31.18 6.72 20.01 CQC-5 2.97 17.73 6.82 6.32 9.10 8.20 CQC-30 3.61 0.00 7.29 0.14 8.49 0.95 MAHEGM-30 3.61 0.00 7.29 0.14 8.49 0.95 TAHEGM 3.61 0.00 7.28 0.00 8.41 0.00 2 SPEM 6.10 − 11.90 − 13.50 − CQC-30 6.11 0.16 11.90 0.00 13.60 0.74 MAHEGM -30 6.11 0.16 11.90 0.00 13.60 0.74 TAHEGM 6.10 0.00 11.90 0.00 13.50 0.00 3 SPEM 4.44 − 8.81 − 10.10 − CQC-30 4.44 0.00 8.82 0.11 10.20 0.99 MAHEGM -30 4.44 0.00 8.82 0.11 10.20 0.99 TAHEGM 4.44 0.00 8.81 0.00 10.10 0.00
下载: 导出CSV
表 3 不同方法的计算时间
Table 3 The computation time of different methods
Method Time/s Total time/s Case 1 Case 2 Case 3 SPEM 34993.39 34993.39 34993.39 104980.17 CQC-1 2.31 − − 2.31 CQC-5 13.47 − − 13.47 CQC-30 335.22 128.22 128.22 591.66 MAHEGM -30 17.17 0.22 0.22 17.61 TAHEGM 7579.91 0.017 0.017 7579.944 注: − 表示根据计算结果不需要计算.
下载: 导出CSV
-
[1] Zhong WX. Review of a high-efficiency algorithm series for structural random responses. Progress in Natural Science, 1996, 6(3): 257-268
[2] Soong TT, Grigoriu M. Random Vibration of Mechanical and Structural Systems. Englewoods Cliffs: Prentice Hall, 1993
[3] 欧进萍, 王光远. 结构随机振动. 北京: 高等教育出版社, 1998 Ou Jinping, Wang Guangyuan. Structural Random Vibration. Beijing: Higher Education Press, 1998 (in Chinese)
[4] Li J, Chen JB. Stochastic Dynamics of Structures. Wiley, 2010
[5] Lutes LD, Sarkani S. Random Vibrations: Analysis of Structural and Mechanical Systems. Amsterdam: Elsevier, 2004
[6] 方同. 工程随机振动. 北京: 国防工业出版社, 1995 Fang Tong. Random Vibration of Engineering. Beijing: National Defense Industry Press, 1998 (in Chinese)
[7] Falsone G, Settineri D. New differential equations governing the response cross-correlations of linear systems subjected to coloured loads. Journal of Sound and Vibration, 2011, 330(12): 2910-2927 doi: 10.1016/j.jsv.2010.12.020
[8] Tootkaboni M, Graham-Brady L. Stochastic direct integration schemes for dynamic systems subjected to random excitations. Probabilistic Engineering Mechanics, 2010, 25(2): 163-171
[9] Spanos PD, Zeldin BA. Monte Carlo treatment of random fields: a broad perspective. Applied Mechanics Reviews, 1998, 51(3): 219-237
[10] Rubinstein RY, Kroese DP. Simulation and the Monte Carlo Method (third ed). New Jersey: John Wiley & Sons, 2016
[11] Fuente E. An efficient procedure to obtain exact solutions in random vibration analysis of linear structures. Engineering Structures, 2008, 30(11): 2981-2990 doi: 10.1016/j.engstruct.2008.04.015
[12] 陈科, 於孝朋, 郑红梅等. 车辆悬架未确知动力学模型与动力响应研究. 应用力学学报, 2017, 34(4): 647-653 (Chen Ke, Yu Xiaopeng, Zheng Hongmei, et al. Study on the unascertained-dynamic model and dynamic response of vehicle suspension. Journal of Applied Mechanics, 2017, 34(4): 647-653 (in Chinese) [13] Zhu S, Li Y, Togbenou K, et al. An efficient optimization algorithm to study the stochastic responses of vehicle-bridge coupling systems. Computing in Science & Engineering, 2018, 21(3): 6-17
[14] Chen N, Zhu S, Li Y. The case study of pseudo excitation method combining self-adaptive Gauss integration in random vibration analysis. Shock and Vibration, 2019, 2019(3): 1-11
[15] Zhu S, Li Y. Random characteristics of vehicle-bridge system vibration by an optimized pseudo excitation method. International Journal of Structural Stability and Dynamics, 2020, 20(5): 2050069 doi: 10.1142/S0219455420500698
[16] Rosenblueth E. A basis for aseismic design. [PhD. Thesis]. Urbana: University of Illinois, 1951
[17] Falsone G, Muscolino G. Cross-correlation coefficients and modal combination rules for non-classically damped systems. Earthquake Engineering & Structural Dynamics, 1999, 28(12): 1669-1684
[18] Kiureghian AD, Nakamura Y. CQC modal combination rule for high‐frequency modes. Earthquake Engineering & Structural Dynamics, 2010, 22(11): 943-956
[19] Alkhaleefi AM, Ali A. An efficient multi-point support-motion random vibration analysis technique. Computers & Structures, 2002, 80(22): 1689-1697
[20] 陈凯, 钱基宏. 随机振动问题的广义坐标合成法. 计算力学学报, 2012, 29(2): 171-177 Chen Kai, Qian Jihong. The generalized coordinates synthesis method for the random cibration analysis. Journal of Computational Mechanics, 2012, 29(2): 171-177 (in Chinese)
[21] 陈塑寰, 刘中生, 赵又群. 振型一阶导数的高精度截尾模态展开法. 力学学报, 1993, 25(4): 427-434 (Chen Suhuan, Liu zhongsheng, Zhao Youqun. An accurate modal superposition method for computing mode shape derivatives in structure dynamics. Journal of Theoretical and Applied Mechanics, 1993, 25(4): 427-434 (in Chinese) [22] 谢壮宁. 风致复杂结构随机振动分析的一种快速算法—谐波激励法. 应用力学学报, 2007, 24(2): 263-266 (Xie Zhuangning. New rapid algorithm for wind-induced random vibration of complex structure. Journal of Applied Mechanics, 2007, 24(2): 263-266 (in Chinese) doi: 10.3969/j.issn.1000-4939.2007.02.020 [23] 林家浩, 钟万勰. 关于虚拟激励法与结构随机响应的注记. 计算力学学报, 1998, 15(2): 217-223 (Lin Jiahao, Zhong Wanxie. Some notes on PEM and structural random response analysis. Journal of Computational Mechanics, 1998, 15(2): 217-223 (in Chinese) [24] Lin J, Yan Z, Zhang Y. Accurate and highly efficient algorithms for structural stationary/non-stationary random responses. Computer Methods in Applied Mechanics & Engineering, 2001, 191(1/2): 103-111
[25] 林家浩, 张亚辉. 随机振动的虚拟激励法. 北京: 科学出版社, 2006 Lin Jiahao, Zhang Yahui. Pseudo Excitation Method for Random Vibration. Beijing: Science Press, 2004 (in Chinese)
[26] 林家浩, 张亚辉, 赵岩. 虚拟激励法在国内外工程界的应用回顾与展望. 应用数学和力学, 2017, 38(1): 1-32 (Lin Jiahao, Zhang Yahui, Zhao Yan. The pseudo-excitation method and its industrial applications in China and abroad. Applied Mathematics and Mechanics, 2017, 38(1): 1-32 (in Chinese) doi: 10.1007/s10483-016-2152-6 [27] 王腾飞, 肖绯雄, 贾宏宇等. 基于多维多点虚拟激励法对转向架构架随机振动疲劳寿命预. 振动与冲击, 2020, 39(22): 192-197 (Wang Tengfei, Xiao Feixiong, Jia Hongyu, et al. Fatigue life prediction of bogie frames under random vibration based on the multi-dimensional and multi-support pseudo excitation method. Journal of Vibration and Shock, 2020, 39(22): 192-197 (in Chinese) [28] 赵银庆. 基于虚拟激励的车体结构随机疲劳问题研究. 大连: 大连理工大学, 2010 Zhao Yingqing. Resarch on random fatigue analysis of the vehicle body based on pseudo excitation method. [PhD Thesis]. Dalian: Dalian University of Technology, 2010 (in Chinese)
[29] 王之宏. 风荷载的模拟研究. 建筑结构学报, 1994, 15(1): 44-52 (Wang Zhihong. Simulation of wind loading. Journal of Buliding Structures, 1994, 15(1): 44-52 (in Chinese) doi: 10.3321/j.issn:1000-6869.1994.01.001 [30] R. 克拉夫, J. 彭津. 结构动力学, 第2版. 北京: 高等教育出版社, 2006 Clough R, Penzien J. Structural Dynamics (Version 2). Beijing: Higher Education Press, 2006 (in Chinese)
[31] 苏成, 徐瑞. 非平稳激励下结构随机振动时域分析法. 工程力学, 2010, 42(3): 512-520 (Su Cheng, Xu Rui. Time-domain method for dynamic reliability of structural systems subjected to non-stationary random excitations. Chinese Journal of Theoretical and Applied Mechanics, 2010, 42(3): 512-520 (in Chinese) [32] 苏成, 黄志坚, 刘小璐. 高层建筑地震作用计算的时域显式随机模拟法. 建筑结构学报, 2015, 36(1): 13-22 (Su Cheng, Huang Zhijian, Liu Xiaolu. Time-domain explicit random simulation method for seismic analysis of tall buildings. Journal of Building Structures, 2015, 36(1): 13-22 (in Chinese) [33] 徐瑞, 苏成. 结构非平稳随机响应分析的快速虚拟激励法. 计算力学学报, 2010, 27(5): 822-827 (Xu Rui, Su Cheng. Fast pseudo-excitation method in structural non-stationary stochastic response analysis. Journal of Computational Mechanics, 2010, 27(5): 822-827 doi: 10.7511/jslx20105013 [34] Davenport AG, The spectrum of horizontal gustiness near the ground in high winds. Quarterly Journal of the Royal Meteorological Society, 1961, 88(376): 194-211
[35] Kaimal JC, Wyngaard JC, Izumi Y, et al. Spectral characteristics of surface-layer turbulence. Quarterly Journal of the Royal Meteorological Society, 1972, 98(417): 563-589 doi: 10.1002/qj.49709841707
[36] 黄本才. 结构抗风分析原理及应用. 上海: 同济大学出版社, 2001 Huang Bencai. Principle and Application of Structural Wind Resistance Analysis. Shanghai: Tongji University Press, 2001 (in Chinese)
[37] 芮珍梅, 陈建兵. 加性非平稳激励下结构速度响应的FPK方程降维. 力学学报, 2019, 51(3): 922-931 (Rui Zhenmei, Chen Jianbing. Dimension reduction of FPK equation for velocity response analysis of structures subjected to additive nonstationary excitations. Journal of Theoretical and Applied Mechanics, 2019, 51(3): 922-931 (in Chinese) doi: 10.6052/0459-1879-18-411 [38] Liu RN, Kang YM, Fu YX. Stochastic resonance and bifurcation of order parameter in a coupled system of underdamped Duffing oscillators. International Journal of Bifurcation and Chaos, 2019, 29(8): 1950108
[39] Liu RN, Kang YM. Observing stochastic resonance in periodic potential systems driven by alpha stable Lévy noise. Physics Letter A, 2018, 382: 1656-1664
计量
- 文章访问数: 885
- HTML全文浏览量: 364
- PDF下载量: 91
