一種新型高精度半顯式基于模型的積分算法

傅 博， 張付泰， 陳 瑾

（長安大學(xué) 建筑工程學(xué)院，陜西 西安 710061）

直接逐步積分方法被廣泛應(yīng)用于求解離散化的結(jié)構(gòu)動(dòng)力學(xué)運(yùn)動(dòng)方程。根據(jù)位移差分方程的表達(dá)式不同，可以將積分算法分為顯式和隱式積分算法。在求解非線性結(jié)構(gòu)動(dòng)力學(xué)問題時(shí)，顯式積分算法無需進(jìn)行迭代，節(jié)省了計(jì)算時(shí)間，因此具有較高的計(jì)算效率。根據(jù)積分時(shí)間步長選取是否有限制，積分算法又可以分為無條件穩(wěn)定和條件穩(wěn)定積分算法。當(dāng)求解具有多自由度的復(fù)雜結(jié)構(gòu)動(dòng)力學(xué)問題時(shí)，條件穩(wěn)定積分算法可能需要極小的積分步長以滿足穩(wěn)定，因此極大地降低了計(jì)算效率。綜上，如果一個(gè)算法同時(shí)兼具顯式和無條件穩(wěn)定，那么將具有非常高的計(jì)算效率。但是，通常顯式積分算法是條件穩(wěn)定，而無條件穩(wěn)定算法往往是隱式。

為了同時(shí)滿足顯式和無條件穩(wěn)定這兩個(gè)特性，近年來學(xué)者提出一類新型積分算法，因?yàn)槠浞e分參數(shù)與模型質(zhì)量、剛度和阻尼相關(guān)，故被稱為基于模型的積分算法［1］。基于模型積分算法計(jì)算效率高，主要用于混合試驗(yàn)［2］、實(shí)時(shí)混合模擬［3］、子結(jié)構(gòu)振動(dòng)臺(tái)試驗(yàn)［4］及倒塌模擬［5］中?；谀Ｐ偷姆e分算法根據(jù)其速度差分方程的不同又可以分為雙顯式和半顯式算法，雙顯式指的是其速度、位移差分方程均為顯式，半顯式指的是僅位移差分方程是顯式，但速度差分方程是隱式。Chang［6］最早提出一種基于模型的Chang 算法，該算法為半顯式算法，具有二階精度。隨后，Chen 和Ricles［7］提出了另一種基于模型的CR算法，該算法為雙顯式，同樣具有二階精度。Gui［8］等基于CR算法的表達(dá)式，提出一族雙顯式的Gui算法，Gui 算法與γ=1/2 的Newmark 族算法具有相同的數(shù)值特性，具有二階精度。Fu［9］等提出一族與Newmark 族算法［10］具有相同數(shù)值特性的雙顯式GCR算法，該族算法的精度不超過二階。

綜上，與傳統(tǒng)積分算法相比，基于模型的積分算法雖然在計(jì)算效率上有很大優(yōu)勢(shì)，但是在精度上和常規(guī)方法類似。為此本文提出一種高精度的基于模型的積分算法，該算法既保持了基于模型的積分算法的顯式和無條件穩(wěn)定的數(shù)值特性，同時(shí)具有更高的數(shù)值精度。

1 新算法的建立

1.1 已有積分算法簡介

對(duì)于線彈性的單自由度（SDOF）結(jié)構(gòu)體系，離散時(shí)間系統(tǒng)下結(jié)構(gòu)運(yùn)動(dòng)方程可表示為

式中：m，c，k分別為質(zhì)量，阻尼系數(shù)和剛度，c=2mωξ，其中ξ為阻尼比，ω為圓頻率，分別為i+1時(shí)刻的加速度，速度和位移；Fi+1為i+1時(shí)刻的外界激勵(lì)。

積分算法可用于求解式（1），以經(jīng)典的Newmark族積分算法為例，其速度、位移差分方程如下：

式中：Δt是積分時(shí)間步長，β，γ是積分參數(shù)，直接影響積分算法的數(shù)值特性。常用的Newmark 算法包括常平均加速度法（CAA，β=1/4，γ=1/2），線性加速度法（LA，β=1/6，γ=1/2），顯式算法（NE，β=0，γ=1/2）。由于β和γ與結(jié)構(gòu)模型無關(guān)，所以Newmark算法屬于模型不相關(guān)的積分算法。

與模型不相關(guān)的積分算法不同，基于模型的積分算法的積分參數(shù)與模型相關(guān)，以CR 算法為例，其速度、位移差分方程如下：

式中：α1，α2為與模型相關(guān)的積分參數(shù)，其值為

類似地，Chang 算法的速度、位移差分方程如下：

式中：α1，α2亦為與模型相關(guān)的積分參數(shù)，其值為

本文提出一種新型高精度半顯式基于模型的積分算法，該算法采用Chang算法相同的速度、位移差分方程，但是具有不同的積分參數(shù)。下文將對(duì)新算法的積分參數(shù)進(jìn)行推導(dǎo)。

1.2 離散控制理論方法

采用離散控制理論的方法來確定基于模型的積分算法的積分參數(shù)，該方法在文獻(xiàn)［7-9］已經(jīng)成功應(yīng)用。根據(jù)離散控制理論，式（7），（8）和式（1）的Z 變換形式如下：

聯(lián)立式（11）～（13），得

其中，分子分母系數(shù)如表1所示。

表1 離散傳遞函數(shù)的分子和分母系數(shù)Tab.1 Numerator and denominator coefficients of the discrete transfer function

一種稱為“零極點(diǎn)匹配”的離散化方法在控制理論中用來從連續(xù)域的極點(diǎn)映射到離散域的極點(diǎn)，這種精確的映射法則規(guī)定如下：

式中：Δt為采樣周期，可以令其與積分算法的時(shí)間步長相等。式（15）的映射法則不太適合實(shí)際應(yīng)用，因此通常采用近似的映射法則。

本文新算法采用2 階pade 近似其映射規(guī)則如下：

采用2階pade近似與式（15）的逼似程度遠(yuǎn)優(yōu)于用于CR算法和Chang算法的1階pade近似，對(duì)比情況如圖1 所示，可以看出，2 階pade 近似的效果要遠(yuǎn)優(yōu)于1階pade近似，尤其當(dāng)sΔt＜-0.5時(shí)。

圖1 兩種Pade近似與指數(shù)函數(shù)對(duì)比Fig.1 Comparisons of two Pade approximations with exponential function

1.3 新算法積分參數(shù)推導(dǎo)（SDOF體系）

采用2階pade近似的離散傳遞函數(shù)極點(diǎn)如下：

根據(jù)式（14）特征方程的零點(diǎn)，特征方程可表達(dá)為

將表1中離散傳遞函數(shù)的分母系數(shù)d2，d1，d0和式（18），同時(shí)代入式（19），可以聯(lián)立求解得到積分參數(shù)α1和α2的值為

1.4 新算法積分參數(shù)推導(dǎo)（MDOF體系）

對(duì)于具有n個(gè)自由度的多自由度（MDOF）線彈性系統(tǒng)，新算法的速度、位移微分方程和運(yùn)動(dòng)方程如下：

假定阻尼矩陣是經(jīng)典阻尼，由于結(jié)構(gòu)振型的正交性，可將式（22）， （23）和（24）寫成模態(tài)坐標(biāo)系下的表達(dá)式，即

因此，可以進(jìn)一步得到基于模型的積分參數(shù)矩陣α1和α2如下：

2 數(shù)值特性分析

2.1 精度分析

首先將式（1）， （7）和（8）寫成遞推矩陣的形式：

式中：A為放大矩陣，其值為

放大矩陣A的特征方程按照|A—λI|=0進(jìn)行計(jì)算如下：

式中：λ為放大矩陣A的特征值；A1、A2、A3的表達(dá)式為

離散時(shí)間系統(tǒng)的位移數(shù)值解與連續(xù)時(shí)間系統(tǒng)的位移精確解的差值稱為局部截?cái)嗾`差［11］，表示如下：

假設(shè)位移任何階導(dǎo)數(shù)均連續(xù)可微，則式（36）的有限階的泰勒級(jí)數(shù)展開如下所示：

式中：x(l)為位移的第l階微分，Tl的表達(dá)式如下所示：

這里假設(shè)L=5，可以得到新算法的局部截?cái)嗾`差表達(dá)式如下：

式中：

從系數(shù)d1，d2，d3可以看出，新算法具有4 階精度，而Chang、CR、CAA、LA、NE 等算法均為2 階精度。為了進(jìn)一步驗(yàn)證本文提出的新算法具有4 階精度，用無阻尼線彈性單自由度結(jié)構(gòu)的自由振動(dòng)來給出新算法的絕對(duì)誤差收斂速率，如圖2所示。取，計(jì)算時(shí)間tn=1s，計(jì)算次數(shù)N=10，100，1 000，10 000。

2.2 穩(wěn)定性分析

本文方法放大矩陣的特征值可由式（35）求得

本文方法的譜半徑如圖3所示。

圖3 新算法的譜半徑Fig.3 Spectra radius of new algorithm

由圖3可知，當(dāng)ξ＞0時(shí)，譜半徑ρ＜1，表明本文方法在有阻尼情況下無條件穩(wěn)定；當(dāng)ξ=0 時(shí)，譜半徑ρ=1，當(dāng)時(shí)，ε=0，意味著有兩個(gè)重根λ1=λ2=—1，表明在該特例下算法不穩(wěn)定，而對(duì)于其他Ω值，本文方法均穩(wěn)定。因?yàn)檎鎸?shí)結(jié)構(gòu)均存在阻尼，因此，本文方法仍可視作無條件穩(wěn)定的算法。

2.3 周期延長和振幅衰減

算法的周期延長（PE）和振幅衰減（AD）可以反映積分算法的數(shù)值準(zhǔn)確度。算法的周期延長PE定義如下：

圖4 振幅衰減的定義Fig.4 Definition of amplitude decay (AD)

圖5 多種積分算法的周期延長PE對(duì)比Fig.5 Comparison of period elongation (PE) of various integration algorithms

3 數(shù)值算例

3.1 無阻尼線彈性SDOF結(jié)構(gòu)的自由振動(dòng)

采用本文方法、Chang 和CR 算法對(duì)無阻尼線彈性SDOF 結(jié)構(gòu)的自由振動(dòng)進(jìn)行求解，采用時(shí)間步長為Δt=0.02s 和Δt=0.05s 進(jìn)行計(jì)算。該結(jié)構(gòu)質(zhì)量m=20kg，k=2 000N·s-1， 因此結(jié)構(gòu)的自振頻率ω=10rad·s-1。結(jié)構(gòu)的初始位移x0=0m，初始速度結(jié)構(gòu)模型如圖6所示。

圖6 無阻尼線彈性SDOF結(jié)構(gòu)的簡圖Fig.6 Schematic diagram of a linear elastic SDOF system without damping

無阻尼線彈性SDOF 結(jié)構(gòu)的自由振動(dòng)的位移解析解［12］如下：

圖7 為3 種積分算法采用兩種不同積分步長計(jì)算得到的位移時(shí)程曲線。由圖7可以看出，本文方法的計(jì)算結(jié)果與精確解非常接近，而CR、Chang算法的計(jì)算結(jié)果與精確解存在較為明顯的差異，尤其是當(dāng)Δt=0.05s時(shí)，由周期延長引起的相位誤差較大。

圖7 位移時(shí)程曲線Fig.7 Time history curves of displacement

進(jìn)一步，采用以下兩個(gè)指標(biāo)來衡量積分算法的誤差：

式中：xN和xR分別代表數(shù)值（Numerical）結(jié)果和參考（Reference）結(jié)果；N是樣本數(shù)目。NEE和NRMSE分別對(duì)幅值和相位誤差比較敏感。表2 列出了位移誤差指標(biāo)。

表2 位移誤差指標(biāo)Tab.2 Error indices of displacement%

由表2可以看出，相同時(shí)間步長的情況下，Chang算法的NEE小于CR 算法，二者的NRMSE接近，而本文方法的NEE和NRMSE均遠(yuǎn)小于CR 和Chang 算法；本文方法在較大時(shí)間步長（0.05s）時(shí)，誤差還不到CR 和Chang算法在較小時(shí)間步長（0.02s）時(shí)的8%。

3.2 兩層鋼框架有限元模型的地震響應(yīng)分析

選取某兩層鋼框架辦公樓，該建筑的層高為4m，跨度為9m，側(cè)向力由抗彎框架承擔(dān)，重力由倚靠柱承擔(dān)，抗彎框架與倚靠柱之間通過剛性樓板假定進(jìn)行關(guān)聯(lián)，其結(jié)構(gòu)如圖8所示。

圖8 兩層鋼框架結(jié)構(gòu)示意圖Fig.8 Schematic diagram of two-story steel frame structure

基于Matlab/Simulink［13］，對(duì)該結(jié)構(gòu)進(jìn)行了有限元建模：材料密度為7.8×103kg·m—3，采用理想彈塑性材料本構(gòu)關(guān)系，彈性模量2×1011Pa，屈服強(qiáng)度為3.45×108Pa。采用基于剛度的纖維梁單元模擬梁柱構(gòu)件［14］，單元的積分點(diǎn)數(shù)為5，截面纖維數(shù)為24；采用彈性梁柱單元模擬倚靠柱并考慮P-Δ 效應(yīng)。該結(jié)構(gòu)模型共包括31 個(gè)節(jié)點(diǎn)，30 個(gè)單元以及83 個(gè)自由度。采用集中質(zhì)量矩陣，阻尼矩陣采用瑞利阻尼，假定前兩階模態(tài)阻尼比為2%。

為了驗(yàn)證基于Matlab/Simulink所建立的有限元模型的可靠性，在OpenSEES中對(duì)該框架進(jìn)行建模：材料本構(gòu)采用uniaxialMaterial ElasticPP，梁、柱均采用dispBeamColumn單元，倚靠柱采用elasticBeamColumn單元，其余模型參數(shù)與Matlab/Simulink模型保持一致。Matlab/Simulink模型與OpenSEES模型前5階周期對(duì)比如表3所示。

表3 Matlab/Simulink 模型與OpenSEES 模型前5 階周期對(duì)比Tab.3 Comparison of first 5th order periods of the Matlab/Simulink model and the OpenSEES model

從表3 中可以看出，基于Matlab/Simulink 建立的有限元模型的前5 階周期與基于OpenSEES 建立的有限元模型的各階周期非常接近，誤差最大僅為1.118 5%，從而說明本文基于Matlab/Simulink 建立的有限元模型的可靠性。

將El Centro NS， 1940作為結(jié)構(gòu)的地震動(dòng)輸入，采用CR、Chang及本文方法進(jìn)行時(shí)程分析，同時(shí)采用Δt=0.001s的CAA算法作為參考解。圖9給出了不同時(shí)間步長（Δt=0.005s， 0.01s 和 0.02s）時(shí)各算法的首層位移時(shí)程曲線。

圖9 首層位移時(shí)程曲線Fig.9 Time history curves of displacement at first story

從圖9和表4可以看出，相同時(shí)間步長下，CR算法和Chang算法的計(jì)算結(jié)果很接近，誤差指標(biāo)也非常接近，而本文方法的結(jié)果與參考解更契合，NEE和NRMSE均小于CR算法和Chang算法，其中本文方法的NRMSE優(yōu)勢(shì)非常明顯，不及CR 算法和Chang 算法的5%。即便在較大時(shí)間步長（Δt=0.02s）時(shí)，本文方法的NEE也不到3%，NRMSE不到0.3%，這說明本文方法在較大時(shí)間步長時(shí)也可以達(dá)到較高精度。表5 給出了本算例中各算法在不同時(shí)間步長的計(jì)算時(shí)間，使用的計(jì)算機(jī)配置為 CPU Intel Core i5-6 300HQ @2.30GHz，內(nèi)存8G。

表4 位移誤差指標(biāo)Tab.4 Error indices of displacement

表5 各算法的計(jì)算時(shí)間對(duì)比Tab.5 Comparison of computation time of different algorithms

由表5 可以看出，當(dāng)時(shí)間步長相同時(shí)，本文方法與CR算法、Chang算法的計(jì)算時(shí)間接近，但是計(jì)算時(shí)間明顯小于隱式CAA 算法。結(jié)合表4、5 可知，Δt=0.02s時(shí)本文方法的NEE仍小于Δt=0.01s時(shí)CR算法和Chang 算法的NEE；Δt=0.02s 時(shí)本文方法的NRMSE仍遠(yuǎn)小于Δt=0.01s 時(shí)CR 算法和Chang 算法的NRMSE，這說明Δt=0.02s時(shí)本文方法的計(jì)算精度高于Δt=0.01s 時(shí)CR 算法和Chang 算法。Δt=0.01s 時(shí)CR 算法和Chang 算法分別耗時(shí)33.90s 和33.27s，而Δt=0.02s 時(shí)本文方法耗時(shí)僅為17.19s。這說明本文方法不僅節(jié)約了約50%的計(jì)算時(shí)間，還達(dá)到了更高的計(jì)算精度。從另一個(gè)角度來說，在保證相近的計(jì)算精度的前提下，本文方法的計(jì)算效率要高于CR算法和Chang算法。

3.3 多自由度非線性質(zhì)量-彈簧系統(tǒng)的動(dòng)力學(xué)問題

基于模型的積分算法非常重要的特性是無條件穩(wěn)定性和顯式表達(dá)式。無條件穩(wěn)定性意味著無需為了保證算法的穩(wěn)定性而可以選擇較大的時(shí)間步長，顯式格式則意味著求解非線性問題時(shí)無需迭代。因此，基于模型的積分算法計(jì)算效率高也對(duì)應(yīng)于上述兩個(gè)方面：一是較大的時(shí)間步長，從而計(jì)算步數(shù)較少；二是每一積分時(shí)間步的計(jì)算時(shí)間短，因?yàn)闊o需迭代。本節(jié)將選取圖10中具有N個(gè)自由度的非線性質(zhì)量-彈簧系統(tǒng)來說明本文方法在計(jì)算效率和計(jì)算精度方面的優(yōu)勢(shì)。

圖10 質(zhì)量-彈簧系統(tǒng)示意圖Fig.10 Schematic diagram of Mass-Spring system

在圖10的質(zhì)量-彈簧系統(tǒng)中，mi=150kg，ki=5.5×105N·m-1（i=1，2，3···N），x?g=10sin（10t） m·s—2（t=0···5s）， 阻尼矩陣采用瑞利阻尼，假定前兩階模態(tài)阻尼比為2%，采用式（46）的非線性剛度，則

式中：ki，t為t時(shí)刻的第i個(gè)彈簧的剛度；xi，t為t時(shí)刻的第i個(gè)質(zhì)量的位移；xi—1，t為t時(shí)刻的第i—1 個(gè)質(zhì)量的位移。考慮N=50，100，400，700，1 000 等5 種不同自由度數(shù)的工況，選取Δt=0.001s 的Newmark 顯式算法作為參考解，對(duì)比Δt=0.05s的CAA算法、振型疊加法、CR 算法、Chang 算法和本文方法的計(jì)算結(jié)果。表6 給出了各算法的計(jì)算時(shí)間和x50的誤差指標(biāo)。

表6 各算法的計(jì)算時(shí)間和位移誤差指標(biāo)對(duì)比Tab.6 Comparison of computation time and error indices of displacement of different algorithms

從表6 中可以看出，CAA 算法的計(jì)算時(shí)間遠(yuǎn)大于本文方法和其他方法，這是因?yàn)镃AA算法需要迭代，CAA 算法的NEE和NRMSE均大于本文方法，說明CAA 算法的計(jì)算精度和計(jì)算效率均遜于本文方法；振型疊加法的計(jì)算時(shí)間相對(duì)于本文方法和其他方法有明顯的優(yōu)勢(shì)，但是該方法在求解非線性問題時(shí)有較大的誤差，雖然在100自由度工況時(shí)，其NEE略小于本文方法，但是其NRMSE是本文方法的5倍，而在其余4 種自由度工況下，振型疊加法的NEE、NRMSE均大于本文方法；與CR 算法和Chang 算法相比，本文方法計(jì)算時(shí)間略長，這是因?yàn)樵谟?jì)算積分參數(shù)矩陣的過程中需要額外的計(jì)算時(shí)間，但是在計(jì)算精度上，本文方法的NEE和NRMSE基本上都要小于CR算法和Chang算法，只有在1 000 自由度工況下，CR 算法的NEE略小于本文方法，但是其NRMSE是本文方法的2.5 倍。綜上，本文方法在計(jì)算精度上有明顯的綜合優(yōu)勢(shì)。

4 結(jié)論

本文提出了一種新型高精度半顯式基于模型的積分算法。新算法與半顯式Chang 算法具有相同的速度、位移差分方程，采用二階Pade近似極點(diǎn)映射方法生成新算法的積分參數(shù)。對(duì)新算法的精度、穩(wěn)定性和周期延長和振幅衰減等數(shù)值特性進(jìn)行分析，發(fā)現(xiàn)本文方法具有四階精度，周期延長極小，并且沒有能量耗散。通過3個(gè)具有代表性的數(shù)值算例，進(jìn)一步論證了本文方法的優(yōu)越性。

作者貢獻(xiàn)聲明：

傅博：論文修改，數(shù)據(jù)分析；

張付泰：論文撰寫，編程計(jì)算；

陳瑾：論文修改，數(shù)據(jù)校核。