具有穩(wěn)定數(shù)值解的三維諧振子

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

1 引言

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

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

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

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

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

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

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

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

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

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

2 三維諧振子

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

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

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

lj=‖rj-r‖

(1)

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

(2)

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

(3)

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

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

(4)

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

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

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

(5)

圖1 量化描述方式

Fig.1 Quantitative description

定常力項(xiàng)與彈簧力項(xiàng)都屬于保守力，合稱為保守力項(xiàng)

F01=F01e01=F0+F1

(6)

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

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

(7)

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

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

(8)

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

(9)

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

e01×Nj=0

(10)

式中0為零向量。此時(shí)，增加一個(gè)隨機(jī)的微量擾動(dòng)力，令

e01×Nj=0+10- neB

(11)

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

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

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

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

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

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

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

(1)v=icsV

(2)v=bcsV[v]

(3)r=icsR

(4)r=bcsR[r]

(5)a=Force[r,v]

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

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

(12)

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

(13)

這個(gè)組合的估算周期值為兩個(gè)周解的最小值

Tj=min.[Ts h,Ts p]

(14)

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

Δt=STmin.{Tj}

(15)

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

迭代過(guò)程包含8個(gè)步驟

(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為一個(gè)位置松弛函數(shù)

Δrb=λbΔr

(16)

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

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

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

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

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

3.1 耗散彈簧擺模型

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

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

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

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

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

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

圖2 初始條件及勢(shì)場(chǎng)

Fig.2 Initial conditions and potential field

表1 定解問(wèn)題

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)動(dòng)軌跡

Fig.3 Trajectory of DSP1

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

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

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

Fig.4 Trajectory of DSP2

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

Fig.5 Trajectory of DSP3

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

Fig.6 Trajectory of DSP4

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

Fig.7 Trajectory of DSP5

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

Fig.8 Trajectory of DSP6

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

Fig.9 Trajectory of DSP7

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

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

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

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

3.2 簡(jiǎn)化弦模型

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

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

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

定常阻尼和彈簧阻尼均設(shè)定為0.5?；谶吔鐥l件和初始條件的設(shè)定，形成三個(gè)定解問(wèn)題DSP8～DSP10。

圖10 DSP1～DSP3能量耗散

Fig.10 Energy dissipation of DSP1～DSP3

圖11 DSP4～DSP7能量耗散

Fig.11 Energy dissipation of DSP4～DSP7

DSP8模擬拉弦動(dòng)作，僅將中間質(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)動(dòng)過(guò)程中最高速度接近1800 m/s，顯然其屬于蠻力型定解條件。

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

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

積分參數(shù)中，松弛系數(shù)為1.0，步長(zhǎng)分率為 1/10，總步數(shù)為2000步，計(jì)算的步長(zhǎng)取兩位有效數(shù)字，均為7.3 μs，模擬總時(shí)長(zhǎng)為14.6 ms。三個(gè)定解問(wèn)題的數(shù)值解依次如圖15～圖17所示。初始時(shí)刻和終止時(shí)刻的弦線均進(jìn)行了顯示。質(zhì)點(diǎn)運(yùn)動(dòng)軌跡投影至xoy和yoz平面上。在多質(zhì)點(diǎn)模型中，質(zhì)點(diǎn)之間的相對(duì)位置關(guān)聯(lián)可以描述弦線的振型。為了在質(zhì)點(diǎn)數(shù)量很少的前提下表示出振型，對(duì)初始時(shí)刻和終止時(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)動(dòng)解(圖15)在位形空間中呈現(xiàn)為有界的近似于圓盤狀的二維空間區(qū)域。從右側(cè)兩質(zhì)點(diǎn)的空心區(qū)域可以看出圓盤面的中心位置質(zhì)點(diǎn)的靜止平衡位置。從xoy投影面可以看出，圓盤面近似垂直于靜止?fàn)顟B(tài)下的弦線。隨著能量的耗散，四個(gè)質(zhì)點(diǎn)的速度逐漸降低，并逐漸向靜止平衡位置靠攏。

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

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

Fig.15 Trajectory of DSP8

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

Fig.16 Trajectory of DSP9

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

Fig.17 Trajectory of DSP10

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

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

3.3 簡(jiǎn)化梁模型

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

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

圖18 DSP8～DSP10能量耗散

Fig.18 Energy dissipation of DSP8～DSP10

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

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

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

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

DSP13模擬梁的扭轉(zhuǎn)行為。將梁的另一個(gè)端部的四個(gè)質(zhì)點(diǎn)的位置以y軸正向?yàn)樾D(zhuǎn)軸按右手螺旋方向旋轉(zhuǎn)45°并固定，作為位置邊界條件。與拉伸和剪切相比，扭轉(zhuǎn)加載本身是幾何非線性的，是對(duì)稱張量難以描述的加載行為。運(yùn)動(dòng)質(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，步長(zhǎng)分率為1/10，總步數(shù)為300步，計(jì)算的步長(zhǎng)取兩位有效數(shù)字后依次為15 μs，16 μs和16 μs。三個(gè)定解問(wèn)題的數(shù)值解依次如圖22～圖24所示。為了顯示變形特征，對(duì)初始時(shí)刻和終止時(shí)刻的簡(jiǎn)分梁的角棱線進(jìn)行了三維曲線擬合，并投影在yox和yoz平面上。

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

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

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

Fig.21 Ics and bcs of DSP13

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

Fig.22 Trajectory of DSP11

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

Fig.23 Trajectory of DSP12

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

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

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

Fig.24 Trajectory of DSP13

圖25 DSP11～DSP13能量耗散

Fig.25 Energy dissipation of DSP11～DSP13

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

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

4 結(jié) 論

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

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

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

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