波動數(shù)值模擬的常加速度顯式算法

波動數(shù)值模擬的常加速度顯式算法孫明社;曲淑英;侯興民【摘要】給出了一種將隱式時域逐步積分算法轉(zhuǎn)換為顯式時域逐步積分算法的方法，避免了求解耦聯(lián)方程組，提高了計算效率.對于暫態(tài)波源作用下的彈性半空間，利用有限單元法劃分網(wǎng)格、建立結(jié)構(gòu)動力方程，并應(yīng)用Fortran語言對中心差分算法和平均常加速度顯式算法編程，求解脈沖荷載作用下的出平面運動.2種算法計算結(jié)果對比表明，常加速度顯式算法可以較好地應(yīng)用于工程波動數(shù)值模擬中.%Explicittimeintegralmethodofnumericalsimulationforengineeringwaveproblemisanimportanttopicinbothnationalandinternationalresearchwork.Anexplicittimeintegralmethodtransformedbycorrespondingim-plicittimeintegralmethod,avoidingsolvingcoupledequationsandimprovingcomputationalefficiency,ispresentedinthispaper.Structuraldynamicequationofelastichalfspaceundertransientwaveisintegratedusingthefiniteel-ementmethod,andtheout-planeresponseunderimpulseloadiscomputedbycentraldifferencemethodandcon-stantaccelerationexplicitmethodrespectively.Comparisonbetweenthecalculatedresultsofthetwomethodsshowsthattheconstantaccelerationexplicitmethodcanbeusedinengineeringwavenumericalsimulationeffectively.【期刊名稱】《煙臺大學學報（自然科學與工程版）》【年(卷)，期】2015(000)001【總頁數(shù)】5頁(P61-65)【關(guān)鍵詞】常加速度顯式算法;中心差分算法;出平面運動;工程波動【作者】孫明社;曲淑英;侯興民【作者單位】煙臺大學土木工程學院，山東煙臺264005;煙臺大學土木工程學院，山東煙臺264005;煙臺大學土木工程學院，山東煙臺264005【正文語種】中文【中圖分類】TU470在土木工程和地震工程領(lǐng)域有很多問題都可以歸結(jié)為波動問題，如強震地面運動、土-結(jié)構(gòu)相互作用、結(jié)構(gòu)的無損檢測與探傷等問題.因而，研究工程波的傳播與振動具有非常重要的意義.通過建立力學模型，工程波動問題一般就可轉(zhuǎn)化為偏微分方程的求解.由于實際工程中許多工程結(jié)構(gòu)建立的動力方程為非線性，很難獲得精確的解析結(jié)果，因而，數(shù)值方法便得到廣泛應(yīng)用.時域逐步積分法是結(jié)構(gòu)動力方程求解中的一種有效方法.根據(jù)是否需要求解耦聯(lián)方程組，時域逐步積分法又可分為顯式算法和隱式算法.隱式時域逐步積分算法的研究成果較多，如常平均加速度方法、Houblt方法、Newmark方法、Gurtin方法、Wilson-0方法、Park方法以及a方法等［1-5］.但是隨著結(jié)構(gòu)自由度數(shù)目的增大，求解耦聯(lián)方程組的工作量巨大，由于顯式時域逐步積分方法無需求解耦聯(lián)方程組，因而具有明顯的優(yōu)勢.廖振鵬、李小軍、周正華、侯興民等［6-20］對顯式時域逐步積分算法都進行了研究.本文在國內(nèi)夕卜學者研究成果的基礎(chǔ)上探討了平均常加速度隱式算法，即通過矩陣級數(shù)展開轉(zhuǎn)變?yōu)轱@式算法，最后將該顯式算法在工程波動數(shù)值模擬中進行應(yīng)用，并與文獻［21］的中心差分算法比較，說明常加速度顯式算法可以較好的應(yīng)用于工程波動數(shù)值模擬中.1隱式算法顯式化推導(dǎo)1.1基本方法對于時域逐步積分方法的廣義線性加速度算法，每向前計算一步都需求解耦聯(lián)線性方程組式中:M為質(zhì)量矩陣;C為阻尼矩陣;火為剛度矩陣;{u}i+1為ti+1時刻的位移;{}i+1為ti+1時刻的廣義荷載;a，b為常數(shù)，不同的加速度取不同的值.將式(1)中左側(cè){u}i+1的系數(shù)矩陣提出矩陣，可得對式(2)矩陣求逆矩陣，有將式(1)進行矩陣求逆、級數(shù)展開、級數(shù)截斷最終得到式(9)，完成了隱式算法轉(zhuǎn)化為顯式算法.式(9)可不必求解耦聯(lián)線性方程組，提高了計算速度，節(jié)省了計算時間.用上述方法，以平均常加速度方法為例，推導(dǎo)其隱式顯式化過程，并進行出平面運動的波動數(shù)值模擬.1.2平均常加速顯式化當Newmark邛逐步積分方法中的參數(shù)y、B取不同值時，該方法退化為平均常加速度法、線性加速度法、中心差分法等.當y=0.5、B=0.25時，為平均常加速法，即式(1)中的參數(shù)取值:a=4,b=2.公式可記作將式(10)與式(1)對比，應(yīng)用前述推導(dǎo)過程，可以得到平均常加速方法對應(yīng)的顯式化公式如下應(yīng)用式(13)、(14)可以分別得到離散時間點上的速度、加速度式(12)、(13)、(14)為平均常加速算法的顯式化公式，若給定初始時刻的初始位移{u}0、初始速度仇}0、初始加速度{"u}0，即循環(huán)使用公式(12)、(13)、(14)，得到所有離散時間點上的位移、速度、加速度.平均常加速方法顯式化的收斂條件，由式(7)可知取向量的“1”范數(shù)，求解上式，可得對于單自由度無阻尼體系E=0,帶入上式，得2出平面波動數(shù)值模擬2.1出平面問題如圖1,在直角坐標系oxyz下，彈性半空間內(nèi)，剪切模量p1=1，介質(zhì)密度P1=1，不考慮阻尼.波源為沿y軸方向深度hs處作用的暫態(tài)荷載：圖1出平面波源問題Fig.1Problemofout-planewave暫態(tài)波源荷載時間分布函數(shù)為三角形脈沖T(t):2.2平面波動數(shù)值模擬波動數(shù)值模擬的精度需要將有限單元的尺寸劃分的與波長相比足夠小，因而在建立動力方程時才可以將每一有限單元內(nèi)慣性體積力視作不隨空間變化的恒定量.為了實現(xiàn)動力方程的空間解耦，采用了集中質(zhì)量有限元模型.采用有限單元法將連續(xù)彈性半空間介質(zhì)離散化，用有限個單元組合體代替原本的連續(xù)介質(zhì).用平面正方形單元將彈性半空間介質(zhì)離散化，假定各個單元只在公共節(jié)點上相互鉸接.離散的正方形單元邊長即空間離散步長為:k=Ay=0.05;單元節(jié)點編號如圖2;單元剛度［15］如矩陣(23).圖2正方形單元節(jié)點編號Fig.2Nodenumberofsquareunit將彈性半空間沿y軸向下劃分n個網(wǎng)格;y軸左右兩側(cè)，沿x軸正負方向各劃分m個網(wǎng)格，見圖3.圖3彈性半空間網(wǎng)格離散Fig.3Discretemeshofelastichalfspace觀測點位移峰值如表1.有限元網(wǎng)格劃分m=40,n=40;時步At=0.0075s;暫態(tài)波源作用位置hs=0即表面荷載;考慮二階透射邊界.采用Fortran程序語言進行中心差分算法和平均常加速隱式轉(zhuǎn)換顯式算法編程計算，取(0,0)(0.5,0)(1.0,0)(0,0.5)(0,1.0)5個觀測點，得到暫態(tài)波源作用下的彈性半空間計算結(jié)果.相同時步下中心差分算法和常加速度顯式算法不同觀測點的位移曲線如圖4.表1觀測點位移峰值Tab.1Displacementpeakofobservationpoint觀測點中心差分算法常加速度顯式算法(0,0)0.3460.346(0.5,0)0.1940.194(1.0,0)0.1370.137(0,0.5)0.2230.223(0,1.0)0.1750.175圖4觀測點位移曲線Fig.4Displacementcurveofobservationpoint將式(12)中系數(shù)逆矩陣取不同展開階數(shù)下的各觀測點位移曲線如圖5.圖5不同階數(shù)下的觀測點位移曲線Fig.5Displacementcurveofobservationpointunderdifferentordernumber3結(jié)論通過對比中心差分和常加速顯式方法的計算結(jié)果分析，得到如下結(jié)論.相同的網(wǎng)格數(shù)目、時步、邊界條件下，2種算法的出平面波動的主要波形完全吻合，峰值相同.曲線后段雖然存在一定的差別，但是，工程中我們主要關(guān)注主要波形，因而可以滿足工程要求.考慮阻尼的情況下，中心差分算法的顯式優(yōu)勢不明顯，而常加速度顯式算法仍為高效率的顯式算法，具有一定優(yōu)勢.逆矩陣在單位矩陣處級數(shù)展開，階數(shù)不同，但觀測點的位移曲線基本重合，說明常加速顯式算法只需取前幾階就可滿足精度要求，且計算效率較高.參考文獻：NewmarkNM.Amethodofcomputationforstructuraldynamics[J].ProcASCE，1959，85(3):69-94.BatheKJ,WilsonEL.Stabilityandaccuracyanalysisofdirectintegrationmethod[J].EarthqEngandStructDynam,1973(1):283-291.HiberHM,HughesTJR,TaylerRL.Impronednumericaldissipationfortimeintegrationalgorithinsinstructuraldynamics[J].EathqEngandStructDynam,1977,5(3):283-292.HiberHM,HughesTJR.Collocationdissipationandovershootfortimeintegrationschemesinstructuraldynamics[J].EathqEngandStructDynam,1978,6(1):99-177.SubbaraJK,DokainishMA.Asurveyofdirecttimeintegrationmethodsincomputationalstructuraldynamics(II)[J].ComputersandStructures,1989,32(6):1387-1401.廖振鵬.近場波動問題的有限元解法[J].地震工程與工程振動，1984,4(2):1-14.李小軍，廖振鵬，杜修力.有阻尼體系動力問題的一種顯式差分解法[J].地震工程與工程振動,1992,12(4):74-79.鐘萬勰結(jié)構(gòu)動力方程的精細時程積分法[J].大連理工大學學報，1994,34(2):131-136.［9］ 李小軍.地震工程中動力方程求解的逐步積分方法［J］.工程力學，1996,13(2):110-118.［10］ 周正華，李山有，侯興民.阻尼振動方程的一種顯式直接積分方法［J］.世界地震工程，1999,15(1):41-44.［11］ 杜修力，王進廷.阻尼彈性結(jié)構(gòu)動力計算的顯式差分法［J］.工程力學，2000,17(5):37-43.［12］ 王進廷，杜修力有阻尼體系動力分析的一種顯式差分法［J］.工程力學，2002,19(3):109-112.［13］ 陳學良，金星，陶夏新.求解加速度反應(yīng)的顯式積分格式研究［J］.地震工程與工程振動，2006,26(5):60-67.［14］ 高小科，鄧子辰，黃永安.基于三次樣條插值的精細積分法［J］振動與沖擊，2007,26(9):75-82.［15］ 廖振鵬，劉恒，謝志南.波動數(shù)值模擬的一種顯式方法-—維波動［J］.力學學報，2009,41(3):350-359.［16］ 劉恒，廖振鵬.波動數(shù)值模擬的一種顯式方法-二維波動［J］.力學學報，2010,42(6):1104-1116.［17］ 謝志南，廖振鵬.波