BAT子彈藥折疊翼展開沖擊非光滑模型的隱式積分算法

張洪銘，顧曉輝，孫 麗

(1.南京理工大學(xué) 機(jī)械工程學(xué)院，南京 210094；2.江蘇科技大學(xué) 機(jī)械工程學(xué)院，江蘇 鎮(zhèn)江 212100)

折疊翼結(jié)構(gòu)能有效縮小中近程導(dǎo)彈、制導(dǎo)炮彈的儲運(yùn)空間，利于增強(qiáng)車輛、戰(zhàn)機(jī)的運(yùn)載能力，提高武器系統(tǒng)的戰(zhàn)術(shù)性能，在現(xiàn)代飛行器的設(shè)計中獲得了廣泛的應(yīng)用[1]。不同于固定翼結(jié)構(gòu)近似剛性的連接，折疊翼展開機(jī)構(gòu)中普遍包含鉸鏈、限位擋塊、鎖定插銷及驅(qū)動裝置?；顒蛹拇嬖诮档土私M件間的連接強(qiáng)度和系統(tǒng)的結(jié)構(gòu)剛度，在迅速展開過程中，驅(qū)動加載的瞬間以及翼面展開到位時刻的碰撞、鎖定動作，均會在部件間產(chǎn)生沖擊載荷，存在使構(gòu)件被破壞的可能性。另外，翼面的展弦比較大時，展開末段的碰撞沖擊會引起翼面振顫，不恰當(dāng)?shù)脑O(shè)計可能引發(fā)翼面與彈體共振，加劇結(jié)構(gòu)破壞的程度或影響飛行器的運(yùn)動軌跡。

因此，針對折疊翼展開的動態(tài)過程，沖擊載荷特性的分析與優(yōu)化是提升展開性能的關(guān)鍵?，F(xiàn)有研究對折疊翼展開過程進(jìn)行了廣泛研究。林三春等[2]結(jié)合隱式動力學(xué)仿真分析和試驗(yàn)研究的結(jié)果，指出柵格翼展開沖擊的主要影響因素是阻尼力矩與展開到位時的角速度，氣動力和過載為次要因素；唐霄漢等[3]詳細(xì)討論了折疊翼機(jī)構(gòu)中常見的鉸鏈關(guān)節(jié)在沖擊過程中的內(nèi)力分析，討論了材料采用動態(tài)本構(gòu)模型時，應(yīng)變率對等效塑性應(yīng)變的影響；甄文強(qiáng)等在文獻(xiàn)[4]中建立了考慮氣動和摩擦阻力的數(shù)學(xué)模型，并與動力學(xué)仿真結(jié)果進(jìn)行了對比分析；更進(jìn)一步地，其在文獻(xiàn)[5]中進(jìn)行了發(fā)射環(huán)境下折疊翼展開過程的研究，為折疊翼展開機(jī)構(gòu)的工程設(shè)計提供了參考。上述研究對折疊翼展開過程中的沖擊載荷、動力學(xué)響應(yīng)進(jìn)行了分析，關(guān)于沖擊載荷的優(yōu)化分析，李莉[6]、蔡德詠[7]等分別建立了縱向和橫向折疊翼的剛體模型，并均以展開到位瞬時的翼片動能為優(yōu)化目標(biāo)，對展開系統(tǒng)進(jìn)行了優(yōu)化。趙俊鋒等[8]以剛?cè)狁詈夏Ｐ蜑榛A(chǔ)，采用序列二次規(guī)劃法(sequential quadratic programming,SQP)對驅(qū)動力-時間曲線進(jìn)行優(yōu)化以降低彈翼展開沖擊；韓雪峰等[9]則采用正交試驗(yàn)對折疊翼展開機(jī)構(gòu)進(jìn)行優(yōu)化。郭宏偉等[10]在提出沖擊加速度、力矩計算方法的同時，給出了兩種降低展開沖擊的設(shè)計方案。

現(xiàn)有研究較多關(guān)注折疊翼展開到位時，翼面對沖擊載荷的響應(yīng)模擬，對于接觸/碰撞沖擊進(jìn)行建模、計算的內(nèi)容相對較少。實(shí)際上，準(zhǔn)確計算折疊翼展開過程中的沖擊力，是對展開機(jī)構(gòu)進(jìn)行優(yōu)化設(shè)計、可靠性分析的基礎(chǔ)和前提。折疊翼展開過程中包含的接觸/碰撞現(xiàn)象，屬于典型的非光滑、非線性過程，其動力學(xué)建模和求解相對困難和復(fù)雜。目前，描述接觸/碰撞過程的動力學(xué)模型主要有兩類[11-12]：基于接觸體幾何約束的非光滑模型[13-14]和基于接觸界面變形量的連續(xù)接觸力模型[15-17](或稱為罰函數(shù)模型[18])。在包含非光滑過程的動力學(xué)模型中，由不等式來描述系統(tǒng)中的單邊約束，對不等式的處理是建立模型的難點(diǎn)。針對該難點(diǎn)，非光滑模型通過微分包含理論，對接觸/碰撞、摩擦等單邊約束進(jìn)行表達(dá)和處理；連續(xù)接觸力模型是另一種被廣泛研究的接觸/碰撞動力學(xué)模型，采用接觸力約束代替幾何約束，將接觸力模型引入到動力學(xué)方程中進(jìn)行計算。相較于連續(xù)接觸力模型，非光滑模型采用局部剛性接觸避免了計算剛性微分方程的困難，計算精度較高[19]，可以處理碰撞過程中能量不協(xié)調(diào)的問題[20]。另外，基于線性互補(bǔ)問題(linear complementarity problem，LCP)的非光滑動力學(xué)模型，不僅使得力學(xué)問題具有更加自然、精確的數(shù)學(xué)描述，而且可以運(yùn)用互補(bǔ)問題解的存在性和唯一性理論分析力學(xué)問題的結(jié)構(gòu)模型，構(gòu)建更加有效和強(qiáng)健的算法。同時，互補(bǔ)形式的非光滑動力學(xué)模型能夠自動判斷觸-離、滯-滑等狀態(tài)的切換，避免了進(jìn)行事件檢測時所需的大量計算資源[21]。

基于上述研究結(jié)果，以某型智能反坦克(brainpower anti-tank,BAT)子彈藥的折疊翼機(jī)構(gòu)為研究背景，建立該縱向折疊翼展開過程的非光滑動力學(xué)模型，并推導(dǎo)模型的線性互補(bǔ)形式；文中提出了一種新型隱式積分算法用于求解線性互補(bǔ)形式的非光滑動力學(xué)模型，并使用基于所提隱式積分算法的時間步進(jìn)算法，來計算翼面展開到位時碰撞沖擊力量值的大小。將模型計算結(jié)果與數(shù)值仿真結(jié)果進(jìn)行了對比，驗(yàn)證文中非光滑模型的計算精度與求解效率。結(jié)果表明，線性互補(bǔ)形式的非光滑動力學(xué)模型能夠有效描述接觸/碰撞過程中的單邊約束，并且文中提出的隱式積分算法相較于Moreau中點(diǎn)法具有更高的計算精度。

1 折疊翼展開過程的非光滑動力學(xué)模型

1.1 折疊翼構(gòu)型及展開特性

某型BAT子彈藥外形如圖1所示，由運(yùn)載平臺釋放后展開折疊翼和尾翼，配合紅外/聲傳感器復(fù)合導(dǎo)引，滑翔至目標(biāo)上空實(shí)施攻擊。彈體中部裝有四片與彈體縱軸垂直、展開后呈十字型分布的折疊式彈翼。

圖1 某型BAT子彈藥外形圖Fig.1 Shape diagram of BAT submunitions

本文的研究內(nèi)容主要關(guān)注縱向彈載折疊翼展開過程中，碰撞沖擊的建模與計算，因此，忽略聲傳感器支撐桿部分。單組折疊翼組件的裝配示意如圖2所示，展開系統(tǒng)由彈體基座、限位擋塊、轉(zhuǎn)軸、翼面組件、連桿、滑塊以及驅(qū)動彈簧共同組成。圖2中，Ox為彈軸方向，Oy為彈徑方向，A、B、C為三處鉸接點(diǎn)。BAT子彈藥折疊翼展開動作過程為：被運(yùn)載平臺拋散釋放后，在展開過程的初始時刻，翼面處于收攏狀態(tài)，圖2中φ1為0°；隨后滑塊在壓縮彈簧的驅(qū)動下沿滑軌向彈頭方向運(yùn)動，通過連桿對翼面組件施加轉(zhuǎn)動力矩，使得翼面繞鉸接點(diǎn)A逆時針轉(zhuǎn)動；翼面展開到90°位時與限位擋塊之間發(fā)生碰撞，同時鎖定機(jī)構(gòu)啟動，對展開到位的翼面進(jìn)行鎖定防止回轉(zhuǎn)，展開動作結(jié)束。

圖2 單組折疊翼展開機(jī)構(gòu)Fig.2 Schematic diagram of folding wing deployment mechanism

依據(jù)設(shè)計指標(biāo)，折疊翼的展開動作需要在0.2 s內(nèi)完成，即翼面組件需要在短時間內(nèi)完成大范圍轉(zhuǎn)動，轉(zhuǎn)動速度高、動能大，在展開動作末段翼面與限位擋塊、鎖定銷之間會發(fā)生相對劇烈的碰撞沖擊。沖擊過程持續(xù)時間短、作用力大，使得包含碰撞沖擊的折疊翼的展開動作是一個典型的非光滑動力學(xué)過程。準(zhǔn)確計算碰撞力量值對展開機(jī)構(gòu)優(yōu)化、機(jī)構(gòu)強(qiáng)度校核均具有重要意義。后續(xù)內(nèi)容建立能夠準(zhǔn)確描述碰撞過程的非光滑動力學(xué)模型，和用于解算非光滑模型數(shù)值計算方法，用于計算翼面與限位擋塊首次碰撞力的大小。

1.2 折疊翼非光滑動力學(xué)建模

對圖2所示折疊翼多體系統(tǒng)，選取最大坐標(biāo)數(shù)目，在絕對笛卡爾坐標(biāo)下描述各部件的位形。在未發(fā)生接觸/碰撞現(xiàn)象的自由運(yùn)動階段，系統(tǒng)內(nèi)僅有光滑的等式雙邊約束，基于拉格朗日乘子法將約束方程與運(yùn)動微分方程結(jié)合，對含有n個自由度、m個雙邊約束(m

(1)

(2)

(3)

通常取α=β∈[5,50]系統(tǒng)響應(yīng)最快，在每個時間步內(nèi)，系統(tǒng)的位置約束與速度約束同時被滿足，使得積分誤差得到部分抑制。式(3)即為光滑多體系統(tǒng)的完備動力學(xué)方程。

多體系統(tǒng)內(nèi)部構(gòu)件間發(fā)生接觸/碰撞時，系統(tǒng)拓?fù)錁?gòu)型發(fā)生改變，除原有雙邊約束外，系統(tǒng)內(nèi)增加了額外的單邊約束使構(gòu)件運(yùn)動狀態(tài)發(fā)生突變。對于多體系統(tǒng)內(nèi)所有潛在發(fā)生正、斜碰撞的接觸面，借助Lagrange乘子法向多體系統(tǒng)的運(yùn)動微分方程中添加法向、切向約束反力得到

(4)

式中：λN、λT分別為未知的法向碰撞力和切向摩擦力；WN與WT分別為法向、切向距離約束函數(shù)對廣義坐標(biāo)q的Jacobian矩陣的轉(zhuǎn)置

(5)

(6)

式中,gN(q,t)、gT(q,t)分別為接觸面法向、切向距離函數(shù)列陣。

由于碰撞過程發(fā)生時間短、作用力大，碰撞前后速度發(fā)生突變，考慮將式(6)中力-加速度關(guān)系在碰撞過程時間Δt上積分，寫成速度-沖量形式。在碰撞起始時刻ta至結(jié)束時刻tb上，式(4)兩側(cè)同時關(guān)于時間積分得到

(7)

(8)

式中，ΛN和ΛT分別為法向、切向上的單邊約束力沖量。將系統(tǒng)的幾何約束方程對時間微分，代入tb時刻的廣義速度，得到目標(biāo)構(gòu)件在碰撞結(jié)束時刻的速度約束方程

(9)

將碰撞結(jié)束時刻接觸面間法向、切向上的相對速度用廣義速度表示為

(10)

(11)

非光滑動力學(xué)角度下的互補(bǔ)理論是對接觸/碰撞等非線性過程中物理量所遵從規(guī)律的描述。基于互補(bǔ)理論，對系統(tǒng)內(nèi)滿足gN(q,t)=0條件的接觸對，法向碰撞力沖量ΛN和切向摩擦余量的沖量ΛTS，分別與各自對應(yīng)方向上，碰撞前后相對速度的合成量構(gòu)成互補(bǔ)關(guān)系

0≤ΛN⊥ξN≥0

(12)

0≤ΛTS⊥ξT≥0

(13)

式中：符號⊥表示兩個向量垂直；切向摩擦余量的沖量ΛTS=μΛN-ΛT；ξN、ξT分別為法向、切向合成速度

(14)

(15)

式中,eN、eT分別為法向、切向上的恢復(fù)系數(shù)對角矩陣。

由式(7)、式(9)、式(12)～式(15)共6個方程共同構(gòu)成了基于互補(bǔ)接觸定律的非光滑動力學(xué)模型

(16)

2 非光滑模型的數(shù)值計算方法

形如式(16)的微分代數(shù)方程組(differential algebraic equations,DAEs)從速度-沖量角度描述了多剛體系統(tǒng)內(nèi)部構(gòu)件之間發(fā)生接觸/碰撞現(xiàn)象時，多體系統(tǒng)的運(yùn)行規(guī)律。式(16)闡明了單邊約束力沖量與合成速度之間存在互補(bǔ)關(guān)系，但未顯式給出兩者間的數(shù)量關(guān)系，無法直接求解。本章推導(dǎo)式(16)的LCP形式以求解非光滑過程，求解主要涉及兩部分內(nèi)容，即成對互補(bǔ)量的計算和DAEs的迭代求解。通過確定相鄰時間點(diǎn)上廣義速度和廣義坐標(biāo)的隱式積分算法，從而對形如式(16)的DAEs采用時間步進(jìn)算法進(jìn)行求解。

2.1 非光滑模型的線性互補(bǔ)形式

采用庫倫模型描述接觸面件的動/靜摩擦力，包含庫倫摩擦模型的非光滑動力學(xué)模型是一個微分包含[22]，在切向相對速度為零時，摩擦力是一個多值函數(shù)。為了便于判斷接觸狀態(tài)的轉(zhuǎn)換，對于法向和切向上相對距離、相對速度為零的接觸對，將摩擦余量的沖量分解為左、右分量

ΛL=μΛN-ΛT

(17)

ΛR=μΛN+ΛT

(18)

則

ΛR=2μΛN-ΛL

(19)

相應(yīng)地將切向合成速度分解為左、右分量之差

ξT=ξR-ξL

(20)

將式(17)代入式(16)第一式、式(20)代入式(16)第四式，聯(lián)立式(14)、式(19)后寫成矩陣形式有

(21)

將式(21)等號兩側(cè)內(nèi)容對換并按分塊矩陣簡寫為

(22)

(23)

將式(23)代入式(22)第二個方程，得到Λ的表達(dá)式

(24)

將式(24)代入式(22)第三個方程中，得到互補(bǔ)量y和x之間的標(biāo)準(zhǔn)線性互補(bǔ)形式

(25)

式中：

2.2 非光滑模型的隱式積分算法

(26)

基于上述討論，構(gòu)造非光滑動力學(xué)模型在速度-沖量層面的時間步進(jìn)算法，求解描述折疊翼展開系統(tǒng)中含有接觸/碰撞過程的微分代數(shù)方程組。時間步進(jìn)算法的流程如圖3所示，已知展開系統(tǒng)在t時刻的狀態(tài)，判斷該時刻是否有接觸對在切向、法向上的距離為零：若非零，則采用Newmark隱式迭代算法求解系統(tǒng)平穩(wěn)運(yùn)行至下一時刻的坐標(biāo)和速度；若單邊約束為零，表明系統(tǒng)內(nèi)部構(gòu)件間發(fā)生接觸/碰撞，動力學(xué)方程組由式(5)變?yōu)槭?18)，解出一對互補(bǔ)的單邊約束力沖量與合成速度后，按式(23)中的迭代算法計算非光滑過程下一時刻的廣義坐標(biāo)和速度，而后重新判斷單邊距離約束是否為零，從而進(jìn)入迭代，直至總計算時間的終點(diǎn)時刻tF。

3 BAT子彈藥折疊翼展開沖擊分析

前兩章內(nèi)容建立了縱向折疊翼展開過程的非光滑動力學(xué)模型，并構(gòu)造了非光滑模型的線性互補(bǔ)形式，該形式使用文中提出的隱式積分算法進(jìn)行求解時，相較于Moreau中點(diǎn)法能夠獲得更加精確的碰撞力值。本章通過一個工程算例，驗(yàn)證基于隱式積分算法的時間步進(jìn)算法對非光滑動力學(xué)模型進(jìn)行數(shù)值計算的精度和計算效率，并對計算結(jié)果進(jìn)行分析。

某型BAT子彈藥折疊翼展開動作原理如圖4所示。翼面解鎖后展開的主要動力來源于滑塊底部的彈簧組件，處于壓縮狀態(tài)的彈簧推動滑塊沿滑桿向前運(yùn)動，滑塊帶動連桿推動翼面繞鉸點(diǎn)A轉(zhuǎn)動，翼面與限位擋塊碰撞后展開動作完成。在折疊翼從收攏至展開的過程中，驅(qū)動彈簧始終處于壓縮狀態(tài)。圖4系統(tǒng)中各量的取值如表1所示。表1中：K為驅(qū)動彈簧的剛度系數(shù)；L0和L分別為驅(qū)動彈簧的初始長度和壓縮剩余長度；H為0.075 m。圖4展開機(jī)構(gòu)中滑移鉸和柱鉸均處于良好潤滑條件，摩擦力量值相較于碰撞力很小，故而不考慮系統(tǒng)中摩擦力的影響。翼面以及擋塊使用的材料為7A04-T6鋁合金，硬度HB150。

圖4 展開機(jī)構(gòu)原理Fig.4 Schematic diagram of folding wing expansion principle

表1 BAT子彈藥折疊翼機(jī)構(gòu)參數(shù)Tab.1 Parameters of folding wing mechanism of BAT submunitions

圖4中彈簧的驅(qū)動力-時間關(guān)系如圖5所示。力的量值為負(fù)表示彈簧處于壓縮狀態(tài)，為正表示拉伸狀態(tài)。

圖5 彈簧驅(qū)動力-時間曲線Fig.5 Spring force-time curve

對圖4中的折疊翼展開系統(tǒng)，建立與彈體固連的慣性笛卡爾坐標(biāo)系Oxy，在圖示二維平面內(nèi)選取最大坐標(biāo)數(shù)目，令q=[x1,y1,θ,x2,y2,φ,x3,y3,γ]T為廣義坐標(biāo)，建立折疊翼展開過程的非光滑動力學(xué)模型。折疊翼展開系統(tǒng)的雙邊約束函數(shù)為

由于該系統(tǒng)為定常系統(tǒng)，雙邊約束函數(shù)中不顯含時間，則可將系統(tǒng)的雙邊約束方程寫成Φ(q)=0；在翼面展開到位并與限位擋塊之間發(fā)生碰撞階段的單邊約束函數(shù)為

將對應(yīng)項代入式(21)～式(22)中，得到折疊翼展開過程的非光滑動力學(xué)模型。使用文中提出的隱式積分算法，取碰撞過程的時間步Δt=0.005 s，解算從收攏至展開到位的完整過程中，翼面自由端的角位移和角速度在不同時刻的值，結(jié)果如圖6所示。

圖6 翼面角位移-時間曲線和角速度-時間曲線Fig.6 Angular displacement-time curve and angular velocity-time curve of folded wing

為了驗(yàn)證2.2節(jié)中基于隱式積分算法的數(shù)值算法的計算精度與穩(wěn)定性，分別考察文中所提算法的計算結(jié)果與連續(xù)力模型計算結(jié)果的吻合程度，以及計算結(jié)果對幾何約束的相對誤差。采用與數(shù)值計算相同的參數(shù)，使用連續(xù)接觸力模型對折疊翼展開過程進(jìn)行計算(連續(xù)力模型的計算在RecurDyn?軟件中實(shí)現(xiàn))。將不同展開時間下，兩種模型計算得到的翼面角位移、角速度分別進(jìn)行對比，結(jié)果如圖7所示。

圖7中曲線對比結(jié)果表明：數(shù)值計算結(jié)果與仿真結(jié)果吻合；首次碰撞發(fā)生后，基于隱式積分算法的數(shù)值計算方法所得出的最大反彈角位移與反彈角速度，略微大于連續(xù)接觸力模型的結(jié)果。原因在于：連續(xù)力模型假設(shè)物體界面在相互接觸/碰撞的過程中可以發(fā)生微小的穿透，且碰撞力由接觸力和阻尼力兩部分組成，兩者均為界面相互侵入量的函數(shù)。連續(xù)力模型中的阻尼力在接觸/碰撞過程中起耗能作用，因此使用罰函數(shù)模型算得的仿真結(jié)果略小于數(shù)值計算的結(jié)果。

(a) 角位移對比

圖4所示系統(tǒng)在展開過程中，翼面的質(zhì)心與鉸點(diǎn)A之間的距離理論上應(yīng)該為固定值，但由于數(shù)值計算存在積分誤差，使得該固定值的計算結(jié)果在理論值附近有微小偏移。采用該固定值數(shù)值計算結(jié)果與理論值之間的相對誤差，來度量文中所提算法的精度具有較高的合理性。基于隱式積分算法的數(shù)值計算結(jié)果與理論值之間的相對誤差如圖8所示。相對誤差值的數(shù)量級在10-8量級，且整個展開過程中誤差均值穩(wěn)定并且分布均勻，表明文中所提的基于隱式積分算法求解的非光滑動力學(xué)模型能夠準(zhǔn)確、穩(wěn)定地描述含有碰撞過程的折疊翼展開動作。

圖8 文中所提方法求解結(jié)果的相對誤差Fig.8 The relative error between the solution result of the proposed method and the exact solution

基于隱式積分算法的時間步進(jìn)算法能夠更準(zhǔn)確地計算接觸力。使用文中提出方法、Moreau中點(diǎn)算法、連續(xù)力模型算法，分別對圖4中的折疊翼展開過程進(jìn)行計算，依據(jù)力與沖量之間的關(guān)系，將文中提出的隱式積分法和Moreau中點(diǎn)算法的碰撞沖量結(jié)果換算成碰撞力，三種算法在所有時間步上的碰撞力隨展開時間的變化如圖9中曲線所示。觀察圖9中的放大圖，比較三種方法算得的首次碰撞力的大小，文中提出的隱式積分算法結(jié)果與RecurDyn計算結(jié)果吻合，Moreau中點(diǎn)算法的計算結(jié)果偏大；將文中提出的隱式積分法計算結(jié)果與連續(xù)接觸力模型的計算結(jié)果(由RecurDyn實(shí)現(xiàn))相比較，結(jié)果表明：隱式積分法基于非光滑動力學(xué)模型，碰撞過程中接觸對表面之間不存在相互侵入過程，高速、低載荷情形下計算結(jié)果更接近真實(shí)物理過程，得到的碰撞力量值更加準(zhǔn)確。

圖9 三種方法對折疊翼展開過程碰撞力的計算結(jié)果Fig.9 The calculation results of impact forces of folding wings by three methods

4 結(jié) 論

對含有接觸非線性的折疊翼展開系統(tǒng)內(nèi)部構(gòu)件間存在的接觸/碰撞過程，在現(xiàn)有的速度-沖量形式非光滑動力學(xué)模型基礎(chǔ)上，推導(dǎo)了模型的線性互補(bǔ)形式。基于線性互補(bǔ)形式的非光滑動力學(xué)模型，借助互補(bǔ)原理能夠以較少的計算量檢測和判斷觸-離、滯-滑狀態(tài)的切換，節(jié)省事件檢測所需的計算資源。

文中提出了新的隱式積分算法，用于迭代求解線性互補(bǔ)形式的非光滑動力學(xué)模型，并基于該隱式積分算法構(gòu)造了基于速度和沖量的時間步進(jìn)算法。文中所提積分算法基于隱式歐拉法，用前一時刻的真實(shí)速度和后一時刻的近似速度，加權(quán)構(gòu)造后一時刻的位移，更新后一時刻各狀態(tài)量，算出后一時刻的真實(shí)速度。相較于Moreau中點(diǎn)算法，文中提出的隱式積分算法具有更高的計算精度，使得基于該算法的時間步進(jìn)算法在求解包含單點(diǎn)碰撞的非光滑過程時，具有更高的效率。文末算例檢驗(yàn)了集成新型積分算法的時間步進(jìn)算法的精度與穩(wěn)定性。