(中國科學院 武漢巖土力學研究所 巖土力學與工程國家重點試驗室，武漢 430071)

1 引言

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

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

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

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

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

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

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

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

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

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

2 三維諧振子

2.1 運動的量化描述

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

在任一時刻t，質(zhì)點與相鄰質(zhì)點之間的相互作用僅與距離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 運動微分方程

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

(4)

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

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

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

(5)

圖1 量化描述方式

Fig.1 Quantitative description

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

F01=F01e01=F0+F1

(6)

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

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

(7)

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

F3為彈簧耗散項，是與質(zhì)點關(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為零向量。此時，增加一個隨機的微量擾動力，令

e01×Nj=0+10- neB

(11)

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

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

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

一個物理系統(tǒng)在建模完成后，可以得到系統(tǒng)的物理參數(shù)，包括質(zhì)點的靜態(tài)坐標和質(zhì)量、彈簧的剛度和原長，以及質(zhì)點與彈簧的拓撲關(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ì)點及與其關(guān)聯(lián)的任一個彈簧形成一個組合。設(shè)質(zhì)點的質(zhì)量為m，初始速度的模為v，初始加速度的模為a(a≠0)，彈簧剛度為k，初始長度為l。簡振模式的周期估算公式為

(12)

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

(13)

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

Tj=min.[Ts h,Ts p]

(14)

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

Δt=STmin.{Tj}

(15)

式中ST為步長分率。當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)為蛙跳法的標準流程，步驟(2)和(8)為邊界條件的施加，bcsRRelax為一個位置松弛函數(shù)

Δrb=λbΔr

(16)

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

3 數(shù)值試驗

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

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

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

3.1 耗散彈簧擺模型

三維耗散彈簧擺模型與三維守恒系統(tǒng)彈簧擺模型是一致的，即一個圓柱螺旋拉伸彈簧連接兩個質(zhì)點。質(zhì)點的質(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ì)點固定在參考系的原點，這相當于模型的邊界條件。其位置bcsR與速度bcsV恒為零向量。兩個質(zhì)點的bcsA0設(shè)定為(0,0,-9.81) N/kg。

另一個質(zhì)點設(shè)定四種初始條件，依次編號為ICS1～ICS4。如圖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)平衡位置點為(0,0,-0.625) m。

定常阻尼和彈簧阻尼設(shè)定四種組合，再與初始條件進行二次組合，形成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所示。初始時刻和終止時刻的彈簧擺均進行了顯示。質(zhì)點運動軌跡投影至xoy，xoz和yoz平面上。

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

Fig.3 Trajectory of DSP1

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

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

圖4 DSP2運動軌跡

Fig.4 Trajectory of DSP2

圖5 DSP3運動軌跡

Fig.5 Trajectory of DSP3

圖6 DSP4運動軌跡

Fig.6 Trajectory of DSP4

圖7 DSP5運動軌跡

Fig.7 Trajectory of DSP5

圖8 DSP6運動軌跡

Fig.8 Trajectory of DSP6

圖9 DSP7運動軌跡

Fig.9 Trajectory of DSP7

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

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

三維諧振子的數(shù)值解是穩(wěn)定的。DSP1～DSP2(圖10)和DSP4～DSP7(圖11)6個定解問題，其能量在運動過程中均持續(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的定常阻尼和彈簧阻尼均為零值，其能量在運動過程中依然在減少。DSP3屬于能量守恒系統(tǒng)，其能量減少主要是因為數(shù)值解法中的邊界松弛技術(shù)。在數(shù)值解法中，對不動質(zhì)點(即懸掛支點)的位移松弛，消減了彈簧的彈性勢能，起到了一定的體系消能的作用。

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

3.2 簡化弦模型

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

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

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

定常阻尼和彈簧阻尼均設(shè)定為0.5?；谶吔鐥l件和初始條件的設(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ì)點由(0,2.5,0)拉至(0,2.7,0.1)并作為固定位置邊界條件，所有質(zhì)點的初始位置如圖13所示。icsV，bcsV和bcsA0 均為零向量。初始狀態(tài)下，系統(tǒng)的能量達到4602.5 J；質(zhì)點在運動過程中最高速度接近1800 m/s，顯然其屬于蠻力型定解條件。

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

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

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

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

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

圖15 DSP8運動軌跡

Fig.15 Trajectory of DSP8

圖16 DSP9運動軌跡

Fig.16 Trajectory of DSP9

圖17 DSP10運動軌跡

Fig.17 Trajectory of DSP10

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

由三維諧振子描述的簡化弦模型的數(shù)值解是穩(wěn)定的。DSP8～DSP10(圖18)的能量在運動過程中均持續(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 簡化梁模型

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

梁的建模采用正方形截面鋼材。截面尺寸為10 cm×10 cm，長度50 cm，材料密度為7.9 g/cm3，材料彈模為200 GPa。將梁離散成四面體單元之后，為了得到三維諧振子的物理參數(shù)，將單元體的體積分配至與單元體關(guān)聯(lián)的棱，再分配至與單元體關(guān)聯(liá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)的頂點；最后，每個頂點和每條棱分配得到體積進行求和。這樣，每一條棱和每一個頂點都有了相應(yīng)的體積量。

圖18 DSP8～DSP10能量耗散

Fig.18 Energy dissipation of DSP8～DSP10

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

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

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

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

DSP13模擬梁的扭轉(zhuǎn)行為。將梁的另一個端部的四個質(zhì)點的位置以y軸正向為旋轉(zhuǎn)軸按右手螺旋方向旋轉(zhuǎn)45°并固定，作為位置邊界條件。與拉伸和剪切相比，扭轉(zhuǎn)加載本身是幾何非線性的，是對稱張量難以描述的加載行為。運動質(zhì)點的初始位置如圖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所示。為了顯示變形特征，對初始時刻和終止時刻的簡分梁的角棱線進行了三維曲線擬合，并投影在yox和yoz平面上。

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

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

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

Fig.21 Ics and bcs of DSP13

圖22 DSP11運動軌跡

Fig.22 Trajectory of DSP11

圖23 DSP12運動軌跡

Fig.23 Trajectory of DSP12

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

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

圖24 DSP13運動軌跡

Fig.24 Trajectory of DSP13

圖25 DSP11～DSP13能量耗散

Fig.25 Energy dissipation of DSP11～DSP13

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

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

4 結(jié) 論

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

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

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

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