结构动力响应计算的精细积分法_55346
在合理的积分步长范围内，它是不会发生稳定性问题的。即使对高
度病态的问题，精细积分法仍具有很高的精度。
2
将状态方程： z(t) Hz(t) F(t) 的通解：
z(t) eH(tt0 )z(t0 ) eHt
t eHsF(s)dst0源自改写为齐次通解zh与特解z p 之和:
z(t) zh (t) z p (t)
假定在每一积分步长 t [tk ,tk1] 内荷载变化是线性的
即： F(t) r1 (t tk )r2
z(t) Hz(t) F(t)
应用待定系数法可求得特解： z p (t) (H1 tI)(H1r2 ) H1(r1 tkr2 ) z p (t) 代入 z(t) T( t)[z(tk ) z p (tk )] z p (t) 中 得：z(tk1) T(t)[z(t k ) H1(r1 H1r2 )] H1(r1 H1r2 tr2 )
F(t)
Fe(t ) m
1
方程 z(t) Hz(t) F(t) 为一阶线性微分方程组，通解为：
z(t) eH(tt0 )z(t0 ) eHt
t eHsF(s)ds
t0
二、精细积分法
钟万勰在2N 类算法计算指数矩阵的基础上，提出了精细积分法，
用于求解热传导、扩散对流问题以及结构动力学问题的暂态过程。
m
N一般取20，于是 m=1,048,576 106 t
4
t tk1 tk
T(t)=eHt (eH )m
令： t , m 2N
m
I
H
(H )2
2!
(H )3
3!
(H )4
4!
m
(I Ta0 )m
令 I Tai (I Ta,i1)2 I 2Ta,i1 T T a,i1 a,i1 (i 1, 2,...N )
结构动力响应计算的精细积分法 （兼谈 matlab 的学习）
一、状态空间(State space)
单自由度运动方程写为： x(t) 2x(t) 2x(t) Fe (t)
m
设状态向量：
z(t)
x(t)
x(t
)
运动方程可改写为： z(t) Hz(t) F(t)
其中
0 1
H
2
2
0
使用(2)式进行递推运算时排除了 I 参与加法运算的情况， 避免了 Tai 中的元素大数相减而严重丧失有效数字，从而 保证了各 Tai 以及 T(t ) 的高度精确。
这是 精细积分的关键之处 – 创新点!
6
T(t)=(eH )m
I
H
(H )2
2!
(H )3
3!
(H )4
4!
m
(I Ta0 )m
若荷载按简谐变化，则上式在计算机上给出精确的数值解。
9
3.指数衰减型简谐荷载精细积分格式 (HPD-E:High Precision Direct integration scheme-Exponential form):
假定在每一积分步长 t [tk ,tk1] 内荷载按下式变化：
F(t) et (r1 sint r2 cost) (a) z p (t) et (C1 sin t C2 cost)
可求得特解：C1 (2I H2 )1(r2 Hr1) C2 (2I H2 )1(r1 Hr2 )
z p (t) 代入 z(t) T( t)[z(tk ) z p (tk )] z p (t) 中
得：
z(tk1) T(t)[z(t k ) C1 sintk C2 costk ] C1 sintk1 C2 costk1
则 I TaN (I Ta,N1)2 (I Ta,N2 )4 ... (I Ta,0 )m T(t)
5
由递推公式 I Tai I 2Ta,i1 Ta,i1Ta,i1 (i 1, 2,...N ) ---(1) 可得 Tai 2Ta,i1 T T a,i1 a,i1 (i 1, 2,...N ) ---(2)
3
所以： z(t) zh (t) z p (t)
T( t)=eHt ( t t tk ) zh (t) T( t)c c z(tk ) z p (tk )
z(t) T( t)[z(tk ) z p (tk )] z p (t) 令 t tk1 就得到了积分步长终点处的状态
对给定时间步长 t tk1 tk 引入微小时段 令： t , m 2N
在某一积分步 t [tk ,tk1]中，其齐次通解为：
zh (t) T( t)c，其中T( t)=eHt ( t t tk )
当t tk 时：
z(tk ) zh (tk ) z p (tk ) T( t)c z p (tk )
注：当t tk 时
T( t)=eHt e0 I
所以：c z(tk ) z p (tk )
若荷载严格按线性变化，则上式在计算机上给出精确的数值解。
8
2.简谐荷载精细积分格式 (HPD-S:High Precision Direct integration scheme-Sinusoidal form):
假定在每一积分步长 t [tk ,tk1] 内荷载变化是简谐的
即：F(t) r1 sint r2 cost z p (t) C1 sin t C2 cost
上述展开式中略去了 5以及更高阶的项，其量级为： O( 5 ) 1030O(t5 )，已是一般计算机的舍入误差范围
之内。对于线性等步长时域积分而言，T(t )只需计算 一次，计算量是很小的。
7
1.线性荷载精细积分格式 (HPD-L:High Precision Direct integration scheme-Linear form):
可求得特解： C1 [(I H)2 2I]1[(I H)r1 r2 ] C2 [(I H)2 2I]1[(I H)r2 r1]
z p (t) 代入 z(t) T( t)[z(tk ) z p (tk )] z p (t) 中
得： z(tk1) T(t)[z(tk ) eatk (C1 sin tk C2 costk )] etk1 (C1 sin tk1 C2 costk1)
力学中的结构动力学响应与优化
力学中的结构动力学响应与优化力学是研究物体静态和动态力学性质的学科，而结构动力学响应与优化则是力学中的一个重要分支，通过分析结构体在外部力作用下的波动响应，找到最优的结构设计方案。
一、结构动力学响应在力学中，结构动力学响应是指结构体在受到外部力作用后所产生的振动与变形情况。
结构动力学响应可以分为静力响应和动力响应两种情况。
1. 静力响应静力响应是指结构体在受到稳定作用力后的平衡状态。
通过分析材料的力学性质和结构体的几何形状，可以计算出结构体在受力状态下的内力和变形情况。
静力响应的分析方法通常采用力平衡方程和材料本构关系进行计算。
2. 动力响应动力响应是指结构体在受到动态作用力或振动载荷时的响应情况。
动力响应的分析需要考虑结构的惯性和阻尼特性。
通过求解结构的振动方程，可以得到结构体在不同频率下的振动模态和共振情况。
动力响应的分析方法通常采用有限元法、模态分析等数值计算方法。
二、结构动力学优化结构动力学优化是在给定一定的约束条件下，通过调整结构体的形状、材料和结构参数，使得结构体在外部力作用下具有更好的响应性能。
结构动力学优化可以分为静力优化和动力优化两种情况。
1. 静力优化静力优化是指通过调整结构体的形状和几何参数，以使结构体在受力状态下具有更小的应力和变形。
静力优化的目标可以是最小化结构的重量、最大化结构的刚度或满足特定的结构性能要求。
静力优化的方法有拓扑优化、形状优化和尺寸优化等。
2. 动力优化动力优化是指通过调整结构体的参数和材料特性，以使结构体在受到动态作用力或振动载荷时具有更好的阻尼特性和振动响应控制能力。
动力优化的目标可以是最小化结构的振动幅值、最大化结构的振动模态频率或实现特定的振动控制要求。
动力优化的方法有结构参数优化、材料优化和阻尼控制优化等。
结构动力学响应与优化在工程领域具有广泛的应用。
例如，在建筑工程中，通过分析房屋结构在地震作用下的动力响应，可以设计出具有良好抗震性能的建筑物；在航空航天工程中，通过优化飞机结构的动力响应特性，可以提高飞机的飞行稳定性和安全性。
计算结构动力响应的分段精细时程积分方法
计算结构动力响应的分段精细时程积分方法
王超;李红云;刘正兴
【期刊名称】《计算力学学报》
【年(卷),期】2003(020)002
【摘要】利用钟万勰等发展的精细时程积分方法,提出了解线性定常结构动力系统响应的分段精细时程积分方法.它能适用于各种激励作用下系统的动力响应.对于载荷项采用线性和两次多项式进行拟合,采用精细时程积分方法和叠代方法对动力响应进行计算.与传统的离散积分方法如纽马克方法和威尔逊方法等及状态方程直接积分方法进行数值比较,本方法具有很高的精度和计算效率.
【总页数】5页(P175-178,203)
【作者】王超;李红云;刘正兴
【作者单位】上海交通大学,建筑工程与力学学院,上海,200030;上海交通大学,建筑工程与力学学院,上海,200030;上海交通大学,建筑工程与力学学院,上海,200030【正文语种】中文
【中图分类】O332
【相关文献】
1.求解一维饱和土固结方程的时程精细积分方法 [J], 蓝林华;富明慧
2.基于精细时程积分方法的时域格林函数计算 [J], 童晓旺;任慧龙;李辉;单鹏昊
3.饱和两相介质动力固结问题时域求解的精细时程积分方法 [J], 李亮;高超;吴利华
4.大型结构动力响应的状态方程的Krylov精细时程积分法 [J], 陈臻林
5.结构动力响应的精细时程积分并行算法 [J], 李红云;沈为平
因版权原因，仅展示原文概要，查看原文内容请购买。
受迫振动之动力响应(含精细积分知识与程序)
式（9.20）对应的线性插值是一种粗糙的近似，也可采用其他插值形式，如： （1）多项式 （2）指数函数 （3）正弦或余弦函数 （4）上述函数的乘积 等。 对于齐次方程， F (t ) = 0 ，由式(9.15)得 f = 0 ，故计算出上式中的指数矩阵 T，就可 得到齐次微分方程的时程积分。即 V1 = TV0 , V2 = TV1 , L,V n = TVn −1 为求出其中的指数矩阵，把指数函数写成如下形式 exp( Hτ ) = [exp( H∆t / s)] s (9.22) (9.23)
(9.4)
其中 θ 是常数，其值介于 1/2 和 1 之间。 考虑到阻尼矩阵的关系式（9.2），在积分区间上的前后两个时刻存在动力学方程如 下 d 2 u0 du M + (αM + βK ) 0 + Ku0 = f 0 (9.5) 2 dt dt d 2u du M 21 + (αM + βK ) 1 + Ku1 = f1 (9.6) dt dt 联立上述方程，可得下面三个类推关系： 1 M + (β + θ∆t )K u1 = θ∆tf1 + (1 − θ )∆tf 0 + α + ∆ θ t (9.7) du0 1 1 α + θ∆t M + (β + θ∆t − ∆t )K u0 + θ M dt du1 1 1 − θ du 0 = (u1 − u0 ) − (9.8) dt θ∆t θ dt d 2u1 1 du1 du 0 1 − θ d 2u 0 = − (9.9) − θ∆t dt dt θ dt 2 dt 2 特别地，当 θ = 1 / 2 时，这种方法就是纽马克法中的“ β = 1 / 4 ”法，也是最常用的积 分方法。
结构动力学第五章
i = 0 ,1,2 ,L
而这种离散化正符合计算机存贮的特点。 • 与运动变量的离散化相对应，体系的运动微分方程也不一定要 求在全部时间上都满足，而仅要求在离散时间点上满足，这相 当于放松了对运动变量的约束。
采用等时间步长离散时，ti = iΔ t ,i = 1,2 ,3, L
&& & mui + cui + kui = Pi
• 根据是否需要联立求解耦联方程组，逐步积分法可分为两 大类：
– 隐式方法：逐步积分计算公式是耦联的方程组，需联立求 解，计算工作量大，通常增加的工作量与自由度的平方成正 比，例如Newmark-β 法、Wilson -θ 法。 – 显式方法：逐步积分计算公式是解耦的方程组，无需联立求 解，计算工作量小，增加的工作量与自由度成线性关系，如 中心差分方法（无阻尼时）。 • 下面先介绍分段解析算法，然后重点介绍两种常用的时域逐步积 分法—中心差分法和Newmark-β 法，同时也介绍Wilson -θ 法，最后介绍非线性问题分析方法。
5.2 分段解析法 （Piecewise Exact Method）
分段解析法对外荷载进行离散化处 理，假设在ti≤t≤t i+1时段内 P
实际荷载
P(τ ) = Pi + α iτ
Pi+1 Pi
插值荷载：P(τ)
α i = ( Pi +1 − Pi )/Δti
如果荷载P( t )采用计算机采样，即 离散数值采样，则以上定义可认为 是“精确”的。 • 分段解析法一般适用于单自由度体系动 力反应分析，对于多自由度体系，有时 可以采用等效方法在满足一定近似的条 件下将多自由度体系化为单自由度问题 进行分析，这时 也可以采用分段解析 法完成体系的动力反应分析。
混凝土结构静力和动力响应分析技术
混凝土结构静力和动力响应分析技术一、前言混凝土结构是建筑领域中最常见的结构形式之一，其性能直接关系到建筑的耐久性、稳定性和安全性。
因此，对混凝土结构的静力和动力响应分析技术的研究和应用具有重要的理论和实践意义。
本文将从混凝土结构的静力和动力响应分析基础、分析方法和应用实例三个方面进行详细阐述。
二、混凝土结构静力和动力响应分析基础1.混凝土材料的力学性质混凝土材料的力学性质是混凝土结构静力和动力响应分析的基础。
混凝土的材料性质包括弹性模量、泊松比、强度等。
其中，弹性模量是指材料在弹性阶段内应力和应变之比，泊松比是指材料在一定应力下垂直于应力方向的应变与应力方向应变之比，强度是指材料在破坏前所能承受的最大应力。
这些参数的确定需要进行试验和计算。
2.混凝土结构的力学模型混凝土结构的力学模型是指将结构抽象为一些理想化的杆件或板件，以便于进行力学分析。
混凝土结构的力学模型可以分为线性和非线性两种。
线性模型是指结构在弹性阶段内的力学行为可以用弹性理论描述，非线性模型是指结构在破坏前后力学行为都呈现出非线性特性。
3.混凝土结构静力分析方法混凝土结构静力分析方法根据结构的力学模型可以分为刚度法和力法。
刚度法是指通过计算结构的刚度矩阵和荷载向量，从而求解结构的内力和位移。
力法是指通过计算结构的受力平衡方程和变形方程，从而求解结构的内力和位移。
在实际工程中，通常采用有限元方法进行混凝土结构静力分析。
4.混凝土结构动力分析方法混凝土结构动力分析方法是指通过计算结构在地震、风等自然荷载下的响应，从而评估结构的抗震性能。
混凝土结构动力分析方法可以分为等效静力法和动力时程分析法。
等效静力法是指通过把地震荷载等效为静力荷载，从而进行静力分析。
动力时程分析法是指通过计算结构在时间上的响应，从而求得结构的内力和位移。
三、混凝土结构静力和动力响应分析方法1.混凝土结构静力响应分析方法混凝土结构静力响应分析方法可以采用有限元法进行计算。
海洋工程结构动力响应分析考核试卷
6.海洋工程结构动力响应分析中，多物理场耦合效应在浅水区域的结构中更为显著。（）
7.在进行海洋工程结构动力响应分析时，不需要考虑结构的疲劳累积损伤。（）
8.海洋工程结构动力响应分析中，地震载荷对浮式结构的影响比对固定结构的影响小。（）
A.时间历程分析法
B.长期统计分析法
C.疲劳累积损伤法
D.稳态响应法
17.以下哪些类型的海洋工程结构在动力响应分析中主要受到海流力的影响？（）
A.浮式结构
B.海底管线
C.潮汐能发电装置
D.锚链结构
18.在动力响应分析中，以下哪些方法可以用来模拟结构在复杂环境下的响应？（）
A.多物理场耦合分析法
B.多尺度分析法
B.海域地理条件
C.结构设计参数
D.海浪载荷
14.在海洋工程结构动力响应分析中，下列哪种方法主要用于研究结构的疲劳寿命问题？（）
A.疲劳累积损伤法
B.极限状态法
C.响应谱法
D.稳态响应法
15.下列哪种类型的海洋工程结构在动力响应分析中主要受到风载荷的影响？（）
A.深水油气平台
B.浮式风力发电装置
C.海底管线
D.实验研究
12.海洋工程结构动力响应分析中，以下哪些因素会影响结构的稳态响应？（）
A.结构的固有频率
B.载荷的频率
C.结构的阻尼比
D.海域的环境条件
13.以下哪些类型的海洋工程结构在动力响应分析中主要受到潮汐力的影响？（）
A.浮式结构
B.潮汐能发电装置
C.塔式结构
D.海底管线
14.在海洋工程结构动力响应分析中，以下哪些方法可以用来提高分析的精度？（）
于开平-结构动力学第十五讲
xt t
( xt t xt ) (1 ) xt (1 )txt t 2
K K a0 M a1C
将它们同时代入第三个方程，只剩下待求时刻的位移，整理得 Kxt t Qt t
Qt t Qt M (a6 xt a2 xt a3 xt ) C (a1 xt a4 xt a5 xt )
x(t ) lim
x(t t ) x(t ) xt t xt t 0 t t
1 x x x x xt t t t t t t 2 t t
x(t ) lim
x(t ) x(t t ) xt xt t t 0 t t
2.3 纽马克方法(Newmark method)
对待求的下一时刻的位移、速度和加速度在当前时刻������进行泰勒展开
1 1 x (t t ) x (t ) tx (t ) t 2 x (t ) t 3 x (t ) O(t 4 ) 2 6 1 x (t t ) x (t ) tx (t ) t 2 x (t ) O(t 3 ) 2 x (t t ) x (t ) x (tn ) O(t ) x(t t ) x(t ) tx (t ) O(t 2 ) t
(2) 确定初始值
x 0 , x 0 , x0
x0 = M 1 (Q (0) Cx0 - Kx0 )
1 1 2
(3) 选择时步长∆������ , 使它满足∆������ < ∆������������������ = ������������ /������(������������ 为系统的最小周期)
结 构 动 力 学
第五章 结构动力学中常用的数值算法
结构力学教学课件-10-2结构动力响应-文档资料
m y ( t ) c y ( t ) ky ( t ) F ( t ) P 其中 F ( t ) m ( t )
P
基底作水平运动时，对体系的作用效果相当 于在运动质量上沿加速度相反方向施加一个 等效水平力，振动方程的形式不变。
10.3 单自由度体系的自由振动
1 1 2 2 1
2
（1）低阻尼和无阻尼 （2）临界阻尼 （3）超阻尼
1 , 0
1 1
10.3 单自由度体系的自由振动
10.3.2低阻尼和无阻尼体系
2 1 1
当 （或 1 0 ） , 为正实数，有 1
1 i1 2 i2
得振动方程的通解：
i t i t t 1 1 y ( t ) e ( G e G e ) 1 2
利用欧拉公式
t 或 y ( t ) C e s i n ( t ) 1
振动方程的建立
10.2.1单自由度体系的力学模型
3) 阻尼体系
F c y ( t ) D
c阻尼常数，表示质点 m以 单位速度移动时所受阻 力 [力]长度时间
1
质量 — —弹簧 — —阻尼系统
振动方程的建立
• 以质量为m的隔离体作为研究对象，所受力为重力G和动 载荷 FP (t ) ,弹性回复力 F S ( t )和粘滞阻尼力 FD (t )，以及假想 作用于其上的惯性力 FI (t )。
上节课内容回顾
• • • • • • 动荷载的定义 动荷载的分类 动力自由度的概念 动力自由度的离散方法 频率、振型和阻尼的概念 单自由度体系的力学模型
上节课内容回顾
• • • • • 动荷载的定义 动荷载的分类 动力自由度的概念 动力自由度的离散方法 频率、振型和阻尼的概念 定义 在振动过程的任一时刻，确定体系全部质量 位置或变形状态所需的独立参数个数，称为体系 的自由度（degree of freedom，简记为DOF）。 振动自由度的求法附加支杆
结构动力学中的常用数值方法(精)
第五章 结构动力学中的常用数值方法5．1．结构动力响应的数值算法....0()(0)(0)M x c x kx F t x a x v ⎧++=⎪⎪=⎨⎪=⎪⎩当c 为比例阻尼、线性问题→模态叠加最常用。
但当C 无法解耦，有非线性存在，有冲击作用（激起高阶模态，此时模态叠加法中的高阶模态不可以忽略）。
此时就要借助数值积分方法，在结构动力学问题中，有一类方法称为直接积分方法最为常用。
所识直接是为模态叠加法相对照来说，模态叠加法在求解之前，需要对原方程进行解耦处理，而本节的方法不用作解耦的处理，直接求解。
（由以力学，工程中的力学问题为主要研究对象的学者发展出来的）中心差分法的解题步骤1. 初始值计算（1） 形成刚度矩阵K ，质量矩阵M 和阻尼矩阵C 。
（2） 定初始值0x ，.0x ，..0x 。
（3） 选择时间步长t ∆，使它满足cr t t ∆<∆，并计算 021()a t =∆，112a t =∆，202a a = （4） 计算...00101122tx x x x a a -∆=-+（5） 形成等效质量阵01M a M a C -=+ （6） 对M -阵进行三角分解T M LDL -= 2．对每一时间步长（1） 计算时刻t 的等效载荷201()()t t t t t Q Q Ka M x a Ma C x--∆=---- （2） 求解t t +∆时刻的位移 ()Tt t t LDL x Q -+∆=（3） 如需要计算时刻t 的速度和加速度值，则.1()t t t t t x a x x +∆-∆=- ..0(2)t t t t t t x a x x x +∆-∆=-+若系统的质量矩阵和阻尼矩阵为对角阵时，则计算可进一步简化。
纽马克法的解题步骤1.初始值计算（1）形成系统刚度矩阵K ，质量矩阵M 和阻尼矩阵C （2）定初始值0x ，.0x ，..0x 。
（3）选择时间步长t ∆，参数γ、σ。
结构动力方程的2种精细时程积分
结构动力方程的2种精细时程积分范宣华;陈璞;慕文品【摘要】Similarities and differences in solving dynamic equations between precise time-integration method with augmented matrix (PTI-AM) and extended precise time-integration method (EPTI) were analyzed. The explicit, discrete and recursive expressions for both methods were deduced, respectively, with the evolutionary random excitations in a general combined form of polynomial, exponential, and sinusoid/cosine functions. Both recursive expressions can be transformed into polynomial functions corresponding to the integral steps. With the same number of terms in the Taylor series, the recursive expression for EPTI contains additional high-order terms besides all the terms in PTI-AM. Therefore, EPTI has higher precision than PTI-AM does. If those additional high-order terms are neglected, the two methods have an identical discrete and recursive expression. In this respect, the two methods are essentially the same despite of different programming realization. An engineering example was presented, showing that the computing precision of both methods was as high as 10-significant-figures, and the computing time of EPTI was over 1 order of magnitude less than that of PTI-AM.%分析了增维精细时程积分和扩展精细时程积分2种方法在求解结构动力方程中的异同.在演变随机激励为多项式、指数函数以及正弦/余弦(虚指数)函数组合的一般形式下,分别推导出了2种方法所对应的显式离散递推表达式.结果表明:2种积分方法所对应的显式递推格式最终都转化为积分步长的多项式函数,并且在相同泰勒级数展开项的条件下,扩展精细积分除包含增维精细积分的所有递推项外,还包含一些高阶小项,理论上具有更高的精度；忽略高阶小项,2种方法尽管算法实现不同,离散递推格式完全一致；工程实例计算表明,二者计算精度都可以达到10位以上有效数字,扩展精细积分计算时间较增维精细积分少一个数量级.【期刊名称】《西南交通大学学报》【年(卷),期】2012(047)001【总页数】6页(P109-114)【关键词】细时程积分;演变随机激励;动力学方程;计算精度【作者】范宣华;陈璞;慕文品【作者单位】北京大学力学与空天技术系,北京100871;中国工程物理研究院总体工程研究所,四川绵阳621900;北京大学力学与空天技术系,北京100871;北京大学力学与空天技术系,北京100871【正文语种】中文【中图分类】O324;O241.4工程结构受到的外载作用大多是非平稳随机过程［1］.在进行结构非平稳随机激励响应分析中，通常采用均匀调制演变随机激励进行简化［2］.在实际应用中，简化后的模型因计算量庞大，其应用仍受到限制.文献［3-4］中提出的虚拟激励法将非平稳随机振动分析转化为确定性时间历程分析，在求得虚拟时程响应后，再转换成求解随机激励响应的功率谱，在保证求解精度的前提下大大降低了计算量.对于均匀调制演变随机激励f(t)可表示为式中:g(t)为慢变均匀调制函数;x(t)为平稳随机过程;t为时间.在采用虚拟激励法后，相应的虚拟激励(t)变为式中:Sxx(ω)为自功率谱密度;ω为圆频率.结构运动方程变为式中:y(t)为虚拟激励下的响应;M和K分别为结构质量、阻尼和刚度矩阵;E为虚拟激励沿各自由度方向分布的非时变向量.求解式(3)的方法有:中心差分法、Newmark法、Wilson-θ法等［5］，这些方法存在对积分步长敏感、精度不高等缺点.近年发展起来的精细时程积分法［6-15］，通过利用指数加法定理求解指数矩阵，使得计算精度得到提高.文献［10］中引入混合精细积分法(HPD-M)，对式(3)进行降阶求解，并对Duhamel积分解析求解，取得了很好的效果.但解析解中包含矩阵求逆，导致计算量增大，并且数值稳定性难以保证.基于混合精细积分发展起来的增维精细积分［8-9，15］和扩展精细积分［11-12］不但成功地避开了矩阵求逆带来的不稳定问题，还进一步提高了数值求解精度.本文对以上2种高精度数值积分方法进行了比较，针对演变随机激励的一般形式，推导出2种方法所对应的离散递推表达式.1 混合精细积分借助y(t)的一阶导数，式(3)可变为一阶非齐次微分方程式中:v0为v(t)在时刻t=0的初始向量;A、r、v(t)分别为利用Duhamel积分，式(4)的通解可表示为式中:τ为积分变量.对式(5)进行离散变换，得到v(t)第k步到第k+1步的离散递推关系式中:η为积分步长.利用指数函数的加法定理［10］，将积分步长内细分为2M份，在精细区段进行泰勒展开，然后利用增量存储得到exp(Aη)的精细解.增量存储指泰勒对指数函数精细区段展开后，对展开项采用M次平方运算，在每次运算时，将泰勒展开的第一项(单位阵)与后面几项小量分别存储，从而保持计算精度.对于式(6)中的Duhamel积分项，当虚拟激励为多项式、指数函数、正余弦函数或其组合形式时，可解析求出.由于解析解包含矩阵求逆运算，使数值求解变得不稳定.2 增维精细积分和扩展精细积分增维精细时程积分［8-9］又称为齐次扩容方法，通过增加矩阵维数把包含时间的外部激励纳入了精细积分的细化过程，使得运动方程变为一阶齐次方程，避免了矩阵求逆.该方法利用了常用的均匀调制函数的封闭性，即当均匀调制函数为多项式和指数函数(正余弦函数可以表示成复指数形式)乘积时，经过对时间有限次求导后的高阶导数，可表示成原函数及其低阶导数的实线性函数(一般而言，可以证明当多项式最高次数为k时，经过2(k+1)次求导后，乘积函数变成封闭).由此可以利用均匀调制函数及时间导数构造向量Γ(t)，使之满足式中:C0为行向量;D为状态矩阵;˙Γ(t)为状态向量Γ(t)对时间的一阶导数;Γ0为Γ(t)在时刻t=0的初始值.将式(4)和(7)合并，便可得到一阶齐次方程为其中，对于给定的精细区段积分步长η，可得到u(t)第k步到第k+1步的离散递推关系为式(9)可以采用指数函数加法定理直接求解，不用进行矩阵求逆.该法求解的关键在于确定状态向量Γ(t)、行向量C0和矩阵D.文献［8］给出常见调制函数g(t)对应以上向量和矩阵的表达式.扩展精细时程积分［11-12］仍然基于混合精细积分方法，不同的是扩展精细积分方法对式(5)中Duhamel积分项进行处理时，不是采用解析求解方法，而是直接采用精细积分方法处理，避免了解析解带来的矩阵求逆问题，保证了数值求解精度. 虚拟激励的一般形式为式中:n为多项式的次数;cn为相应多项式的系数;ρ为指数系数.将式(10)代入式(6)，利用多项式展开(t－τ)n，同时将指数矩阵进行泰勒展开成多项式并对多项式进行积分，最终可以得到精细区段离散递推关系式为式中:β= －ρ－iω;tk+1为k+1步对应的时刻t;I为单位矩阵.对于精细区段，泰勒展开4～5项就足够了.有了精细区段的离散递推关系，再利用加法定理和增量存储方法［11-12］，可求得给定时间步长内的递推关系.可以看出，以上2种积分方法均以混合精细积分方法为基础及避开矩阵求逆为目的，在求解过程中均用矩阵指数的加法定理和矩阵指数的泰勒级数展开定理.文献［6-11］中表明数值求解过程是稳定的，求得结果精度高.2种算法的思路和编程实现又完全不同.增维精细积分是以增大矩阵维数为代价来避开矩阵求逆，通过引入新的状态向量将激励纳入矩阵运算，从而去掉了动力方程的非齐次项，该方法需先根据激励形式推导状态矩阵和向量，并且方程维数有所扩大，加上每个区段都要进行一次矩阵指数的精细计算，显著增加了计算量.扩展精细积分方法同样避开Duamel积分的解析解，在精细区段，一般虚拟激励形式下，利用泰勒展开将Duhamel积分项内的指数矩阵变成多项式形式进行求解.由于多项式展开都是在精细区段进行，多项式展开的精度已经超过了计算机浮点运算精度，故多项式展开引入的误差不会对数值解产生影响，该方法同样成功地避免了矩阵求逆，在不扩大矩阵维数的同时又保证了求解的精度，与增维精细积分相比具有明显优势.3 2种积分的显式递推关系一般虚拟激励形式如式(10)所示，该表达式包括多项式、指数函数、正余弦函数及其组合的各种情形.为推导方便，本文中多项式只取两项，对于多项式包含更高此项的情形，推导过程完全一样.取虚拟激励方程为其中:c0和c1为多项式系数.在精细区段η∈(tk，tk+1)中，2种方法有如下递推关系.利用增维精细积分方法，确定Γ(t)、C0和矩阵 D.取Γ(t)第1项后面几项Γ2、Γ3等依次取前一项对时间的导数.取4项后的状态方程满足式(7)，各向量和矩阵表达式为将式(12)代入式(9)中，考虑到精细区段，可以对指数矩阵H进行泰勒展开，取前5项，有增维精细积分编程求解时，依次计算式(13)右侧的各矩阵，再根据加法定理得到所求区段响应值.本文为了与扩展精细积分进行比较，将H及其高次项按照式(8)中H的分块方式展开，并只取u(t)中的v(t)部分响应(Γ(t)部分自然满足)，得到的v(t)离散递推关系为式(14)已经对增维精细积分方法利用矩阵分块进行了降维处理，可以避免矩阵维数增大带来的数值求解时间过长问题，进而提高计算效率.这种分块降维方法在文献［15］中也有所体现.对于扩展精细积分，直接将虚拟激励代入到式(11)，经过化简最终可得式(15)中同样考虑精细区段η∈(tk，tk+1)，并利用了关系式tk+1=tk+η，对出现的指数矩阵采用泰勒展开，在泰勒级数展开全取5项的条件下，式(15)中的高阶小项最高次数达到η［11］.从式(14)和式(15)对应的2种方法显式离散递推表达式可以看出，扩展精细积分递推关系仅比增维精细积分递推关系多出了一个高阶小项，因此扩展精细积分理论上具有更高的数值求解精度.在忽略高阶小项的情况下，2种积分方法所对应的离散递推表达式是完全一致的.需要说明，以上离散递推关系仅仅考虑了精细区段，对于给定的积分步长，2种方法均是利用加法定理进行增量存储得到积分步长内的递推关系，因而在采用相同精细份数的情况下，积分步长内的递推关系也必然是相同的.4 算例分析以一幢70层高楼建筑在地面水平激励(如地震)作用下的顶层结构响应为例进行说明.初始设计时，为获取整体结构大致动力学特性，将高楼建筑各层等效成一个质量单元m，将楼层间的弯曲特性用一个弹簧-阻尼系统(k和 c)描述［16］.经过简化处理得到图1所示的70个自由度体系［12］，地面水平激励采用演变随机加速度激励¨xg(t)进行描述.图1 70个自由度体系Fig.1 A 70-freedom degree system取¨xg(t)=g(t)x(t)，其中，调制函数g(t)取为式中:S0、ωg、ζg分别为结构谱常量、固有频率和固有阻尼比，计算时分别取为:以上参数均采用量纲为1的形式.采用增维精细时程积分法(PTI-AM)和扩展精细时程积分法(EPTI)计算t∈［0，40］区间的顶层位移时变方差，取ω∈［0，200］，Δω =0.5，Δt=0.5.2种方法各个时刻的计算结果和CPU计算时间见表1.从表1可以看出，2种方法得到的计算结果完全相同(二者结果具有10位以上相同有效数字，限于长度在此只列出5位)，充分说明了二者具有很高的数值计算精度.另一方面，既然二者离散递推表达式在决定数值精度的低阶部分完全一致，其数值计算结果也必然相同.EPTI的计算时间比PTI-AM少一个数量级以上，具有更高的效率，这主要是因为PTI-AM在每个频点上都需计算1次增维的矩阵指数.表1 顶层位移时变方差结果和计算时间Tab.1 Results of time-varying variance on top-layer displacement and computing time算法 t/s 8.0 16.0 32.040.0CPU/s PTI-AM 1504.59631.07924.16801.9 700.2 EPTI1504.59631.07924.16801.9 13.05 结束语本文通过理论推导和数值算例，比较验证了增维精细积分、扩展精细积分在数值实现和计算效率方面的差异.结果表明，2种积分方法尽管在算法实现上不同，但在相同泰勒展开项的前提下，二者的离散解析表达式在略去高阶小项的情况下却是完全一致的，数值算例也表明二者的求解结果在10位有效数字内完全相同，扩展精细积分求解效率明显高于増维精细积分方法.参考文献:【相关文献】［1］李斌，刘学毅.客运专线铁道车辆随机振动特性［J］.西南交通大学学报，2010，45(2):191-196.LI Bin，LIU Xueyi.Random vibration property of highspeed railway vehicle in passeenger dedicated line［J］.Journal of Southwest Jiaotong University，2010，45(2):191-196.［2］TOC W S.Response statistics of discrete structures to non-stationary random processes［J］.Journal of Sound Vibration，1986，105(6):217-231.［3］林家浩.随机振动的虚拟激励法［M］.北京:科学出版社，2004:42-54.［4］王波，张海龙，徐丰.考虑拉索局部振动的斜拉桥地震时程响应分析［J］.西南交通大学学报，2008，43(4):441-446.WANG Bo，ZHANG Hailong，XU Feng.Seismic timehistory response analysis of cable-stayed bridge with local cable vibration［J］.Journal of Southwest Jiaotong University，2008，43(4):441-446.［5］赵丽滨，王寿梅.结构动力分析中时间积分方法进展［J］.力学与实践，2001，23(2):10-15.ZHAO Libin， WANG Shoumei. Progressoftime integration methods in structural dynamics analysis［J］.Mechanics and Engineering，2001，23(2):10-15.［6］钟万勰.结构动力方程的精细时程积分法［J］.大连理工大学学报，1994，34(2):131-136.ZHONG Wanxie.On precise time-integration method for structural dynamics［J］.Journal of Dalian University of Technology，1994，34(2):131-136.［7］ZHONG Wanxie.Review of a high-efficiency algorithm series for structural responses ［J］.Progress in Natural Science，1996，6(3):257-268.［8］张文首，林家浩，刘婷婷，等.结构非平稳随机响应的增维精细时程积分［J］.振动与冲击，2006，25(5):18-20.ZHANG Wenshou，LIN Jiahao，LIU Tingting，et al.Precise time-integration method with augmented matrix for analysis of non-stationary random responses［J］.Journal of Vibration and Shock，2006，25(5):18-20.［9］张文首，林家浩，于骁.受演变随机激励结构响应的增维精细时程积分法［J］.大连理工大学学报，2007，47(2):160-164.ZHANG Wenshou，LIN Jiahao，YU Xiao.Precise timeintegration method with augmented matrix for structural responses to evolutionary random excitations［J］.Journal of Dalian University of Technology，2007，47(2):160-164. ［10］LIN Jiahao， SHEN Weiping， WILLIAMSFW.Accurate high speed coinputation of nonstationary random seismic response［J］.Engineering Structures，1997，19(7):586-593.［11］谭述君，钟万勰.非齐次动力方程Duhamel项的精细积分［J］.力学学报，2007，39(3):374-381.TAN Shujun，ZHONG Wanxie.Precise integration method for duhamelterms arising from non-homogenous dynamic systems［J］.Chinese Journal of Theoretical and Applied Mechanics，2007，39(3):374-381.［12］幕文品.受演变随机激励结构响应的扩展精细积分方法［D］.北京:北京大学，2009:39-42. ［13］富明慧，刘祚秋，林敬华.一种广义精细积分法［J］.力学学报，2007，39(5):672-677.FU Minghui，LIU Zuoqiu，LIN Jinghua.A generalized precise time step integration method ［J］. Chinese Journal of Theoretical and Applied Mechanics，2007，39(5):672-677. ［14］富明慧，林敬华.精细积分法在非线性动力学问题中的应用［J］.中山大学学报，2008，47(3):1-5.FU Minghui， LIN Jinghua. Precise time step integration method in nonlinear dynamics［J］. Acta Scientiarum Naturalium Universitatis Sunyatseni，2008，47(3):1-5. ［15］张继峰，邓子辰.结构动力方程的增维分块精细积分法［J］.振动与冲击，2008，27(12):88-90.ZHANG Jifeng，DENG Zichen.Dimensional increment and partitioning precise integration method for structural dynamic equation［J］.Journal of Vibration and Shock，2008，27(12):88-90.［16］高淑英，沈火明.线性振动教程［M］.北京:中国铁道出版社，2003:1-5.。