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

背景技術：

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

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

技術實現要素：

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

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

1)確定積分參數γ、β(t)和步長Δt：

所述積分參數γ設為：γ＝1/2，

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

所述步長Δt滿足：

其中，

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

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

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

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

5)計算t+Δt時刻的有效剛度

6)計算t+Δt時刻的有效外界激振力

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

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

9)使t＝t+Δt；

10)重復步驟1)至9)，計算下一時刻的系統(tǒng)響應，最終求得所需時程范圍內的所有振動響應。

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

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

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

其中，表示變質量動力系統(tǒng)t時刻加速度。

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

其中，為t+Δt時刻的有效剛度。

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

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

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

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

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

附圖說明

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

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

具體實施方式

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

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

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

積分參數γ設為：γ＝1/2，

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

步長Δt滿足：

其中，

其中，式(3)中，

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

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

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

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

步驟五，計算t+Δt時刻的有效剛度計算公式如下：

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

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

其中，表示變質量動力系統(tǒng)t時刻加速度，

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

步驟七，求解t+Δt時刻的位移x_t+Δt，計算公式如下：

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

步驟八，求解t+Δt時刻的速度和加速度計算公式如下：

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

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

步驟十，重復上述步驟一到九，計算下一時刻的系統(tǒng)響應，最終求得所需時程范圍內的所有振動響應。

實施例

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

表1變質量系統(tǒng)參數

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

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

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