(中國科學(xué)院 武漢巖土力學(xué)研究所 巖土力學(xué)與工程國家重點(diǎn)試驗(yàn)室，武漢 430071)

1 引言

在經(jīng)典力學(xué)中，動力學(xué)可以劃分為質(zhì)點(diǎn)動力學(xué)和剛體動力學(xué)兩個分支[1]。質(zhì)點(diǎn)動力學(xué)研究質(zhì)點(diǎn)的運(yùn)動規(guī)律，即質(zhì)點(diǎn)的時間演化規(guī)律[1]，也即質(zhì)點(diǎn)位置與時間的函數(shù)關(guān)系，是一個未知函數(shù)。質(zhì)點(diǎn)動力學(xué)的理想實(shí)現(xiàn)途徑是，(1) 將物理系統(tǒng)簡化描述為由牛頓第二定律和力函數(shù)兩部分組成的質(zhì)點(diǎn)動力系統(tǒng)； (2) 設(shè)置系統(tǒng)參數(shù)和定解條件，采用一定積分方法，得到未知函數(shù)的閉合形式解或者數(shù)值解。然而，除了一維諧振子[2,3]之外，單擺[2,3]和彈簧擺等大多數(shù)動力學(xué)系統(tǒng)都是非線性的運(yùn)動微分方程，很難求解[1]。

彈簧擺是一個二維或三維的諧振子，是一個非線性微分方程組，屬于一個不可積系統(tǒng)，得不到閉合形式解[4,5]。其非線性主要體現(xiàn)在多個自由度之間的運(yùn)動耦合和能量轉(zhuǎn)換[6-10]。

基于微振幅假定，一些學(xué)者推導(dǎo)出彈簧擺系統(tǒng)的形式解，如文獻(xiàn)[4，11]推導(dǎo)出級數(shù)形式解，文獻(xiàn)[10，12]推導(dǎo)出指數(shù)形式解。

由于彈簧擺的運(yùn)動解存在隨著系統(tǒng)參數(shù)變化的混沌現(xiàn)象，許多學(xué)者[6,12-15]在定常參數(shù)空間或時變參數(shù)空間中定位出穩(wěn)定平衡解、不穩(wěn)定平衡解、周期運(yùn)動解、準(zhǔn)周期運(yùn)動解、混沌運(yùn)動解和不穩(wěn)定運(yùn)動解之間的分岔點(diǎn)。

為了顯示特定解在相空間的吸引子幾何形態(tài)，很多學(xué)者都對特定模態(tài)或定常參數(shù)的彈簧擺系統(tǒng)進(jìn)行數(shù)值求解，李銀山等[14]用Gear法得到數(shù)值解，張揚(yáng)鍵[7]用 Runge -Kutta 法得到數(shù)值解。文獻(xiàn)[6，16]基于數(shù)值解分析，得到系統(tǒng)外部激勵參數(shù)最為敏感的結(jié)論。李云龍等[5]基于數(shù)值解，分析了系統(tǒng)的能量在擺動模式與振動模式之間不斷相互轉(zhuǎn)化的特征。鄭建龍等[8]用高階 Runge -Kutta 法求解降階的非線性動力學(xué)系統(tǒng)，繪制了彈簧擺系統(tǒng)內(nèi)共振耦合發(fā)生的區(qū)域圖。

彈簧擺對于系統(tǒng)參數(shù)和初始條件的極度敏感性，使得三維諧振子難以應(yīng)用于真實(shí)物理系統(tǒng)的動力學(xué)模擬。但是，非線性動力系統(tǒng)的吸引子是可以控制和改變的[17]。文獻(xiàn)[10，18]通過解析方法證明了懸掛支點(diǎn)的豎向激勵有可能使彈簧擺的擺動運(yùn)動消失。Eissa等[19]分析了位移阻尼器和速度阻尼器對激勵系統(tǒng)進(jìn)行穩(wěn)態(tài)控制的可行性。李欣業(yè)等[9]發(fā)現(xiàn)選擇適當(dāng)?shù)姆答伩刂茀?shù)(如加入時滯項)可以使不穩(wěn)定的零解變得穩(wěn)定。Butikov[20]設(shè)計了一個分段函數(shù)形式的摩擦阻尼器，在外部激勵條件下可以實(shí)現(xiàn)比速度阻尼器更好的扭擺系統(tǒng)穩(wěn)態(tài)控制。牛江川等[21]研究了分?jǐn)?shù)階 PID 控制器對單自由度線性振子自由振動的控制性能。

從微分方程數(shù)值解法角度，一維諧振子、一維Duffing振子[22]、四維蔡氏振子[23]和三維Lorenz振子等都可以得到相當(dāng)精確的數(shù)值解[24]。但數(shù)值解法中應(yīng)用的諸多積分技術(shù)難以應(yīng)用于大規(guī)模的數(shù)值計算。在彈性力學(xué)中，彈性體模型是質(zhì)點(diǎn)力學(xué)中彈簧諧振子的三維推廣[1]。近年來，為了實(shí)現(xiàn)材料特大變形、破碎甚至發(fā)生相變等問題的動力學(xué)模擬，無網(wǎng)格法[25,26]、物質(zhì)點(diǎn)法[27,28]和近場動力學(xué)法[29,30]等數(shù)值計算方法均繞開了材料的物理構(gòu)形的描述。

對于一維諧振子，馮康等[31]檢驗(yàn)出蛙跳法屬于哈密頓體系辛算法，而 Runge -Kutta 法和 Adams 法解得的相空間軌跡存在失真現(xiàn)象。然而，對于一維弱非諧振子，蛙跳法就不再是辛算法。由于一維諧振子向三維推廣時必然導(dǎo)致幾何非線性，所以對于哈密頓體系的三維諧振子，蛙跳法不是辛算法[31]。要使用蛙跳法對三維諧振子進(jìn)行數(shù)值積分，脫離哈密頓體系應(yīng)該是一種可行的途徑。事實(shí)上，很多物理過程遵循達(dá)朗貝爾原理、最小能量原理或最小勢能原理；這意味著它們不屬于哈密頓體系，而屬于耗散動力系統(tǒng)。加入瑞利耗散函數(shù)[2]的一維阻尼諧振子雖然屬于耗散動力系統(tǒng)，但其無法擺脫對系統(tǒng)參數(shù)的依賴性[2]。

即使是一維諧振子，在物理問題動力學(xué)模擬的應(yīng)用過程中，也會出現(xiàn)系統(tǒng)不穩(wěn)定現(xiàn)象。劉白伊酈等[32]在應(yīng)用線性諧振子模擬晶格熱運(yùn)動時，采用了聲子熱浴方法維持系統(tǒng)穩(wěn)定狀態(tài)；并指出 Nosé -Hoover 熱浴引入的耗散項在處理非熱運(yùn)動時，可能會使溫度(能量)過分的衰減。

通過適當(dāng)?shù)目刂坪透淖?，三維諧振子也有可能脫離對系統(tǒng)參數(shù)和初始條件的依賴。將三維諧振子表示為耗散動力系統(tǒng)而非哈密頓體系，由三維諧振子描述的物理模型，有可能得到穩(wěn)定的數(shù)值解。為此，本文引入耗散項，設(shè)計一種三維諧振子，通過數(shù)值試驗(yàn)來驗(yàn)證運(yùn)動解的穩(wěn)定性，以期為弦振動問題和材料力學(xué)非線性大變形問題的動力學(xué)模擬提供一條可行的實(shí)現(xiàn)途徑。

2 三維諧振子

2.1 運(yùn)動的量化描述

如圖1所示，質(zhì)點(diǎn)運(yùn)動在靜止參照系Oijk中進(jìn)行量化描述。質(zhì)點(diǎn)的位置矢量為r。質(zhì)點(diǎn)與多個相鄰質(zhì)點(diǎn)之間有兩兩相互作用，相鄰質(zhì)點(diǎn)的位置矢量為rj。質(zhì)點(diǎn)指向相鄰質(zhì)點(diǎn)的方向基矢量為ej(j=1,2,3,…)。

在任一時刻t，質(zhì)點(diǎn)與相鄰質(zhì)點(diǎn)之間的相互作用僅與距離lj相關(guān)

lj=‖rj-r‖

(1)

相互作用采用彈簧力來描述，彈簧力函數(shù)遵循虎克定律

(2)

式中Nj(j=1,2,3，…)為彈簧力矢量的模，eN j為彈簧力的方向基矢量，k為彈簧剛度，lj0為彈簧原長，lj為彈簧長度。彈簧力對應(yīng)的彈性勢函數(shù)為

(3)

2.2 運(yùn)動微分方程

三維諧振子的運(yùn)動微分方程由牛頓第二定律和力函數(shù)兩部分組成。將牛頓第二定律轉(zhuǎn)化為兩個一階常微分方程，矢量形式為

(4)

式中v為速度矢量，m為質(zhì)點(diǎn)的質(zhì)量，F(xiàn)為力函數(shù)，由定常力項F0、彈簧力項F1、定常耗散項F2和彈簧耗散項F3四項組成。

F0為定常力項，表示質(zhì)點(diǎn)所受的與時間和位置等其他物理量都無關(guān)的、大小和方向均為恒定的已知力，如重力和恒定外力。

F1為彈簧力項，由質(zhì)點(diǎn)所受的相鄰質(zhì)點(diǎn)的作用力進(jìn)行矢量求和，

(5)

圖1 量化描述方式

Fig.1 Quantitative description

定常力項與彈簧力項都屬于保守力，合稱為保守力項

F01=F01e01=F0+F1

(6)

式中F01為保守力的模，e01為保守力的方向基矢量。

F2為定常耗散項，函數(shù)形式為

(7)

式中μ2(0<μ2<1)為定常阻尼系數(shù)，ev為質(zhì)點(diǎn)速度的方向基矢量，e0為定常力的方向基矢量，e01為保守力的方向基矢量。

F3為彈簧耗散項，是與質(zhì)點(diǎn)關(guān)聯(lián)的所有彈簧的耗散力的矢量求和

(8)

式中F3j(j=1,2,3，…)為每個彈簧的耗散力，是一個分段函數(shù)

(9)

式中μ3(0<μ3<1)為彈簧阻尼系數(shù)。系統(tǒng)可能出現(xiàn)極端情況

e01×Nj=0

(10)

式中0為零向量。此時，增加一個隨機(jī)的微量擾動力，令

e01×Nj=0+10- neB

(11)

式中n為正整數(shù)，建議值為4～10，eB為一個垂直于Nj的隨機(jī)的方向基矢量。微擾力的作用是使質(zhì)點(diǎn)逃離不穩(wěn)定的奇點(diǎn)，副作用是使質(zhì)點(diǎn)無法靜止在穩(wěn)定平衡位置。另一方面，微擾力作用向系統(tǒng)中輸入了隨機(jī)的附加能量，雖然很小，但也會引起系統(tǒng)能量的非規(guī)律性波動或振蕩。從參數(shù)角度，n值越大，系統(tǒng)能量的非正常波動越??；n值越小，系統(tǒng)能量的非正常波動越大。微擾作用的引入，使三維諧振子得不到靜態(tài)平衡解。

2.3 定解問題的數(shù)值解法

應(yīng)用三維諧振子進(jìn)行動力學(xué)模擬，大致可以劃分為三個子過程，即建模、定解和求解。

一個物理系統(tǒng)在建模完成后，可以得到系統(tǒng)的物理參數(shù)，包括質(zhì)點(diǎn)的靜態(tài)坐標(biāo)和質(zhì)量、彈簧的剛度和原長，以及質(zhì)點(diǎn)與彈簧的拓?fù)潢P(guān)聯(lián)。

定解就是設(shè)定數(shù)值解的條件，得到一個定解問題。定解條件又可以分為阻尼系數(shù)、邊界條件(bcs)和初始條件(ics)。邊界條件主要包括位置邊界條件(bcsR)、速度邊界條件(bcsV)和定常力邊界條件(bcsF0或bcsA0)。初始條件主要包含位置初始條件(icsR)和速度初始條件(icsV)。阻尼系數(shù)包含定常阻尼系數(shù)μ2和彈簧阻尼系數(shù)μ3。設(shè)定定解條件之后，可以得到定解問題。

定解問題的數(shù)值求解采用蛙跳法。積分過程可以劃分為初始化、積分參數(shù)設(shè)定和迭代三個子過程。

初始化過程是指迭代變量(位置r、速度v和加速度a)的賦值。初始化過程的流程為

(1)v=icsV

(2)v=bcsV[v]

(3)r=icsR

(4)r=bcsR[r]

(5)a=Force[r,v]

其中，[]表示一個函數(shù)，bcsV[]表示速度邊界條件的賦值函數(shù)，bcsR[]表示位置邊界條件的賦值函數(shù)，F(xiàn)orce[]表示力函數(shù)。

積分參數(shù)包括步長和迭代總步數(shù)。步長的估算方法如下。在模型中，提取任意一個質(zhì)點(diǎn)及與其關(guān)聯(lián)的任一個彈簧形成一個組合。設(shè)質(zhì)點(diǎn)的質(zhì)量為m，初始速度的模為v，初始加速度的模為a(a≠0)，彈簧剛度為k，初始長度為l。簡振模式的周期估算公式為

(12)

單擺模式的周期估算公式為

(13)

這個組合的估算周期值為兩個周解的最小值

Tj=min.[Ts h,Ts p]

(14)

式中下標(biāo)j表示一種組合的編號。當(dāng)a=0時，忽略單擺模式的周期。估算模型中所有組合的周期，得到周期解集合{Tj}。積分步長估算公式為

Δt=STmin.{Tj}

(15)

式中ST為步長分率。當(dāng)v=0時，式(12,13)可以退化為簡振周期和單擺周期。

迭代過程包含8個步驟

(1)v=v+0.5aΔt

(2)r=bcsR[r]

(3) Δr=vΔt

(4) Δr=bcsRRelax[Δr]

(5)r=r+Δr

(6)a=Force[r,v]/m

(7)v=v+0.5aΔt

(8)v=bcsV[v]

其中，Δr為位移矢量，步驟(1),(3),(5)～(7)為蛙跳法的標(biāo)準(zhǔn)流程，步驟(2)和(8)為邊界條件的施加，bcsRRelax為一個位置松弛函數(shù)

Δrb=λbΔr

(16)

式中λb(0<λb<1)為松弛系數(shù)，Δr為一個時步內(nèi)邊界質(zhì)點(diǎn)在無外部約束條件下的位移矢量，Δrb為外部約束許可的邊界點(diǎn)位移矢量。這是一種邊界松弛技術(shù)，允許約束邊界產(chǎn)生適度的振蕩。作用是使與質(zhì)點(diǎn)關(guān)聯(lián)的彈簧力減少，從而適當(dāng)降低系統(tǒng)的能量，類似于應(yīng)力松弛技術(shù)。

3 數(shù)值試驗(yàn)

建立耗散彈簧擺模型、簡化弦和簡化梁三個模型，進(jìn)行動力學(xué)模擬試驗(yàn)。試驗(yàn)?zāi)康氖球?yàn)證三維諧振子數(shù)值解的運(yùn)動穩(wěn)定性。

運(yùn)動解就是質(zhì)點(diǎn)位置量在位形空間的運(yùn)動軌跡。但是，通常質(zhì)點(diǎn)瞬時速度量也是一個很重要的物理量。為了彌補(bǔ)位形空間無法表示質(zhì)點(diǎn)速度的缺陷，將質(zhì)點(diǎn)速度表示為連續(xù)梯度條，基于梯度條對球體進(jìn)行著色的方式表達(dá)質(zhì)點(diǎn)速度的時間演化過程。

質(zhì)點(diǎn)的運(yùn)動解不易直接用來量化評估諧振子系統(tǒng)的運(yùn)動穩(wěn)定性。為此，計算出質(zhì)點(diǎn)的動能和靜力勢能，以及彈簧的彈簧勢能，求和后得到系統(tǒng)的總能量。采用能量的時程變化，基于里雅普諾夫運(yùn)動穩(wěn)定性定義，來討論運(yùn)動解的穩(wěn)定性。

3.1 耗散彈簧擺模型

三維耗散彈簧擺模型與三維守恒系統(tǒng)彈簧擺模型是一致的，即一個圓柱螺旋拉伸彈簧連接兩個質(zhì)點(diǎn)。質(zhì)點(diǎn)的質(zhì)量均為1.2105 kg。圓柱螺旋拉伸彈簧來源于GB T 2088-2009《普通圓柱螺旋拉伸彈簧尺寸及參數(shù)》，材料直徑為1 mm，彈簧直徑為10 mm，彈簧節(jié)距為5 mm，有效長度(彈簧原長)為500 mm，材料彈模為200 GPa，泊松比為 0.3，根據(jù)Dym[33]的計算公式得到彈簧剛度為 0.095 N/mm。

一個質(zhì)點(diǎn)固定在參考系的原點(diǎn)，這相當(dāng)于模型的邊界條件。其位置bcsR與速度bcsV恒為零向量。兩個質(zhì)點(diǎn)的bcsA0設(shè)定為(0,0,-9.81) N/kg。

另一個質(zhì)點(diǎn)設(shè)定四種初始條件，依次編號為ICS1～I(xiàn)CS4。如圖2所示，ICS1和 ICS2 的彈簧水平放置，初始長度為250 mm和750 mm。ICS3和ICS4的彈簧豎直放置，初始長度為750 mm和1000 mm。四種初始條件的icsV均設(shè)定為零向量。四種初始條件的能量依次為2.969 J，2.969 J， 11.880 J 和0.000 J。模型最小勢能為-6.680 J，對應(yīng)的靜態(tài)平衡位置點(diǎn)為(0,0,-0.625) m。

定常阻尼和彈簧阻尼設(shè)定四種組合，再與初始條件進(jìn)行二次組合，形成7個定解問題列入表1。

積分參數(shù)中，松弛系數(shù)為1.0，步長分率為 1/50，總步數(shù)為1000步，計算的步長取兩位有效數(shù)字為13 ms(DSP1～DSP4)和14 ms(DSP5～DSP7)，模擬總時長為13 s和14 s。定解問題的數(shù)值解依次如 圖3 ～圖9所示。初始時刻和終止時刻的彈簧擺均進(jìn)行了顯示。質(zhì)點(diǎn)運(yùn)動軌跡投影至xoy，xoz和yoz平面上。

DSP1經(jīng)過近200步的無規(guī)則振蕩后，進(jìn)入一個近似于圓錐擺的準(zhǔn)周期運(yùn)動。運(yùn)動軌跡在xoy投影平面內(nèi)近似于圓周運(yùn)動，在xoz投影面內(nèi)近似于水平的往復(fù)運(yùn)動。如表1所示，DSP1只施加了定常耗散項，即只對重力勢能進(jìn)行耗散。在xoz平面的近水平運(yùn)動表明，定常耗散項消減了質(zhì)點(diǎn)的重力勢能的振蕩，使重力勢能趨于穩(wěn)定。

圖2 初始條件及勢場

Fig.2 Initial conditions and potential field

表1 定解問題

Tab.1 Definite solution problems

IDicsμ2μ3DSP1ICS10.50.0DSP2ICS10.00.5DSP3ICS10.00.0DSP4ICS10.50.5DSP5ICS20.50.5DSP6ICS30.50.5DSP7ICS40.50.5

圖3 DSP1運(yùn)動軌跡

Fig.3 Trajectory of DSP1

DSP2經(jīng)過近170步的無規(guī)則振蕩后，運(yùn)動軌跡開始呈現(xiàn)出一定的規(guī)律性。運(yùn)動軌跡在空間上呈缽狀，在xoy投影平面呈現(xiàn)為一個旋轉(zhuǎn)的橢圓，這種運(yùn)動特征可以歸為準(zhǔn)周期運(yùn)動。如表1所示，DSP2只施加了彈簧耗散項，即只對彈性勢能進(jìn)行耗散。彈簧耗散項消減了質(zhì)點(diǎn)的彈性勢能的振蕩，使彈性勢能趨于穩(wěn)定。

DSP3未施加耗散項，其運(yùn)動軌跡是二維的，完全平行yoz投影平面，在xoy和xoz投影平面內(nèi)呈一條直線，雖然位形空間是三維的。對比前100步，DSP1～DSP3在yoz投影平面的運(yùn)動軌跡是類似的，不同之處在xoz投影平面，定常阻尼使質(zhì)點(diǎn)向x軸負(fù)向偏移，彈簧阻尼使質(zhì)點(diǎn)向x軸正向偏移。兩個耗散項均可以將質(zhì)點(diǎn)的運(yùn)動軌跡由二維推向三維。從運(yùn)動速度角度來看，DSP1和DSP2在200步之后均小于1.6 m/s，而DSP3在靜止平衡位置附近的速度始終超過2.7 m/s。質(zhì)點(diǎn)在運(yùn)行至靜止平衡位置時，速度最大。兩個耗散項均通過將質(zhì)點(diǎn)推離靜止平衡位置的方式，降低了系統(tǒng)的動能。

圖4 DSP2運(yùn)動軌跡

Fig.4 Trajectory of DSP2

圖5 DSP3運(yùn)動軌跡

Fig.5 Trajectory of DSP3

圖6 DSP4運(yùn)動軌跡

Fig.6 Trajectory of DSP4

圖7 DSP5運(yùn)動軌跡

Fig.7 Trajectory of DSP5

圖8 DSP6運(yùn)動軌跡

Fig.8 Trajectory of DSP6

圖9 DSP7運(yùn)動軌跡

Fig.9 Trajectory of DSP7

DSP4～DSP6是三種不同初始條件下的數(shù)值解。計算結(jié)果表明，在400步之后，三個運(yùn)動解均為類似于圓錐擺運(yùn)動特征的準(zhǔn)周期解。表明三維諧振子可以得到定態(tài)解。其中，如果缺少式(11)描述的微擾作用，DSP6在施加兩個耗散項的情況下，依然會處于一維的振動狀態(tài)，這樣會導(dǎo)致系統(tǒng)出現(xiàn)多解現(xiàn)象。微擾作用使質(zhì)點(diǎn)在不穩(wěn)定奇點(diǎn)附近的一維運(yùn)動產(chǎn)生偏心，從而進(jìn)入三維運(yùn)動狀態(tài)。三維諧振子的微擾描述，是與真實(shí)的物理現(xiàn)象相吻合的，一個真實(shí)的彈簧擺運(yùn)動不可能穩(wěn)定在一個不穩(wěn)定點(diǎn)的奇點(diǎn)附近。

DSP7的數(shù)值解在xoy投影平面內(nèi)的振蕩幅度很小，最大值不超過30 mm，也是一個準(zhǔn)周期解。與DSP6一樣，如果不施加微擾，DSP7的運(yùn)動解僅是一維空間的振蕩運(yùn)動。微擾作用將質(zhì)點(diǎn)推出靜止平衡位置，進(jìn)入類似于旋進(jìn)橢圓的準(zhǔn)周期運(yùn)動。豎直向和水平向微幅振蕩運(yùn)動的切換是彈簧擺內(nèi)共振現(xiàn)象的典型特征。由此，三維諧振子很有可能可以用來描述彈簧擺的內(nèi)共振現(xiàn)象。

三維諧振子的數(shù)值解是穩(wěn)定的。DSP1～DSP2(圖10)和DSP4～DSP7(圖11)6個定解問題，其能量在運(yùn)動過程中均持續(xù)減少,在接近最小能量值時趨于穩(wěn)定。其中，DSP7在0.5 s時能量變化趨于平穩(wěn)，穩(wěn)定速度最快；DSP1～DSP2和DSP4～DSP5的初始能量都為2.969 J，都在2.0 s時能量變化趨于平穩(wěn)；DSP6的初始能量最大，在4.0 s時能量變化才趨于平穩(wěn)。DSP3的定常阻尼和彈簧阻尼均為零值，其能量在運(yùn)動過程中依然在減少。DSP3屬于能量守恒系統(tǒng)，其能量減少主要是因?yàn)閿?shù)值解法中的邊界松弛技術(shù)。在數(shù)值解法中，對不動質(zhì)點(diǎn)(即懸掛支點(diǎn))的位移松弛，消減了彈簧的彈性勢能，起到了一定的體系消能的作用。

三維諧振子滿足能量最小原理。7個定解問題的最小勢能是-6.680 J(圖2)。由于初始速度均為零，系統(tǒng)的最小能量也是-6.680 J。系統(tǒng)在運(yùn)動過程中能量持續(xù)減少并趨于平穩(wěn)，是能量最小原理的描述內(nèi)容。因此，在三維諧振子的力函數(shù)中，定常耗散項和彈簧耗散項的設(shè)置具備一定的理論基礎(chǔ)。

3.2 簡化弦模型

從一維波動方程提出開始，弦振動問題一直是物理學(xué)的經(jīng)典問題之一。經(jīng)典弦理論采用偏微分方程對弦振動問題進(jìn)行描述，應(yīng)用邊值問題解法進(jìn)行積分求解。本節(jié)將驗(yàn)證三維諧振子應(yīng)用于弦振動問題動力學(xué)模擬的可行性。同時，為了驗(yàn)證數(shù)值解的穩(wěn)定性，定解問題的設(shè)定將以蠻力為特征，不考慮材料彈性屈服，也不模擬真實(shí)的加載過程。

建模的弦材料采用鋼絲。鋼絲的直徑為 1 mm，長度為500 mm，材料密度為7.9 g/cm3，材料彈模為200 GPa。將弦進(jìn)行離散化，生成一個七質(zhì)點(diǎn)六彈簧的簡化弦模型。彈簧長度約為0.083 m，彈簧剛度為1.88 kN/mm。將每段彈簧的分布質(zhì)量二分至與其關(guān)聯(lián)的質(zhì)點(diǎn)，得到質(zhì)點(diǎn)的質(zhì)量約為0.26 g(端部兩質(zhì)點(diǎn))和0.52 g(中間五質(zhì)點(diǎn))。

將端部兩質(zhì)點(diǎn)的位置固定，作為邊界條件(圖12)，其速度設(shè)定為零向量。不考慮預(yù)張拉，兩個邊界質(zhì)點(diǎn)的距離為0.5 m，位置坐標(biāo)分別為(0,0,0)和(0,0.5,0)，放置于Y軸正半軸上。此時，模型的彈性勢能為0 J。

定常阻尼和彈簧阻尼均設(shè)定為0.5。基于邊界條件和初始條件的設(shè)定，形成三個定解問題DSP8～DSP10。

圖10 DSP1～DSP3能量耗散

Fig.10 Energy dissipation of DSP1～DSP3

圖11 DSP4～DSP7能量耗散

Fig.11 Energy dissipation of DSP4～DSP7

DSP8模擬拉弦動作，僅將中間質(zhì)點(diǎn)由(0,2.5,0)拉至(0,2.7,0.1)并作為固定位置邊界條件，所有質(zhì)點(diǎn)的初始位置如圖13所示。icsV，bcsV和bcsA0 均為零向量。初始狀態(tài)下，系統(tǒng)的能量達(dá)到4602.5 J；質(zhì)點(diǎn)在運(yùn)動過程中最高速度接近1800 m/s，顯然其屬于蠻力型定解條件。

DSP9模擬放弦動作，以DSP8的終止?fàn)顟B(tài)為初始條件，并在邊界條件中剔除中間質(zhì)點(diǎn)的位置約束。DSP8和DSP9組合形成一個完整的拉放弦的動力學(xué)模擬。

DSP10模擬重弦。設(shè)定所有質(zhì)點(diǎn)的重力加速度bcsA0為(0,0,-9.81)×106m/s2，icsV和bcsV設(shè)定為零向量，其他的初始條件和邊界條件如圖14所示。中間的五個質(zhì)點(diǎn)的重力均超過 5 kN，總重力超過25 kN，其也屬于蠻力型定解條件。DSP10不考慮材料的彈性屈服。為了方便作圖，將重力勢能向z軸負(fù)方向移動0.06 m，作為重力零勢面。

積分參數(shù)中，松弛系數(shù)為1.0，步長分率為 1/10，總步數(shù)為2000步，計算的步長取兩位有效數(shù)字，均為7.3 μs，模擬總時長為14.6 ms。三個定解問題的數(shù)值解依次如圖15～圖17所示。初始時刻和終止時刻的弦線均進(jìn)行了顯示。質(zhì)點(diǎn)運(yùn)動軌跡投影至xoy和yoz平面上。在多質(zhì)點(diǎn)模型中，質(zhì)點(diǎn)之間的相對位置關(guān)聯(lián)可以描述弦線的振型。為了在質(zhì)點(diǎn)數(shù)量很少的前提下表示出振型，對初始時刻和終止時刻的弦線進(jìn)行了三維曲線擬合。

圖12 有界弦模型

Fig.12 Bounded string model

圖13 DSP8初始條件和邊界條件

Fig.13 Ics and bcs of DSP8

圖14 DSP10初始條件和邊界條件

Fig.14 Ics and bcs of DSP10

DSP8的運(yùn)動解(圖15)在位形空間中呈現(xiàn)為有界的近似于圓盤狀的二維空間區(qū)域。從右側(cè)兩質(zhì)點(diǎn)的空心區(qū)域可以看出圓盤面的中心位置質(zhì)點(diǎn)的靜止平衡位置。從xoy投影面可以看出，圓盤面近似垂直于靜止?fàn)顟B(tài)下的弦線。隨著能量的耗散，四個質(zhì)點(diǎn)的速度逐漸降低，并逐漸向靜止平衡位置靠攏。

DSP9的運(yùn)動解(圖16)在xoy和yoz投影平面內(nèi)均呈現(xiàn)為有界的窄條狀分布，且總體上垂直于y軸，即垂直于靜止平衡狀態(tài)的弦線。質(zhì)點(diǎn)存在沿y軸的振動分量，但振幅較小。弦線的振動主要在xoz平面上。這種特征與真實(shí)的有界弦振動的觀察結(jié)果比較吻合。

圖15 DSP8運(yùn)動軌跡

Fig.15 Trajectory of DSP8

圖16 DSP9運(yùn)動軌跡

Fig.16 Trajectory of DSP9

圖17 DSP10運(yùn)動軌跡

Fig.17 Trajectory of DSP10

DSP10的運(yùn)動解(圖17)在前70個時步高速向靜止平衡位置運(yùn)動，運(yùn)動主要發(fā)生在yoz平面上。在70步之后，五個質(zhì)點(diǎn)在yoz平面上振幅已經(jīng)很小，弦線的形狀在yoz平面上近似于懸鏈線，運(yùn)動轉(zhuǎn)換為近平行于x軸方向的水平振蕩，振幅最大為20 mm，且振幅持續(xù)收窄，在800步之后，振幅最大為10 mm。

由三維諧振子描述的簡化弦模型的數(shù)值解是穩(wěn)定的。DSP8～DSP10(圖18)的能量在運(yùn)動過程中均持續(xù)減少,并趨于穩(wěn)定。其中，DSP9在0.5 ms時能量變化就趨于平穩(wěn)，其初始能量最低，穩(wěn)定速度也最快；DSP10在4.0 ms時能量變化趨于平穩(wěn)；DSP10的能量在持續(xù)耗散，且耗散速率也在持續(xù)減少，顯現(xiàn)出趨于平穩(wěn)的趨勢。

3.3 簡化梁模型

材料大變形、幾何非線性和破裂一直是材料力學(xué)的研究熱點(diǎn)，而梁模型是最為常用的數(shù)值驗(yàn)證模型。現(xiàn)有理論和方法大多采用偏微分方程的邊界問題解法，屬于擬靜力模擬。本節(jié)將驗(yàn)證三維諧振子應(yīng)用于材料力學(xué)進(jìn)行動力學(xué)模擬的數(shù)值穩(wěn)定性。為了驗(yàn)證數(shù)值解的穩(wěn)定性，定解問題的設(shè)定也將以蠻力為特征，不考慮材料彈性屈服，也不模擬真實(shí)的加載過程。

梁的建模采用正方形截面鋼材。截面尺寸為10 cm×10 cm，長度50 cm，材料密度為7.9 g/cm3，材料彈模為200 GPa。將梁離散成四面體單元之后，為了得到三維諧振子的物理參數(shù)，將單元體的體積分配至與單元體關(guān)聯(lián)的棱，再分配至與單元體關(guān)聯(lián)的頂點(diǎn)。分配方案如下。首先，以四面體的四個關(guān)聯(lián)面的面積作為權(quán)重，將單元體的體積分配至四個關(guān)聯(lián)的三角形面；其次，以三角形的三個關(guān)聯(lián)棱的棱長作為權(quán)重，將三角形分配到的體積二次分配至三個關(guān)聯(lián)的棱；然后，采用二分法，將棱分配到的體積平均分配至兩個關(guān)聯(lián)的頂點(diǎn)；最后，每個頂點(diǎn)和每條棱分配得到體積進(jìn)行求和。這樣，每一條棱和每一個頂點(diǎn)都有了相應(yīng)的體積量。

圖18 DSP8～DSP10能量耗散

Fig.18 Energy dissipation of DSP8～DSP10

頂點(diǎn)的體積和材料密度相乘，得到質(zhì)點(diǎn)的質(zhì)量。彈簧原長直接等于棱的長度。棱的體積除以棱長，得到棱的截面面積，通過計算得到彈簧剛度。建模完成后，形成一個以三維諧振子描述的24個質(zhì)點(diǎn)和75條彈簧的簡化梁模型。模型的彈性勢能為0 J。

彈簧阻尼設(shè)定為0.3。bcsA0設(shè)定為零向量，定常阻尼系數(shù)可以忽略，不進(jìn)行設(shè)定?；谶吔鐥l件和初始條件的設(shè)定，形成三個定解問題DSP11～DSP13。三個定解問題中，相同的定解條件有，(1) 模型一端的四個質(zhì)點(diǎn)進(jìn)行位置固定，作為位置邊界條件；(2) bcsV，bcsA0和icsV均為零向量。

DSP11模擬梁的拉伸行為。將梁的另一個端部四個質(zhì)點(diǎn)的位置沿y軸拉長0.1 m并固定，作為位置邊界條件。按彈性力學(xué)的應(yīng)變概念，一些彈簧的線應(yīng)變達(dá)到了1.0，遠(yuǎn)超現(xiàn)有大變形模擬方法的描述范圍，可以稱之為蠻力加載。質(zhì)點(diǎn)的初始位置如圖19所示。

DSP12模擬梁的剪切行為。將梁的另一個端部的四個質(zhì)點(diǎn)的位置直接沿z軸平移0.1 m并固定，作為位置邊界條件，質(zhì)點(diǎn)的初始位置如圖20所示。

DSP13模擬梁的扭轉(zhuǎn)行為。將梁的另一個端部的四個質(zhì)點(diǎn)的位置以y軸正向?yàn)樾D(zhuǎn)軸按右手螺旋方向旋轉(zhuǎn)45°并固定，作為位置邊界條件。與拉伸和剪切相比，扭轉(zhuǎn)加載本身是幾何非線性的，是對稱張量難以描述的加載行為。運(yùn)動質(zhì)點(diǎn)的初始位置如圖21所示。

圖19 DSP11初始條件和邊界條件

Fig.19 Ics and bcs of DSP11

圖20 DSP12初始條件和邊界條件

Fig.20 Ics and bcs of DSP12

積分參數(shù)中，松弛系數(shù)為0.5，步長分率為1/10，總步數(shù)為300步，計算的步長取兩位有效數(shù)字后依次為15 μs，16 μs和16 μs。三個定解問題的數(shù)值解依次如圖22～圖24所示。為了顯示變形特征，對初始時刻和終止時刻的簡分梁的角棱線進(jìn)行了三維曲線擬合，并投影在yox和yoz平面上。

在拉伸作用下，DSP11(圖22)經(jīng)歷了10個時步，16個質(zhì)點(diǎn)均已運(yùn)動至靜止平衡位置附近。而后，質(zhì)點(diǎn)的運(yùn)動主要表現(xiàn)為xoz平面的振蕩。在100時步之后，振蕩的幅值均已小于3.0 mm，簡化梁進(jìn)入穩(wěn)定的位形。在拉伸作用下，簡化梁的變形呈現(xiàn)了非線性特征。從初始時刻和終止時刻的角棱線擬合曲線可以看出，簡化梁產(chǎn)生了扭曲變形。

在剪切作用下，DSP12(圖23)經(jīng)歷了100個時步后質(zhì)點(diǎn)運(yùn)動至靜止平衡位置附近區(qū)域，期間最高運(yùn)動速度接近1400 m/s。在100時步之后，振蕩的幅值均已小于2.0 mm，簡化梁進(jìn)入穩(wěn)定的位形。在剪切作用下，簡化梁的變形是非線性的。四條角棱在兩個投影平面都是曲線，且四組端點(diǎn)重合的角棱線的投影曲線均不重合。

圖21 DSP13初始條件和邊界條件

Fig.21 Ics and bcs of DSP13

圖22 DSP11運(yùn)動軌跡

Fig.22 Trajectory of DSP11

圖23 DSP12運(yùn)動軌跡

Fig.23 Trajectory of DSP12

在扭轉(zhuǎn)作用下，DSP13(圖24)經(jīng)歷了100個時步進(jìn)入靜止平衡位置附近區(qū)域，期間最高運(yùn)動速度接近520 m/s。在100時步之后，振蕩的幅值均已小于0.6 mm，簡化梁進(jìn)入穩(wěn)定的位形。扭轉(zhuǎn)作用本身是一種非線性的作用，終止時刻四條角棱在兩個投影平面雖然都是三維空間曲線，但線形均較為平直。

由三維諧振子描述的簡化梁模型的數(shù)值解是穩(wěn)定的。DSP11～DSP13(圖25)的能量在運(yùn)動過程中均持續(xù)減少,且都在1.0 ms時能量變化趨于平穩(wěn)。其中，DSP11的初始輸入能量最大(25.3 MJ)，在穩(wěn)定運(yùn)動狀態(tài)時系統(tǒng)保持的能量也最大(3.0 MJ)；DSP12的初始輸入能量為11.2 MJ，在穩(wěn)定運(yùn)動狀態(tài)時系統(tǒng)保持的能量為0.3 MJ；DSP13的初始輸入能量最小(1.1 MJ)，在穩(wěn)定運(yùn)動狀態(tài)時系統(tǒng)保持的能量也最小(0.16 MJ)。

圖24 DSP13運(yùn)動軌跡

Fig.24 Trajectory of DSP13

圖25 DSP11～DSP13能量耗散

Fig.25 Energy dissipation of DSP11～DSP13

DSP11～DSP13進(jìn)入穩(wěn)態(tài)的時間均在1.0 ms之內(nèi)，表現(xiàn)出三維諧振子對外部作用的高速響應(yīng)能力，這與物理試驗(yàn)中觀察到的梁材料對外部荷載的響應(yīng)速率是吻合的。

三個定解問題中，系統(tǒng)保持能量的方式都體現(xiàn)為質(zhì)點(diǎn)在平衡位置附近的微幅振動，這種特征符合晶格振動現(xiàn)象的物理描述。將材料的網(wǎng)格剖分按晶格構(gòu)形進(jìn)行離散化，三維諧振子有可能實(shí)現(xiàn)晶格振動現(xiàn)象的三維動力學(xué)模擬。

4 結(jié) 論

本文提出一種三維諧振子，并驗(yàn)證了其應(yīng)用于動力學(xué)模擬時具有良好的運(yùn)動穩(wěn)定性。

陳述了三維諧振子的運(yùn)動描述方式、運(yùn)動微分方程和數(shù)值積分流程。其中，在運(yùn)動微分方程的力函數(shù)中，引入了定常耗散項和彈簧耗散項；在蛙跳積分法中，引入了邊界松弛技術(shù)。

應(yīng)用三維諧振子，建構(gòu)了耗散型彈簧擺、簡化弦和簡化梁等三個模型，設(shè)定了13個定解問題進(jìn)行動力學(xué)模擬，得到運(yùn)動解。結(jié)果表明運(yùn)動解是穩(wěn)定的。

雖然數(shù)值解穩(wěn)定性問題是動力學(xué)研究領(lǐng)域的核心問題，但是，本文只驗(yàn)證了三維諧振子的可實(shí)現(xiàn)性和穩(wěn)定性，并未陳述有界弦問題和懸臂梁問題的理論背景，也未開展與現(xiàn)有理論和方法的對比驗(yàn)證研究。即三維諧振子的應(yīng)用前景未得到充分的陳述。因此，在本文研究的基礎(chǔ)上，將在下一步研究中開展三維諧振子的建模應(yīng)用研究，討論網(wǎng)格剖分形式、系統(tǒng)參數(shù)和積分參數(shù)等在有界弦和懸臂梁等定解問題中對數(shù)值解正確性和準(zhǔn)確性的影響，并開展基于三維諧振子的建模方法與現(xiàn)有理論方法的對比驗(yàn)證研究。