本發(fā)明涉及一種用于變質(zhì)量動(dòng)力系統(tǒng)的數(shù)值計(jì)算方法，屬于數(shù)值計(jì)算仿真技術(shù)領(lǐng)域。

背景技術(shù)：

近年來，由于航空航天技術(shù)和工業(yè)技術(shù)的發(fā)展，變質(zhì)量動(dòng)力系統(tǒng)的研究越來越受到研究人員的重視，如火箭、飛機(jī)和車橋耦合振動(dòng)系統(tǒng)等。變質(zhì)量動(dòng)力系統(tǒng)振動(dòng)特性不同于恒定質(zhì)量動(dòng)力系統(tǒng)，其研究也由來已久，早在1897年，俄國(guó)學(xué)者密歇爾斯基就提出了用于研究變質(zhì)量質(zhì)點(diǎn)動(dòng)力學(xué)的密歇爾斯基方程。

單自由度變質(zhì)量動(dòng)力系統(tǒng)，如圖1所示，圖1中m(t)、c、k分別表示系統(tǒng)質(zhì)量、阻尼和剛度，質(zhì)量m(t)是與時(shí)間相關(guān)的函數(shù)；x和f分別表示系統(tǒng)位移和外界激振力。變質(zhì)量動(dòng)力系統(tǒng)屬于非線性動(dòng)力系統(tǒng)，若直接求解動(dòng)力方程比較困難，有時(shí)甚至無法求解，因此研究人員常采用數(shù)值計(jì)算的方法求其振動(dòng)響應(yīng)。中國(guó)發(fā)明專利“變質(zhì)量動(dòng)力吸振器瞬態(tài)過程仿真方法”(公布號(hào)：CN 103851125A)公布了一種用于求解變質(zhì)量動(dòng)力吸振器瞬態(tài)過程的數(shù)值計(jì)算方法，該方法僅考慮到質(zhì)量變化對(duì)系統(tǒng)速度和加速度的影響，但是并沒有考慮到積分參數(shù)和步長(zhǎng)的影響因素，其計(jì)算精度有待提高。發(fā)表在《力學(xué)學(xué)報(bào)》第44卷第5期上的學(xué)術(shù)論文“與結(jié)構(gòu)動(dòng)特性協(xié)同的自適應(yīng)Newmark方法”在計(jì)算變質(zhì)量動(dòng)力系統(tǒng)振動(dòng)響應(yīng)時(shí)，所建立的振動(dòng)方程需考慮質(zhì)量變化率對(duì)阻尼的影響，所建立的數(shù)學(xué)模型較復(fù)雜。

技術(shù)實(shí)現(xiàn)要素：

本發(fā)明所要解決的技術(shù)問題是克服現(xiàn)有技術(shù)的缺陷，提供一種用于變質(zhì)量動(dòng)力系統(tǒng)的數(shù)值計(jì)算方法，可以對(duì)變質(zhì)量動(dòng)力系統(tǒng)的振動(dòng)過程進(jìn)行仿真計(jì)算。

為解決上述技術(shù)問題，本發(fā)明提供一種用于變質(zhì)量動(dòng)力系統(tǒng)的數(shù)值計(jì)算方法，包括以下步驟：

1)確定積分參數(shù)γ、β(t)和步長(zhǎng)Δt：

所述積分參數(shù)γ設(shè)為：γ＝1/2，

所述積分參數(shù)β(t)定義為：

所述步長(zhǎng)Δt滿足：

其中，

k表示變質(zhì)量動(dòng)力系統(tǒng)的剛度，m_t表示變質(zhì)量動(dòng)力系統(tǒng)t時(shí)刻的質(zhì)量，c表示變質(zhì)量動(dòng)力系統(tǒng)的阻尼；

2)計(jì)算積分變量α₀(t)，α₁(t)，α₂(t)，α₃(t)，α₄(t)，α₅(t)；

3)獲得變質(zhì)量動(dòng)力系統(tǒng)t時(shí)刻的狀態(tài)x_t和其中，x_t表示變質(zhì)量動(dòng)力系統(tǒng)t時(shí)刻位移，表示變質(zhì)量動(dòng)力系統(tǒng)t時(shí)刻速度；對(duì)于初始時(shí)刻t為0的變質(zhì)量動(dòng)力系統(tǒng)狀態(tài)為已知量，在求解動(dòng)力響應(yīng)時(shí)，均已知初始時(shí)刻的位移和速度；對(duì)于t≠0時(shí)刻的狀態(tài)，通過t-1時(shí)刻的計(jì)算獲得；

4)獲得t、t+Δt時(shí)刻的變質(zhì)量動(dòng)力系統(tǒng)參數(shù)，包括：t時(shí)刻變質(zhì)量動(dòng)力系統(tǒng)的質(zhì)量m_t，t+Δt時(shí)刻變質(zhì)量動(dòng)力系統(tǒng)的質(zhì)量m_t+Δt，t時(shí)刻變質(zhì)量動(dòng)力系統(tǒng)的外界激振力f_t，t+Δt時(shí)刻變質(zhì)量動(dòng)力系統(tǒng)的外界激振力f_t+Δt，變質(zhì)量動(dòng)力系統(tǒng)的阻尼c和剛度k；

5)計(jì)算t+Δt時(shí)刻的有效剛度

6)計(jì)算t+Δt時(shí)刻的有效外界激振力

7)求解t+Δt時(shí)刻的位移x_t+Δt；

8)求解t+Δt時(shí)刻的速度和加速度

9)使t＝t+Δt；

10)重復(fù)步驟1)至9)，計(jì)算下一時(shí)刻的系統(tǒng)響應(yīng)，最終求得所需時(shí)程范圍內(nèi)的所有振動(dòng)響應(yīng)。

前述的步驟2)中積分變量的計(jì)算公式如下：

前述的步驟5)有效剛度的計(jì)算公式如下：

前述的步驟6)有效外界激振力的計(jì)算公式如下：

其中，表示變質(zhì)量動(dòng)力系統(tǒng)t時(shí)刻加速度。

前述的步驟7)t+Δt時(shí)刻的位移x_t+Δt的計(jì)算公式如下：

其中，為t+Δt時(shí)刻的有效剛度。

前述的步驟8)t+Δt時(shí)刻的速度和加速度的計(jì)算公式如下：

本發(fā)明所達(dá)到的有益效果：

(1)本發(fā)明的數(shù)值計(jì)算方法在計(jì)算變質(zhì)量動(dòng)力系統(tǒng)振動(dòng)響應(yīng)時(shí)，所采用的系統(tǒng)振動(dòng)方程不需考慮質(zhì)量變化率對(duì)阻尼的影響，因此，構(gòu)建變質(zhì)量動(dòng)力系統(tǒng)振動(dòng)方程比較容易。

(2)本發(fā)明在計(jì)算過程中，不但考慮到質(zhì)量對(duì)系統(tǒng)位移、速度和加速度的影響，而且考慮到積分參數(shù)和步長(zhǎng)的影響因素，使算法計(jì)算精度較好。

(3)本發(fā)明的數(shù)值計(jì)算方法適用性廣，可以用于所有的變質(zhì)量動(dòng)力系統(tǒng)。

附圖說明

圖1為變質(zhì)量動(dòng)力系統(tǒng)單自由度模型；

圖2是采用本發(fā)明的數(shù)值計(jì)算方法與自適應(yīng)Newmark算法、傳統(tǒng)Newmark算法所得計(jì)算結(jié)果對(duì)比圖。

具體實(shí)施方式

下面對(duì)本發(fā)明作進(jìn)一步描述。以下實(shí)施例僅用于更加清楚地說明本發(fā)明的技術(shù)方案，而不能以此來限制本發(fā)明的保護(hù)范圍。

本發(fā)明的用于變質(zhì)量動(dòng)力系統(tǒng)的數(shù)值計(jì)算方法，具體包括以下步驟：

步驟一，確定積分參數(shù)γ、β(t)和步長(zhǎng)Δt：

積分參數(shù)γ設(shè)為：γ＝1/2，

積分參數(shù)β(t)定義為：

步長(zhǎng)Δt滿足：

其中，

其中，式(3)中，

k表示變質(zhì)量動(dòng)力系統(tǒng)的剛度，m_t表示變質(zhì)量動(dòng)力系統(tǒng)t時(shí)刻的質(zhì)量，c表示變質(zhì)量動(dòng)力系統(tǒng)的阻尼。

步驟二，計(jì)算積分變量α₀(t)，α₁(t)，α₂(t)，α₃(t)，α₄(t)，α₅(t)：

步驟三，獲得變質(zhì)量動(dòng)力系統(tǒng)t時(shí)刻(初始時(shí)刻t為0)的狀態(tài)(x_t、)，x_t表示變質(zhì)量動(dòng)力系統(tǒng)t時(shí)刻位移，表示變質(zhì)量動(dòng)力系統(tǒng)t時(shí)刻速度；對(duì)于初始時(shí)刻t為0的狀態(tài)為已知量，在求解動(dòng)力響應(yīng)時(shí)，均已知初始時(shí)刻的位移和速度；對(duì)于t≠0時(shí)刻的狀態(tài)，通過t-1時(shí)刻的計(jì)算獲得。

步驟四，獲得t、t+Δt時(shí)刻的變質(zhì)量動(dòng)力系統(tǒng)參數(shù)，包括：t時(shí)刻變質(zhì)量動(dòng)力系統(tǒng)的質(zhì)量m_t，t+Δt時(shí)刻變質(zhì)量動(dòng)力系統(tǒng)的質(zhì)量m_t+Δt，t時(shí)刻變質(zhì)量動(dòng)力系統(tǒng)的外界激振力f_t，t+Δt時(shí)刻變質(zhì)量動(dòng)力系統(tǒng)的外界激振力f_t+Δt，變質(zhì)量動(dòng)力系統(tǒng)的阻尼c和剛度k。

步驟五，計(jì)算t+Δt時(shí)刻的有效剛度計(jì)算公式如下：

上述計(jì)算公式是在Newmark算法的基礎(chǔ)上添加變質(zhì)量元素以及修改積分變量得到的。

步驟六，計(jì)算t+Δt時(shí)刻的有效外界激振力計(jì)算公式如下：

其中，表示變質(zhì)量動(dòng)力系統(tǒng)t時(shí)刻加速度，

上述計(jì)算公式是在Newmark算法的基礎(chǔ)上添加變質(zhì)量元素以及修改積分變量得到的。

步驟七，求解t+Δt時(shí)刻的位移x_t+Δt，計(jì)算公式如下：

上述計(jì)算公式是在Newmark算法的基礎(chǔ)上添加變質(zhì)量元素以及修改積分變量得到的。

步驟八，求解t+Δt時(shí)刻的速度和加速度計(jì)算公式如下：

上述計(jì)算公式是在Newmark算法的基礎(chǔ)上添加變質(zhì)量元素以及修改積分變量得到的。

步驟九，使t＝t+Δt；

步驟十，重復(fù)上述步驟一到九，計(jì)算下一時(shí)刻的系統(tǒng)響應(yīng)，最終求得所需時(shí)程范圍內(nèi)的所有振動(dòng)響應(yīng)。

實(shí)施例

下面通過與自適應(yīng)Newmark算法以及傳統(tǒng)Newmark算法對(duì)比，來驗(yàn)證本發(fā)明所提出的數(shù)值計(jì)算方法。應(yīng)用本發(fā)明數(shù)值計(jì)算方法，對(duì)圖1所示的單自由度變質(zhì)量系統(tǒng)模型進(jìn)行仿真計(jì)算，圖1中系統(tǒng)質(zhì)量m_t、剛度k和阻尼c如表1所示。

表1變質(zhì)量系統(tǒng)參數(shù)

從表1中可以看出，系統(tǒng)質(zhì)量m_t是隨時(shí)間t變化的變量，在初始時(shí)刻m₀為1。系統(tǒng)所受外界激振力f_t為sin10t，初始位移x₀、速度分別取0。

本發(fā)明的數(shù)值計(jì)算方法、自適應(yīng)Newmark算法以及傳統(tǒng)Newmark算法所得計(jì)算結(jié)果如圖2所示。圖2為計(jì)算步長(zhǎng)為0.01s時(shí)變質(zhì)量系統(tǒng)加速度響應(yīng)，從圖中可以看出，本發(fā)明的數(shù)值計(jì)算方法所得結(jié)果與自適應(yīng)Newmark算法、傳統(tǒng)Newmark算法一致，且處于兩算法所得結(jié)果之間，從而說明了本發(fā)明的數(shù)值計(jì)算方法用于計(jì)算變質(zhì)量動(dòng)力吸振器振動(dòng)響應(yīng)的有效性。圖2中，表示本發(fā)明算法，----表示自適應(yīng)Newmark算法，-----表示傳統(tǒng)Newmark算法。

以上所述僅是本發(fā)明的優(yōu)選實(shí)施方式，應(yīng)當(dāng)指出，對(duì)于本技術(shù)領(lǐng)域的普通技術(shù)人員來說，在不脫離本發(fā)明技術(shù)原理的前提下，還可以做出若干改進(jìn)和變形，這些改進(jìn)和變形也應(yīng)視為本發(fā)明的保護(hù)范圍。