非均质结构非平稳随机响应的快速算法

非均质结构非平稳随机响应的快速算法

陈玉震, 张盛, 陈飙松, 张洪武

(大连理工大学工业装备结构分析国家重点实验室运载工程与力学学部工程力学系,辽宁大连116024)

摘要：将扩展多尺度有限元法应用于非平稳随机振动时域显式法中，实现了对非均质结构非平稳随机响应的快速精确计算。首先，论文阐述了扩展多尺度有限元法基本原理。其次，探讨了时域显式法在非平稳随机响应分析中的优势。时域显式法基于动力响应显式表达式直接在时域进行随机振动分析，较传统反应谱法，具有良好的计算精度和计算效率。最后，根据两算法的特性和优势提出统一的多尺度算法框架，使其适用于非均质结构非平稳随机响应分析的快速求解。数值算例验证了该算法的高效性和精确性。

关键词：非均质结构; 非平稳随机振动； 扩展多尺度有限元法; 基函数； 时域显式法

中图分类号:TB122

文献标志码：A

DOI:10.13465/j.cnki.jvs.2015.19.004

Abstract:A fast computation was implemented for non-stationary random responses of heterogeneous material structures by applying the extended multiscale finite element method (EMsFEM) into the time-domain explicit method(TDEM). Firstly, the fundamental principle of EMsFEM was presented. Then, the advantages of TDEM in non-stationary random response analysis were explored. Based on the explicit expressions of dynamic responses, the random response analysis was performed directly in time domain, it was shown that TDEM has good accuracy and efficiency compared with the response spectrum method. Finally, according to the characteristics and advantages of the two algorithms, a unified multiscale framework was proposed, it was applicable for non-stationary random response analysis of heterogeneous material structures. Numerical examples showed that the proposed method has high efficiency and accuracy.

基金项目：国家自然科学基金资助项目(51275540)；中央高校基本科研业务费资助项目(CDJZR13110001)

收稿日期：2014-08-07修改稿收到日期：2014-10-17

A fast algorithm for non-stationary random responses of heterogeneous material structures

CHENYu-zhen,ZHANGSheng,CHENBiao-song,ZHANGHong-wu(State Key Laboratory of Structural Analysis for Industrial Equipment, Department of Engineering Mechanics, Faculty of Vehicle Engineering and Mechanics, Dalian University of Technology, Dalian 116024, China)

Key words:heterogeneous material structures; non-stationary random response; extended multiscale finite element method (EMsFEM); base function; time-domain explicit method (TDEM)

自然存在以及人工形成的大部分材料都具有非均质特性，例如地下岩土以及航空航天工业中广泛使用的复合材料等，由这些非均质材料构成的结构，在工程实践中每时每刻都要受到各种形式的动态载荷，其中，非平稳随机激励是最为普遍的一种载荷形式，例如地震作用、阵风载荷等，为了确保这些结构的可靠性与安全性，研究非均质结构在非平稳随机激励作用下的动力响应具有非常重要的意义。

在随机振动分析中，功率谱法是求解随机响应的主要方法[1-2]，对于多点(m点)非平稳随机响应问题，为了求解各离散频点处的响应功率谱，功率谱法在每个频点都需要进行约m次时程分析，对于大规模结构，计算量巨大，工程中难以接受。时域显式法[3-4](TDEM)从时域出发，根据线性动力系统输入与输出之间的线性关系，建立各时刻结构响应关于随机激励的显式表达式，从而可以直接利用一阶矩和二阶矩的运算规律计算各时刻结构响应的均值和方差，对于在m个非平稳随机激励作用下的结构，显式表达式的建立只需要进行2m次时程分析，计算量小，工程中容易接受。时域显式法中结构响应已经完全解耦，只需要对所关心的响应求解即可，实现了空间尺度降维。当随机响应的均值和方差为时间的慢变函数时，时域显式法可以在较大时间间隔上计算响应的均值和方差，实现了时间尺度降维。时域显式法仅需要非平稳随机激励时域相关函数信息，所以可以摒弃各种功率谱模型，适用范围更广。时域显式法是一种高效、精确的非平稳随机振动分析方法。

时域显式法在求解非均质结构的非平稳随机响应时，为了捕获结构的微观非均质特性，结构网格需要剖分的非常细致，这使得结构规模巨大，甚至无法计算。多尺度方法[5]是求解非均质结构力学性能的一种有效途径，本文选用扩展多尺度有限元法[6](EMsFEM)，且为了能够求解结构动力响应，与传统扩展多尺度有限元法不同，该方法[7]是以各时间步内等效静力平衡方程为基础，数值构造出满足局部特性微分算子并能反映结构动力性能的多尺度基函数，从而在宏观尺度上对原问题进行求解，很大程度地减少计算量。时域显式法的主要计算来自2m次时程分析，将扩展多尺度有限元法引入时域显式法中，能够很大程度的减少每次时程分析的计算量，从而使得原问题能够快速求解。数值算例表明，扩展多尺度有限法与时域显式法的结合应用，在保证精度的同时大大提高了计算效率，为非均质结构的非平稳随机响应分析提供了一种快速精确的计算方法。

1扩展多尺度有限元法

扩展多尺度有限元法的基本思想是通过局部求解宏观单元域内的平衡方程，数值构造出宏观单元的多尺度基函数(形函数)，这些基函数能反映宏观单元内部的微观非均质特征，从而只需在宏观尺度上对原问题进行求解即可以获得满意的结果，大大减少了计算量。扩展多尺度有限元法的计算流程主要可分为三个部分[8-9]，即多尺度基函数构造、宏观计算、降尺度计算，本节以二维问题为例简要介绍一下这一流程。

图1　四节点宏观单元 Fig.1 Four-node coarse element

1.1多尺度基函数构造

见图1的某平面四节点宏观单元，其内部微观节点的位移可以表示为

(1)

式中:N为宏观单元的多尺度基函数矩阵；u单胞内所有微观节点的位移向量；u′E为宏观单元四个节点的位移向量。它们可以表达成如下公式

u=[u1v1u2v2…unvn]T

N=[RT1RT2…RTn]T

(2)

其中

(3)

式中,n为子网格的微观节点总数，Nixy表示在宏观节点i发生x方向单位位移时，宏观单元的子网格上所有微观节点产生的y向位移值。

1.2宏观计算

当多尺度基函数构造完成后，即可推导得到相应宏观单元的等效刚度阵。见图1，考虑子网格内任意一个细网格单元e，其节点编号为(e1，e2，e3，e4)，单元e的弹性应变能Πe可表达为

(4)

式中:Ke是单元e的常规刚度矩阵；ue是单元e的节点位移向量。

由式(1)和式(2)可得细网格单元e的节点位移向量与相应的宏观单元节点位移向量的关系为

ue=Geu′E,Ge=[RTe1RTe2RTe3RTe4]T

(5)

式中:Ge为细网格单元e的转换矩阵，其表征细网格单元e和对应的宏观单元之间的节点位移映射关系。

将宏观单元内部的所有细网格单元应变能力相加，可得到该宏观单元的等效刚度阵

(6)

式中:p为宏观单元内部细网格单元的总数。

宏观单元的等效刚度阵获取后，就可以组装总体刚度阵，建立宏观层次的求解方程，进而进行宏观尺度的求解。

1.3降尺度计算

宏观计算后，即可利用式(1)得到宏观单元内部任一细网格单元e的节点位移ue，进而可以得到该单元的应力应变

εe=Beue=Teu′E,Te=BeGe

(7)

和

σe=Deεe=Seu′E,Se=DeBeGe

(8)

依次在各宏观单元内计算，即可得到整个结构的微观位移Yk+1、微观应力σk+1和微观应变εk+1。

2非平稳随机振动时域显式法

2.1动力响应显式表达式

对于n个自由度的线性体系，在m个激励作用下结构体系的运动方程可写为

(9)

文献[3-4]从精细积分法出发推导动力响应显式表达式，本文从更一般的时程积分法Newmark法出发推导动力响应显式表达式。记时间步长为Δt，时间步数为l，外载F(t)在时刻t0,t1,…,tl处的值分别记做F0,F1,…,Fl。结构第k+1步的响应可由第k步的响应递推得到，递推公式为

Vk+1=TVk+WFk+1

(10)

式中,

(11)

其中

(12)

(13)

(15)

记做

V0=QF0

(16)

则由递推公式可写出时刻ti=iΔt时刻，结构动力响应Vi的表达式

Vi=Ai,0F0+Ai,1F1+…+Ai,iFi

(i=1,2,…,l)

(17)

式中,Ai,0,Ai,1,…Ai,i的值为

(18)

各时刻响应递推公式展开后形式如下

V1=A1,0F0+A1,1F1

V2=A2,0F0+A2,1F1+A2,2F2

⋮

Vl=Al,0F0+Al,1F1+…+

Al,l-1Fl-1+Al,lFl

(19)

矩阵形式为

V=AR

(20)

式中:V=[VT1VT2…VTl]T，R=[FT0FT1…FTl]T。

第i时刻的响应，可进一步表示为

Vi=AiRi,(i=1,2,…,l)

(21)

其中

(22)

式(17)或式(21)即为结构动力响应显式表达式。

2.2响应均值和方差

结构动力响应显式表达式只与结构相关，与外载性质无关，当外载F(t)为随机激励时，动力响应显式表达式依然成立。此时，随机激励F(t)在各时刻的离散量F0,F0,…,Fl，分别为一随机激励。根据一阶矩与二阶矩的运算规律，由式(21)可得时刻ti处响应Vi的均值向量和协方差矩阵分别为

E[Vi]=AiE[Ri](i=1,2,…,l)

(23)

cov[Vi,Vi]=Aicov[Ri,Ri]ATi

(i=1,2,…,l)

(24)

式中,E[Ri]和cov[Ri,Ri]分别为随机激励向量的均值向量和协方差矩阵，可分别由F(t)的均值函数向量和互相关函数矩阵求出。

式(21)中，结构响应已完全解耦，可只计算所关心的响应量，即

vi=aiRi(i=1,2,…,l)

(25)

在结构随机振动分析中，一般只对结构某些关键部位的响应感兴趣，并不需要求解所有响应的均值和方差。由式(25)即可实现对关键部位响应均值和方差的单独计算，如下式

μvi=E(vi)=aiE(Ri)(i=1,2,…,l)

(26)

(i=1,2,…,l)

(27)

这样实现了对结构的空间尺度降维，大大节省了计算量，提高了计算效率。

为了满足积分精度和激励描述的需要，时间步长为Δt不能太大，因此动力响应显式表达式的推导是在一系列间隔较小的离散时间点上进行的。但当随机响应的均值和方差为时间t的慢变函数时，并不需要在上述每一个时刻计算响应的均值和方差，即可在较大时间间隔上应用式(26)和(27)，这样实现了对结构的时间尺度降维，进一步提高了计算效率。

由上述过程可以看出，时域显式法求解响应均值和方差时，直接采用随机激励的时域相关函数信息，无需用到在时频域中表征随机激励的时变功率谱密度函数信息，可以摒弃非平稳随机激励的各种功率谱模型，适用范围更广。

2.3时域显式法系数矩阵的快速计算

时域显式法在随机响应分析时主要包含两部分计算：一是动力响应显式表达式中系数矩阵的计算；二是基于动力响应显式表达式的随机响应计算。当只求解结构某些关键部位的响应时，时域显式法的计算量主要来自前者。而由式(18)可以发现，整个系数矩阵除第一列Ai,0与第二列Ai,1(i=1,2,…,l)需要单独计算外，其他各列均可由前两列递推得到。

当结构受单点激励作用(例如随机地震作用)时，Ai,0和Ai,1为3n维向量。由式(17)及式(19)可知，Ai,0的值为在t=0时刻施加单位脉冲载荷F0=1，其他时刻载荷为0时，结构在各时刻的响应Vi。而Ai,1的值为在t=Δt时刻施加单位脉冲载荷F1=1，其他时刻载荷为0时，结构在各时刻的响应Vi。此时，仅需要2次时程分析即可得到系数矩阵前两列，进而递推得到整个系数矩阵。

当结构受m(m≥2)点激励作用时，Ai,0和Ai,1为3n×m维矩阵。此时Ai,0每列的值分别为t=0时刻在各激励点源依次施加单位脉冲载荷F0,j=1(j表示第j个激励点源，j=1,2,…,m)，其他时刻载荷为0时，结构在各时刻的响应Vi,j。而Ai,1每列的值分别为t=Δt时刻在各激励点源依次施加单位脉冲载荷F1,j=1(j=1,2,…,m)，其他时刻载荷为0时，结构在各时刻的响应Vi,j。此时，仅需要2m次时程分析即可得到系数矩阵前两列，进而递推得到整个系数矩阵。

3时域显式法的多尺度框架

由第2节内容可知，采用时域显式法计算非平稳随机振动的主要计算量来自2m次特定脉冲激励下的时程分析，提高单次时程分析的计算效率成为提高整个算法效率的关键。对于非均质结构，采用多尺度方法是减小问题规模，提高计算效率的有效途径，为了能够求解非均质结构的动力响应，本文所采用的扩展多尺度方法以各时间步内等效静力平衡方程为基础来构造多尺度基函数，其余过程与第1节相同。为了提高计算精度，本文采用多节点宏观单元策略[10]。本节将详细介绍一下将多尺度方法引入时域显式法，构建时域显式法多尺度框架的基本流程。框架图见图2，基本流程如下:

(1)读取粗细网格信息，基于各时间步等效静力平衡方程计算各宏观单元动力基函数Nd。

(4)构造2m个脉冲激励，1:m个激励在t=0时刻施加，m+1:2m个激励在t=Δt时刻施加，其幅值均为：F=JE，J为一n×m阶常数矩阵，用于定位激励，E为m×m阶单位矩阵。

(5)初始状态计算，根据构造激励计算宏观初始状态向量，并降尺度得到微观初始状态向量。

(6)各时间步内等效载荷为

(28)

此时，载荷直接施加在微观网格上，这将产生不平衡节点力，本文采用位移分解技术[10]来处理，真实位移u被分解成宏观整体响应u′和局部扰动响应u″，u′由宏观等效力驱动，u″则由微观扰动力驱动，本步骤即计算第k步宏观整体等效载荷及各单胞域内微观等效载荷。

(9)判断是否达到最大时间步，完成整个时程分析计算。

(10)提取2m个时程响应结果，组装得到系数矩阵前两列Ai,0与Ai,1(i=1,2,…,l)。

(11)由系数矩阵前两列，递推得到整个系数矩阵A。

(12)通过式(26)或(27)计算响应均值和方差。

图2　时域显式法多尺度框架 Fig.2 The multiscale framework for TDEM

4数值算例

通过2个数值算例来验证本文算法的有效性。为了便于比较，两算例均以虚拟激励法的计算结果作为参考解。

两算例均采用下列形式的均匀调制非平稳激励

F(t)=g(t)f(t)

(29)

式中,g(t)为均匀调制函数，有两种形式

g1(t)=1

(30)

(31)

f(t)的自谱密度均取过滤白噪声模型

(32)

式中,tb=8.5s,tc=20.0s,c1=0.1572,S0=142.75,ωg=19.07,ζg=0.544。

两算例中，结构尺寸及材料属性等数据均为无量纲量。计算所取时间和频率分别为t∈[0,40],Δt=0.5和ω∈[0,50],Δω=0.5。均采用瑞利阻尼，系数为α=0.015和β=0.02。

4.1算例1周期性非均质结构受非平稳水平地面加速度响应分析

见图3(a)T型结构，无量纲尺寸图中已给出，细网格单元尺寸为10×10，粗网格单元尺寸为100×100，粗网格剖分及单胞见图3(b)。结构底部全约束，整体受一非平稳水平地面加速度作用，均匀调制函数由式(31)给定，激励自谱密度由式(32)给定。结构由5种无量纲材料周期性构成，各材料属性见表1，计算右上顶点处的位移均方值。

图3　T型结构模型图 Fig.3 The T-shape structure model

材料密度弹性模量泊松比E150003.0e110.4E230001.0e110.3E338002.0e110.35E432002.1e110.32E548002.5e110.38

粗网格单元边界宏观节点数依次取2、4、6和8，四种情况下虚拟激励法(PEM)、时域显式法(TDEM)以及扩展多尺度时域显式法(EMsTDEMm)的计算结果见图4，其中m表示粗网格单元边界宏观节点数目。以虚拟激励法(PEM)作为数值参考解，其他方法相对数值误差通过L2范数定义，其公式如下：

(33)

式中,va表示PEM参考解，v表示TDEM及EMsTDEMm计算结果；n为时间步数。

图4　端点位移均方值 Fig.4 Mean square of the endpoint displacement

由图4及表2中的L2范数可以看出，TDEM与PEM的计算精度一致。随着边界宏观节点数目的增加，EMsTDEMm的计算精度逐渐趋近PEM，当m=8时，基本一致。

由表2中的计算时间可得，TDEM计算用时远远小于PEM。随着边界宏观节点数目的增加，EMsTDEMm的计算用时逐渐增加，但效率较TDEM依然有极大的提高。

表2　L 2范数及计算用时

4.2算例2 随机非均质悬臂梁结构受多点部分相干随机激励响应分析

见图5(a)T型结构，无量纲尺寸图中已给出，细网格单元尺寸为5×2.5，粗网格单元尺寸为50×25，粗网格剖分及单胞见图5(b)。悬臂梁左端全约束，上边界中点与右端点处施加两部分相干随机激励。结构由3种材料构成(见表3)，且3种材料在结构中完全随机分布，计算右下端点处的位移均方值。

图5　悬臂梁结构模型图 Fig.5 The cantilever structure model

材料密度弹性模量泊松比E130001.0e110.31E240002.0e110.35E350003.0e110.38

两随机激励的均匀调制函数分别由式(31)和(30)给定，激励自谱密度分别取式(32)及S0，两激励相干矩阵为

(34)

粗网格单元边界宏观节点数依次取2、5、8和11，四种情况下虚拟激励法(PEM)、时域显式法(TDEM)以及扩展多尺度时域显式法(EMsTDEMm)的计算结果见图6，其中m表示粗网格单元边界宏观节点数目，L2范数及计算用时见表4。

图6　结构端点位移方差 Fig.6 The variance of the structure’s endpoint

方法自由度L2范数计算时间/sPEM8282—1325.00TDEM82820.000049.67EMsTDEM21100.14025.06EMsTDEM56740.05896.24EMsTDEM812380.01838.48EMsTDEM1118022.39e-712.80

由图6及表4中的L2范数可以看出，TDEM与PEM的计算精度一致。随着边界宏观节点数目的增加，EMsTDEMm的计算精度逐渐趋近PEM，当m=11时，基本一致。

表4中的计算时间可得，TDEM计算用时远远小于PEM。随着边界宏观节点数目的增加，EMsTDEMm的计算用时逐渐增加，但效率较TDEM依然有极大的提高。

此外，本文提出的扩展多尺度时域显式法(EMsTDEMm)在计算单点、多点完全相干或者激励点源不多的部分相干非平稳随机振动时，具有显著优势，计算效率很高；但在处理激励点源较多(m个)的部分相干非平稳随机振动时，需要计算2m个脉冲激励下结构的时程分析，计算量巨大，该方法会受到一定的限制，因此，需要进一步从并行程序设计角度来解决此问题。

5结论

大型复杂结构在非平稳随机激励作用下的响应分析仍是目前面临的一个难点问题，当结构材料具有非均质特性时，分析变得更加困难。针对这一问题，本文从非平稳随机响应分析算法和非均质结构动力响应算法两个方面着手，提出一种适用于非均质结构平稳随机响应分析的快速算法。非平稳随机响应分析算法方面采用时域显式法，该方法从时域出发建立结构动力响应显式表达式，直接通过一阶矩和二阶矩计算结构响应的均值和方差。对于单点非平稳随机响应分析，主要计算仅为2次时程分析，较功率谱法几十至几百次的时程分析，具有较高的计算效率。非均质结构动力响应算法方面，采用扩展多尺度有限元法，该法以时程分析各时间步的等效静力平衡方程为基础，构造能够反映结构动力效应的广义多尺度基函数，从而在宏观尺度上对原动力问题进行求解，较常规有限元法，大大减小了问题规模，提高了计算效率。本文通过采用两种方法相互耦合计算，为非均质复杂结构的非平稳随机响应分析提供了一条可行途径。数值算例表明，所提出算法在保证计算精度的前提下，具有很高的计算效率。

参考文献

[1]Lin J H, Zhang W S, Li J J. Structural responses to arbitrarily coherent stationary random excitations[J]. Computers & structures, 1994, 50: 629-633.

[2]Lin J H, Li J J, Zhang W S, et al. Random seismic responses of multi-support structures in evolutionary inhomogeneous random fields[J]. Earthquake Engineering & Structural Dynamics, 1997, 26: 135-145.

[3]苏成, 徐瑞. 非平稳激励下结构随机振动时域分析法[J]. 工程力学, 2010, 12:77-83.

SU Cheng, XU Rui. Random vibration analysis of structures subjected to non-stationary excitations by time domain method[J]. Engineering Mechanics, 2010,12:77-83.

[4]苏成, 徐瑞, 刘小璐,等. 大跨度空间结构抗震分析的非平稳随机振动时域显式法[J]. 建筑结构学报, 2011, 11:169-176.

SU Cheng, XU Rui, LIU Xiao-lu,et al. Non-stationary seismic analysis of large-span spatial structures by time-domain explicit method[J]. Journal of Building Structures, 2011,11:169-176.

[5]Babuska I, Caloz G, Osborn E. Special finite element methods for a class of second order elliptic problems with rough coefficients[J]. SIAM J. Numer. Anal, 1994,31(4): 945-981.

[6]Zhang H W, Fu Z D, Wu J K. Coupling multiscale finite element method for consolidation analysis of heterogeneous saturated porous media[J]. Advances in Water Resources, 2009, 32(2): 268-279.

[7]Zhang H W, Lu M K, Zheng Y G, et al. General coupling extended multiscale FEM for elasto-plastic consolidation analysis of heterogeneous saturated porous media[J].International Journal for Numerical and Analytical Methods in Geomechanics,2015,39(1):63-95.

[8]张洪武, 吴敬凯, 刘辉,等. 扩展的多尺度有限元法基本原理[J]. 计算机辅助工程, 2010, 19(2): 3-9.

ZHANG Hong-wu, WU Jing-kai, LIU Hui, et al. Basic theory of extended multiscale finite element method[J]. Computer Aided Engineering, 2010, 19(2): 3-9.

[9]张洪武, 吴敬凯, 付振东. 周期性点阵桁架材料力学性能分析的一种新的多尺度计算方法[J]. 固体力学学报, 2011, 32(2): 109-118.

ZHANG Hong-wu, WU Jing-kai, FU Zhen-dong. A new multiscale computational method for mechanical analysis of periodic truss materials[J]. Chinese Journal of Solid Mechanics, 2011, 32(2): 109-118.

[10]Zhang S, Yang D S, Zhang H W, et al. Coupling extended multiscale finite element method for thermoelastic analysis of heterogeneous multiphase materials[J]. Computers and Structures, 2013, 121: 32-49.

第一作者杨洋女，硕士，1988年生

通信作者褚志刚男，博士，副教授，硕士生导师，1978年生