控制棒下落與流體流動(dòng)的耦合狀態(tài)方程及其保辛算法*

趙 珂，陳昌義，席炎炎，黃東威，吳 鋒，鐘萬勰

（1.大連理工大學(xué) 工程力學(xué)系 工業(yè)裝備結(jié)構(gòu)分析國家重點(diǎn)實(shí)驗(yàn)室，遼寧 大連 116024；2.中廣核研究院有限公司，廣東 深圳 518000）

（我刊編委鐘萬勰來稿）

引 言

控制棒是核反應(yīng)堆在發(fā)生緊急情況時(shí)控制停堆的重要部件，控制棒的下落時(shí)間是保證核反應(yīng)安全的關(guān)鍵指標(biāo).近年來，國內(nèi)外許多學(xué)者對控制棒下落時(shí)間的問題做了大量研究：文獻(xiàn)[1]通過商業(yè)軟件ADINA 建立了單根控制棒二維流-固仿真模型，并將ADINA 得到的計(jì)算結(jié)果與未公開的企業(yè)內(nèi)部代碼得到的結(jié)果進(jìn)行比較，落棒時(shí)間誤差在7%以內(nèi)；文獻(xiàn)[2]基于FLUENT 動(dòng)網(wǎng)格算法，建立了單根控制棒三維流體仿真模型，將得到的控制棒落棒時(shí)間與試驗(yàn)結(jié)果進(jìn)行比較，誤差小于6%.采用二維、三維動(dòng)力學(xué)模型的缺點(diǎn)是自由度多、計(jì)算量大，因此計(jì)算效率較低，往往只能分析單根控制棒，不能分析整個(gè)控制棒組件.考慮到控制棒和導(dǎo)向管均具有長細(xì)比大的特點(diǎn)，且控制棒的直徑與導(dǎo)向管的直徑相近，文獻(xiàn)[3]基于一維Navier-Stokes 方程建立了控制棒下落的數(shù)學(xué)模型，得到的結(jié)果能夠很好地匹配試驗(yàn)結(jié)果；文獻(xiàn)[4-5]構(gòu)建了基于一維Navier-Stokes 方程的落棒動(dòng)力學(xué)模型，通過編譯C++程序進(jìn)行計(jì)算.這些研究中常將流體與固體分開進(jìn)行分析[4-5]：首先，假定控制棒的運(yùn)動(dòng)狀態(tài)恒定不變，計(jì)算流體流動(dòng)的流量和壓強(qiáng)分布，得到流體阻力；然后，假定流體為擬穩(wěn)態(tài)流動(dòng)，其速度和壓強(qiáng)狀態(tài)不變，根據(jù)流體阻力計(jì)算控制棒的速度等變量.這種算法本質(zhì)上是動(dòng)力學(xué)求解算法中的算子分裂法[6]，因此只有一階精度，計(jì)算精度較低，需要取非常小的時(shí)間步長[3]才能保證“擬穩(wěn)態(tài)流動(dòng)”近似的精度.此外，控制棒下落過程中，流場中的流速、壓強(qiáng)等會(huì)發(fā)生非線性突變，采用單一的時(shí)間步長，算子分裂算法還存在不收斂問題，本文的研究將致力于解決這些問題.

在控制棒下落過程中，控制棒的下落與流體流動(dòng)是同時(shí)發(fā)生的，因此應(yīng)該采用流固耦合方法，以真實(shí)地描述控制棒的下落過程.基于流固耦合思想，本文首先建立了控制棒運(yùn)動(dòng)方程和導(dǎo)向管內(nèi)流體壓降平衡的非線性微分方程組；然后引入關(guān)鍵恒等變換=v，得到非線性微分方程組的狀態(tài)方程.為精確求解該非線性狀態(tài)方程，本文建立了一種時(shí)間步長自適應(yīng)保辛算法.該方法可以根據(jù)流動(dòng)的實(shí)際狀態(tài)實(shí)時(shí)調(diào)整時(shí)間步長，精確捕捉流動(dòng)狀態(tài)的突變，避免數(shù)值發(fā)散，并可以在流動(dòng)變化平穩(wěn)的時(shí)間段采用較大的時(shí)間步長，提高計(jì)算效率.

1 理論模型

1.1 物理模型

控制棒和導(dǎo)向管的縱剖面圖如圖1所示，其中陰影部分表示控制棒，控制棒下面是導(dǎo)向管，控制棒的直徑小于導(dǎo)向管的直徑，且控制棒的長度大于導(dǎo)向管的長度.導(dǎo)向管和控制棒全部豎直放置，導(dǎo)向管管壁和底部中心處均有流水孔，分別記為側(cè)孔和底部流水孔.導(dǎo)向管和控制棒完全浸在水中，底部流水孔保持進(jìn)水狀態(tài)，而側(cè)孔存在流進(jìn)和流出兩種狀態(tài)的突變.在控制棒下落過程中，導(dǎo)向管固定不動(dòng).初始時(shí)刻，控制棒在導(dǎo)向管側(cè)孔上方，如圖1(a)所示.控制棒由靜止釋放后，在重力、流體作用力等不同力的作用下開始下落，經(jīng)過一段時(shí)間，控制棒底部下落至導(dǎo)向管側(cè)孔下方，如圖1(b)所示.

對圖1中的結(jié)構(gòu)建立坐標(biāo).以導(dǎo)向管底部截面圓心為坐標(biāo)原點(diǎn)o，以導(dǎo)向管中心軸為z軸，向上為正.圖1中，導(dǎo)向管底部Z0的坐標(biāo)為0，對應(yīng)的壓強(qiáng)為PZ0；導(dǎo)向管側(cè)孔的坐標(biāo)為Z1，對應(yīng)的壓強(qiáng)為PZ1；控制棒底部的坐標(biāo)為Z，對應(yīng)的壓強(qiáng)為PZ，Z在下落過程中會(huì)隨時(shí)間而變化；導(dǎo)向管頂部的坐標(biāo)為Z2，對應(yīng)的壓強(qiáng)為PZ2.

圖1(a)中，控制棒底部在導(dǎo)向管側(cè)孔上方.根據(jù)導(dǎo)向管內(nèi)流量的不同可將導(dǎo)向管分成三個(gè)區(qū)間，分別為[Z,Z2]，[Z1,Z]，[Z0,Z1]，對應(yīng)的流量分別記為Q2，Q1_Z，Q0.從導(dǎo)向管側(cè)孔流進(jìn)或流出的流量記為Q1，以流出為正.同樣地，圖1(b)中也根據(jù)導(dǎo)向管內(nèi)流量的不同，將導(dǎo)向管分成三個(gè)區(qū)間，分別為[Z1,Z2]，[Z,Z1]，[Z0,Z]，對應(yīng)的流量記為Q2，Q1_Z，Q0，而導(dǎo)向管側(cè)孔處的流量仍記為Q1，以流出為正.

圖1 控制棒在導(dǎo)向管內(nèi)的下落示意圖Fig.1 The falling diagram of the control rod in the guide tube

為方便論述，導(dǎo)向管頂部以上的控制棒部分所受的力暫不考慮；且假定控制棒在下落過程中保持豎直狀態(tài)，即控制棒沒有橫向偏移和橫向振動(dòng)；反應(yīng)堆內(nèi)的流體滿足不可壓縮條件.

1.2 控制棒的運(yùn)動(dòng)方程

控制棒的運(yùn)動(dòng)受Newton 第二定律控制，其控制方程可以寫成

其中，M是控制棒的質(zhì)量；w是控制棒的位移；表示控制棒的加速度；F是控制棒受到的合力，可表示為

式(2)中，F(xiàn)g為控制棒的重力；Fb為控制棒受到的浮力；frod為控制棒受到的流體摩擦阻力；Fr為控制棒受到的流體壓差阻力；Fc為控制棒受到的機(jī)械摩擦力，由實(shí)驗(yàn)擬合得到；Fs為控制棒受到的彈簧阻力，此力在控制棒快要到導(dǎo)向管底部時(shí)產(chǎn)生，防止由于控制棒速度過大從而與導(dǎo)向管產(chǎn)生強(qiáng)烈碰撞，可表示為

其中，F(xiàn)prl表示彈簧預(yù)緊力，Ks表示彈簧剛度，hs表示控制棒接觸彈簧時(shí)的位移.

在某區(qū)間段[Za,Zb]內(nèi)，摩擦力frod可表示為[7]

其中，ρ為流體密度，vf為流體流速，v為控制棒的下落速度，C2為控制棒的摩擦因數(shù)，l2為控制棒的濕周長，sgn為符號(hào)函數(shù).

在式(2)所示控制棒受到的幾種合力中，重力和浮力很容易計(jì)算，摩擦阻力可通過式(4)來計(jì)算，而壓差阻力Fr的計(jì)算涉及到流體的非定常流動(dòng)，計(jì)算相對復(fù)雜，將在下一小節(jié)詳細(xì)闡述.

1.3 壓降和壓差阻力計(jì)算

導(dǎo)向管內(nèi)的流體流動(dòng)可分為兩種情況，一種為圖1所示Q0對應(yīng)的區(qū)間，流體在變截面的圓形導(dǎo)向管內(nèi)流動(dòng)；第二種為圖1所示Q2對應(yīng)的區(qū)間，流體在變截面的導(dǎo)向管和控制棒之間的環(huán)形通道內(nèi)流動(dòng).由于通道的長細(xì)比很大，流動(dòng)可模擬為一維N-S 方程[1]：

其中，t為時(shí)間，P為壓強(qiáng)，C1為導(dǎo)向管的摩擦因數(shù)，l1為導(dǎo)向管的濕周長，S1為控制棒與導(dǎo)向管之間的環(huán)形面積，Kc為局部形阻系數(shù)，由實(shí)驗(yàn)測定，令環(huán)形面積之間的流量為Q=vfS1.將上式在區(qū)間 [Za,Zb]上積分，可得壓降方程：

引入如下變量：

則壓降方程(6)化簡為

式(7)中的I1～I(xiàn)5均有特定的物理意義：I1表示由于局部加速度導(dǎo)致的壓降變化系數(shù)，稱為慣性壓降系數(shù)；I2表示對流加速度產(chǎn)生的壓降系數(shù)，稱為對流壓降系數(shù)；I3表示導(dǎo)向管受到的流體摩擦力產(chǎn)生的壓降；I4表示控制棒受到的流體摩擦力產(chǎn)生的壓降；I5表示由于導(dǎo)向管形狀改變產(chǎn)生的形阻壓降.根據(jù)求得的壓降，可進(jìn)一步求出控制棒在區(qū)間[Za,Zb]受到的壓差阻力[5]：

其中，Srod表示控制棒的橫截面積.

以上分析討論的是導(dǎo)向管內(nèi)存在控制棒時(shí)的壓降和壓差阻力計(jì)算.當(dāng)導(dǎo)向管內(nèi)沒有控制棒時(shí)，式(8)壓降的計(jì)算要減去控制棒摩擦力產(chǎn)生的壓降部分，即減去I4(Za,Zb)，則壓降方程為

2 控制棒與流體耦合運(yùn)動(dòng)狀態(tài)方程

當(dāng)控制棒底部在導(dǎo)向管側(cè)孔上部或者下部時(shí)，流場分布是不同的，壓降方程的建立以及控制棒所受壓差阻力的計(jì)算也是不同的，因此需分別討論.

2.1 控制棒底部在導(dǎo)向管側(cè)孔之上

首先討論控制棒底部在導(dǎo)向管側(cè)孔上部的情況.如1.1 小節(jié)論述，當(dāng)控制棒底部在導(dǎo)向管側(cè)孔之上時(shí)，將導(dǎo)向管分成三個(gè)區(qū)間，分別為[Z,Z2]，[Z1,Z]，[Z0,Z1].在區(qū)間 [Z,Z2]內(nèi)，由于流體在圓環(huán)內(nèi)流動(dòng)，其兩端壓降計(jì)算采用式(8)，可表示為

式(11)中的最后一項(xiàng)表示由于局部損失產(chǎn)生的壓降.為方便表示，記I11=I1(Z,Z2)，I11的第一個(gè)下標(biāo)表示壓降類型，此處為慣性壓降系數(shù)；第二個(gè)下標(biāo)表示導(dǎo)向管的區(qū)間，此處為區(qū)間[Z,Z2].其他壓降類型和其他導(dǎo)向管區(qū)間也類似表示，所以上式變?yōu)?/p>

在[Z1,Z]和 [Z0,Z1]兩段，流體在圓管內(nèi)流動(dòng)，因此可利用式(10)計(jì)算兩端壓強(qiáng)，分別表示為

式(14)中最后一項(xiàng)表示流體從導(dǎo)向管底部流水孔流進(jìn)導(dǎo)向管時(shí)產(chǎn)生的局部突變壓降，其中表示局部壓降損失系數(shù)，SK表示導(dǎo)向管底部流水孔的面積.在導(dǎo)向管側(cè)孔內(nèi)外存在局部突變的壓降[4]可表示為

其中，C3為局部壓降損失系數(shù)[7]；ζC3是開孔的局部水頭損失系數(shù)，一般取1；nC3為導(dǎo)向管側(cè)孔的數(shù)目；SC3為側(cè)孔的面積.根據(jù)式(12)、(13)、(15)可得

根據(jù)式(14)、(15)可得

式(16)中的PZ1-PZ2表示導(dǎo)向管側(cè)孔外與導(dǎo)向管頂部外側(cè)之間的壓降；式(17)中的PZ0-PZ1為導(dǎo)向管底部外側(cè)與導(dǎo)向管側(cè)孔外之間的壓降.這兩個(gè)壓降的計(jì)算通常通過實(shí)驗(yàn)求得.

式(16)與式(17)可以組成一個(gè)非線性方程組，其中的未知量有Q0，Q1，Q2.觀察圖1(a)中的四個(gè)流量及控制棒的運(yùn)動(dòng)速度v，根據(jù)不可壓縮條件，可以得出如下關(guān)系：

式(16)、(17)與式(1)可以組成三個(gè)方程和三個(gè)未知量的非線性微分方程組.此時(shí)可以使用將控制棒的運(yùn)動(dòng)和流體的流動(dòng)解耦的方法求解.首先，假定控制棒運(yùn)動(dòng)狀態(tài)不變，計(jì)算流體流動(dòng)的流量和壓強(qiáng)，從而得到流體阻力；然后，假定流體的流動(dòng)狀態(tài)不變，根據(jù)流體阻力計(jì)算控制棒的速度和位移.計(jì)算控制棒的運(yùn)動(dòng)狀態(tài)時(shí)，文獻(xiàn)[3]采用了向前差分的方法，文獻(xiàn)[4-5]采用了Taylor 展開的方法.然而，這些傳統(tǒng)方法違背了控制棒下落運(yùn)動(dòng)的真實(shí)情況，只能通過不斷縮小迭代時(shí)間步長來提高精度，從而導(dǎo)致計(jì)算效率較低.

本文將控制棒的運(yùn)動(dòng)和流體的流動(dòng)同時(shí)分析，即通過流-固耦合對控制棒下落進(jìn)行求解.此時(shí)引入關(guān)鍵恒等變換=v，將微分方程的個(gè)數(shù)升級(jí)為四個(gè)，非線性微分方程組變?yōu)?/p>

其中

將式(19)寫成狀態(tài)方程的形式：

式(21)即為描述控制棒運(yùn)動(dòng)與導(dǎo)向管內(nèi)流體流動(dòng)的流-固耦合狀態(tài)方程，其中U為狀態(tài)向量，是待求解的未知量.當(dāng)給定初始時(shí)刻的狀態(tài)向量，即可直接對上述非線性微分狀態(tài)方程進(jìn)行時(shí)程積分，得到各個(gè)時(shí)間節(jié)點(diǎn)的狀態(tài)向量，同時(shí)得到控制棒的位移、速度、流體的流量等所有未知量.對于該方程有許多成熟算法，其中保辛算法[8-11]具有精度高和穩(wěn)定性好的優(yōu)點(diǎn)，在解決動(dòng)力學(xué)[12]問題上有獨(dú)特的優(yōu)勢，并在動(dòng)力學(xué)的最優(yōu)控制問題[13-14]、水波問題[15-16]等領(lǐng)域得到廣泛運(yùn)用.本論文選擇辛Euler 中點(diǎn)格式，此部分將在后續(xù)介紹.

2.2 控制棒底部在導(dǎo)向管側(cè)孔之下

同理，如1.1 小節(jié)中所述，當(dāng)控制棒底部在導(dǎo)向管側(cè)孔下方時(shí)，將導(dǎo)向管分成三個(gè)區(qū)間，分別為[Z1,Z2]，[Z,Z1]，[Z0,Z].

類似于上一小節(jié)的方法，可以得到如下非線性微分方程組：

其中

本小節(jié)與上一小節(jié)的不同主要體現(xiàn)在計(jì)算導(dǎo)向管內(nèi)壓降和控制棒所受壓差阻力上.本小節(jié)中控制棒所在區(qū)間為[Z1,Z2]，[Z,Z1]，在計(jì)算區(qū)間[Z,Z1]內(nèi)導(dǎo)向管的壓降時(shí)，需考慮控制棒所受摩擦阻力產(chǎn)生的壓降；同時(shí)，本小節(jié)中控制棒受到的壓差阻力需將這兩個(gè)區(qū)間上控制棒受到的壓差阻力相加.

同理，將式(24)寫成與式(21)完全相同的狀態(tài)方程，其中f1，f2和f3的表達(dá)式見式(25)，而H(U)的表達(dá)式如下：

至此，我們給出了描述控制棒下落與流體流動(dòng)統(tǒng)一的耦合狀態(tài)方程，無論控制棒在導(dǎo)向管側(cè)孔上方還是下方，耦合方程的形式是相同的，即為f(U)=H(U)，不同之處在于H(U)和f(U)的表達(dá)形式.如果考慮控制棒上部分驅(qū)動(dòng)桿結(jié)構(gòu)和流場分析，也可表示為如式(21)所示的統(tǒng)一的狀態(tài)方程求解.

3 時(shí)間步長自適應(yīng)保辛算法

求解式(21)所示狀態(tài)方程，本文選擇辛Euler 中點(diǎn)格式，在此詳細(xì)介紹.假設(shè)狀態(tài)方程(21)中t0時(shí)刻U0已知，要計(jì)算U1.將狀態(tài)方程變形為

則根據(jù)辛Euler 中點(diǎn)格式[9]有

其中

式(28)是一組非線性微分方程，可以迭代求解，首先假設(shè)=U0，然后按n=0,1,2,··· 開始循環(huán)：

直到與之間的相對誤差小于允許值，則停止計(jì)算，此時(shí)U1=.再以U1為初值，計(jì)算U2，U3,···.計(jì)算到臨界邊界條件的值時(shí)，結(jié)束計(jì)算.

在整個(gè)落棒過程中，導(dǎo)向管內(nèi)的流場是非定常的，甚至出現(xiàn)突變，因此狀態(tài)向量的時(shí)程曲線并非光滑變化，將導(dǎo)致數(shù)值計(jì)算的不穩(wěn)定.具體體現(xiàn)為迭代格式(30)在實(shí)際計(jì)算時(shí)不收斂.一般需要將時(shí)間步長取得足夠小才能保證數(shù)值計(jì)算收斂.然而在大部分時(shí)程積分區(qū)間，狀態(tài)向量的時(shí)程曲線是光滑的，時(shí)間步長并不需要取很小就可以保證收斂.綜合以上兩方面考慮，在保證整個(gè)時(shí)程分析收斂的同時(shí)，盡可能地提高計(jì)算效率，本文引入時(shí)間步長自適應(yīng)算法.該方法將采用一個(gè)穩(wěn)定的、較大的時(shí)間步長，記為Δt0（比如0.01 s）.因?yàn)樵诖蟛糠謺r(shí)間段內(nèi)，流場沒有突變，狀態(tài)向量變化是光滑的，因此取較大的時(shí)間步長 Δt0，可以保證迭代的收斂.當(dāng)流場發(fā)生突變時(shí)，仍然取 Δt0會(huì)出現(xiàn)迭代格式不收斂的情況，據(jù)此可以根據(jù)迭代步數(shù)來判斷流場是否發(fā)生突變.在實(shí)際計(jì)算時(shí)，以給定的允許迭代步數(shù)Ni為依據(jù)，如果某次計(jì)算時(shí)，迭代步數(shù)大于Ni，則判斷該時(shí)間步長內(nèi)，流場發(fā)生了突變，這時(shí)將時(shí)間步長減半為0.5Δt0，重新迭代.重新迭代后，如果迭代步數(shù)仍然超過Ni，則將時(shí)間步長再次減半為(0.5)2Δt0，重新迭代；如果迭代步數(shù)沒有超過Ni就收斂，則結(jié)束本次迭代，進(jìn)入下個(gè)時(shí)間步計(jì)算.當(dāng)給定Δt=Δt0時(shí)，在每個(gè)時(shí)間步的具體計(jì)算流程為：

1）令=U0，n=0.

2）判斷n＜Ni是否成立，如果是，計(jì)算式(30)，得到，進(jìn)入步驟3）；如果否，Δt=進(jìn)入步驟1）.

4）U1=，進(jìn)入下一個(gè)時(shí)間步進(jìn)行計(jì)算.

4 算 例

以本文建立的落棒分析模型分析某小型反應(yīng)堆的控制棒下落過程，并將計(jì)算結(jié)果與在核工業(yè)長期使用的某商業(yè)軟件計(jì)算結(jié)果進(jìn)行對比.整個(gè)控制棒下落組件分為上下兩個(gè)部分，這兩部分由星型架連接.其中下部分主要是24 根控制棒插入到24 根導(dǎo)向管中，上部分是驅(qū)動(dòng)桿和驅(qū)動(dòng)機(jī)構(gòu)，驅(qū)動(dòng)機(jī)構(gòu)包括：熱套管、調(diào)節(jié)器、鉤爪、行程套筒等.小堆的主要輸入?yún)?shù)如表1所示.

表1 小堆的主要輸入?yún)?shù)Table 1 Main input parameters for a small reactor

圖2所示為采用本文模型和商業(yè)軟件計(jì)算得到的控制棒在下落過程中的位移、速度和加速度的時(shí)程曲線對比圖，其中實(shí)線為本文方法的計(jì)算結(jié)果，虛線為商業(yè)軟件的計(jì)算結(jié)果.本文模型采用自適應(yīng)時(shí)間步長計(jì)算，時(shí)間步長為Δt0=0.01 s ；商業(yè)軟件的時(shí)間步長為0.001 s.從圖2可以看出，本文模型采用時(shí)間步長 0.01 s計(jì)算得到的位移、速度和加速度時(shí)程曲線與商業(yè)軟件采用時(shí)間步長為0.001 s計(jì)算得到的位移、速度和加速度時(shí)程曲線高度吻合.計(jì)算結(jié)果不僅驗(yàn)證了本文建立的狀態(tài)方程的正確性，且說明了本文方法的高效性.

圖2 小型反應(yīng)堆下本文模型與商業(yè)軟件關(guān)于位移、速度、加速度隨時(shí)間的變化對比Fig.2 The comparison of time-varying displacements,velocities and accelerations,between the proposed model and the commercial software for the case of a small reactor

以小堆為例，本文模型的初始時(shí)間步長 Δt0分別取為0.001 s，0.005 s，0.01 s，0.015 s，商業(yè)軟件的初始時(shí)間步長取為0.001 s.將本文模型得到兩個(gè)關(guān)鍵參數(shù)T5和T5+T6（T5為進(jìn)入緩沖段時(shí)間，T5+T6為控制棒整體落棒時(shí)間）與商業(yè)軟件得到T5和T5+T6進(jìn)行比較，結(jié)果如表2所示.

表2中，括號(hào)內(nèi)的百分?jǐn)?shù)表示本文模型的計(jì)算結(jié)果與商業(yè)軟件的計(jì)算結(jié)果的相對誤差.從表中可以看到，本文使用的4個(gè)時(shí)間步長所得結(jié)果與商業(yè)軟件對比，相對誤差均小于1%.本文使用的最大時(shí)間步長Δt0=0.015 s，是商業(yè)軟件時(shí)間步長的15 倍.而當(dāng)商業(yè)軟件的時(shí)間步長取為0.002 s 時(shí)，計(jì)算結(jié)果發(fā)散.

表2 本文模型在不同時(shí)間步長下與商業(yè)軟件關(guān)于T5 和T5+T6的比較Table 2 The comparison of T5 and T5+T6 between the commercial software and the proposed model for different initial time steps

5 總 結(jié)

針對反應(yīng)堆內(nèi)控制棒下落問題，本文建立了導(dǎo)向管內(nèi)控制棒下落和流體流動(dòng)的理論模型，重點(diǎn)分析了控制棒受到的摩擦阻力、壓降阻力等作用力.在此基礎(chǔ)上，建立了描述控制棒下落與流體流動(dòng)的耦合非線性方程組，再通過引入關(guān)鍵恒等變換v=，得到了流固耦合的非線性狀態(tài)方程，進(jìn)一步通過高精度的辛Euler 中點(diǎn)格式對得到的狀態(tài)方程進(jìn)行了時(shí)程分析.考慮到在不同時(shí)程分析段內(nèi)導(dǎo)向管中流動(dòng)狀態(tài)的突變，本文將時(shí)間步長自適應(yīng)算法引入到時(shí)程分析中.最后以某小型反應(yīng)堆控制棒下落為例進(jìn)行模型和算法的驗(yàn)證.算例中本文模型選取的最大時(shí)間步長 Δt0是商業(yè)軟件時(shí)間步長的15 倍，關(guān)于落棒時(shí)間T5和T5+T6的相對誤差小于1%，從而證明了本文模型和算法的正確性和高效性.本文建立的流固耦合非線性狀態(tài)方程和時(shí)間步長自適應(yīng)算法不僅為控制棒下落分析提供更符合實(shí)際的計(jì)算思路，還為更加復(fù)雜的圓管內(nèi)流動(dòng)模型或者圓環(huán)內(nèi)流動(dòng)模型提供了高精度、高效率的解決方案.

致謝本文作者衷心感謝大連市“青年科技之星”項(xiàng)目（2018RQ06）對本文的資助.