王少華，閆 賀，張作鵬，刁興中

(清華大學(xué) 核能與新能源技術(shù)研究院 先進(jìn)反應(yīng)堆工程與安全教育部重點(diǎn)實(shí)驗(yàn)室，北京 100084)

控制棒薄壁殼緩沖器是安裝在高溫氣冷堆堆芯石墨結(jié)構(gòu)內(nèi)的重要塑性變形吸能設(shè)備，用于緩沖驅(qū)動(dòng)線斷裂情況下產(chǎn)生的控制棒跌落沖擊。高溫氣冷堆的控制棒運(yùn)行在堆芯側(cè)反射層的石墨孔道中，在極端事故情況下，假設(shè)控制棒驅(qū)動(dòng)線斷裂，由于僅存在氦氣阻力，控制棒的下落近似于自由落體，當(dāng)控制棒下落至孔道底部碰撞到底部的石墨時(shí)會(huì)具有很高的速度。為避免控制棒直接沖擊堆內(nèi)石墨構(gòu)件造成破壞，必須在控制棒運(yùn)行孔道底部石墨上設(shè)置一緩沖吸能裝置，以保護(hù)堆內(nèi)石墨結(jié)構(gòu)。因此，需設(shè)計(jì)且實(shí)現(xiàn)一種可靠的、非能動(dòng)的緩沖器，保證當(dāng)控制棒驅(qū)動(dòng)線斷裂時(shí)底部的石墨不發(fā)生結(jié)構(gòu)性破壞。該緩沖器需滿足以下條件：1) 在滿足使沖擊力峰值盡可能小的同時(shí)，緩沖器必須能吸收最大可能的沖擊動(dòng)能，即當(dāng)約200 kg的控制棒從8 m高處落下時(shí)的動(dòng)能；2) 保證控制棒在與底部石墨發(fā)生碰撞后的回跳盡可能小，避免誤引入正反應(yīng)性；3) 保證控制棒緩沖過(guò)程的平穩(wěn)連續(xù)，緩沖器能自適應(yīng)控制棒的徑向及角度偏差，不發(fā)生側(cè)翻、傾倒；4) 保證在常溫和高溫(550 ℃)條件下緩沖器都能滿足上述要求。

以上條件對(duì)控制棒緩沖器的設(shè)計(jì)提出了苛刻的要求，常見(jiàn)的非能動(dòng)吸能結(jié)構(gòu)包括金屬薄壁殼[1]、復(fù)合材料薄壁殼[2]、泡沫金屬[3]、蜂窩夾層結(jié)構(gòu)[4]等。在反應(yīng)堆工程中，須根據(jù)不同堆型的特點(diǎn)選擇適當(dāng)?shù)目刂瓢艟彌_器，如在FFR實(shí)驗(yàn)快堆中采用了變形套筒式后備減震器[5]、在水力驅(qū)動(dòng)控制棒中采用了水力緩沖器[6]、在壓水堆中曾采用彈簧式緩沖器[7]，綜合考慮高溫氣冷堆的實(shí)際情況，選擇薄壁金屬殼方式是比較適合的。金屬殼材料選擇耐高溫的奧氏體304不銹鋼。金屬殼的形狀采用圓筒型式，這也是控制棒孔道的幾何結(jié)構(gòu)要求的，同時(shí)根據(jù)Jones[8]的總結(jié)，在截面面積相同的情況下，圓形截面的薄壁殼的吸能能力最高，漸進(jìn)屈曲過(guò)程也最為穩(wěn)定。

薄壁殼是一種常見(jiàn)的吸能結(jié)構(gòu)，國(guó)外Alexander[9]、Ari-Gur等[10]及國(guó)內(nèi)王仁等[11]對(duì)此問(wèn)題進(jìn)行了大量研究。HTR-10中也采用薄壁殼結(jié)構(gòu)作為緩沖器，袁碧[12]對(duì)其緩沖效能進(jìn)行了一系列的研究。但HTR-10的緩沖器結(jié)構(gòu)，已不能滿足緩沖200 kg的控制棒從8 m高處跌落的沖擊，因此本文考慮實(shí)際工況，通過(guò)理論分析結(jié)合試驗(yàn)對(duì)緩沖器進(jìn)行設(shè)計(jì)驗(yàn)證。

1 緩沖器設(shè)計(jì)

薄壁殼緩沖器結(jié)構(gòu)如圖1所示，其關(guān)鍵幾何參數(shù)為薄壁殼的壁厚t、長(zhǎng)度L、半徑r。其中，承沖擊環(huán)的作用是將沖擊力均勻加載到薄壁殼，導(dǎo)向外筒用來(lái)限制薄壁殼的徑向偏移，而緩沖器底板可將沖擊力均勻地分散到底部的石墨。

圖1 薄壁殼緩沖器結(jié)構(gòu)

在HTR-10的設(shè)計(jì)[12]中，并未將圖1所示的薄壁殼的直徑與長(zhǎng)度作為重要的考慮因素。事實(shí)上，按照Singace等[13]的考慮偏心效應(yīng)的圓柱殼能量吸收模型，薄壁殼的平均載荷Pavg和緩沖器能吸收的總能量Etotal可由以下公式計(jì)算：

(1)

Pavg/M0=22.27(2r/t)0.5+5.632

(2)

Etotal=ηPavgL

(3)

式中：M0為單位長(zhǎng)度塑性極限彎矩；η為有效長(zhǎng)度系數(shù)，工程上一般取0.75；σ0為304不銹鋼的平均屈服應(yīng)力，可用J-C本構(gòu)模型[14]表示：

(4)

J-C本構(gòu)模型中各參數(shù)的取值列于表1。

表1 J-C本構(gòu)模型的參數(shù)

可看出，當(dāng)吸收的總能量不變時(shí)，增加緩沖器的長(zhǎng)度便能減小所需要的壁厚，同時(shí)減小了薄壁殼的峰值載荷和平均載荷，提高了整體的緩沖效果。但隨高度/半徑比值的增加，薄壁殼屈曲的穩(wěn)定性會(huì)逐漸降低，當(dāng)高徑比達(dá)到一定數(shù)值時(shí)，會(huì)發(fā)生整體的歐拉壓桿式失穩(wěn)，緩沖器吸能效果急劇惡化，因此選擇合適的半徑和高度非常重要。

1.1 薄壁殼半徑的選擇

當(dāng)薄壁殼發(fā)生屈曲時(shí)其水平截面的形狀會(huì)發(fā)生變化，如圖2所示。

圖2 薄壁殼屈曲后水平截面的形狀

由圖2可看出，當(dāng)發(fā)生徑向半波數(shù)為2的非軸對(duì)稱屈曲時(shí)徑向位移最大，若此時(shí)薄壁殼與導(dǎo)向外筒接觸，外筒將會(huì)妨礙薄壁殼的屈曲且額外增加摩擦力，降低緩沖效果。假設(shè)半波數(shù)為2時(shí)薄壁殼的截面變成2條直線，直線長(zhǎng)度為πr，可得出1個(gè)薄壁殼不同導(dǎo)向外筒內(nèi)壁接觸的充分條件：

r≤2R/π

(5)

圖3 半波數(shù)為2的屈曲模態(tài)

式中：r為薄壁殼的半徑；R為導(dǎo)向外筒的內(nèi)徑。通過(guò)式(5)可計(jì)算出薄壁殼的最大半徑。實(shí)際上，薄壁殼不可能屈曲成直線，r≤0.7R的薄壁殼就能滿足屈曲后不與導(dǎo)向外筒的內(nèi)壁發(fā)生較強(qiáng)接觸。利用ABAQUS模擬半波數(shù)為2的屈曲模態(tài)，屈曲后的外型如圖3所示，此時(shí)薄壁殼屈曲后與圓心的最大距離為1.4r，表明半徑r=0.7R的薄壁殼能滿足設(shè)計(jì)需求。

1.2 薄壁殼長(zhǎng)度的選擇

薄壁殼整體是否發(fā)生歐拉失穩(wěn)的臨界長(zhǎng)徑比由多種因素決定，包括載荷的偏心程度、薄壁殼的材料、載荷的加載速度等，無(wú)法獲得一理論解，只能利用試驗(yàn)驗(yàn)證。Hsu等[15]曾使用沖擊錘對(duì)用304不銹鋼制成的薄壁殼進(jìn)行了沖擊試驗(yàn)，根據(jù)其試驗(yàn)結(jié)果，當(dāng)薄壁殼的L/2r達(dá)到6時(shí)開(kāi)始出現(xiàn)整體失穩(wěn)。根據(jù)此試驗(yàn)結(jié)果，可得到緩沖器的總長(zhǎng)度L應(yīng)滿足的條件為L(zhǎng)≤12r，考慮實(shí)際安裝及試驗(yàn)條件本文選擇L=11r。

1.3 薄壁殼壁厚的選擇

圖4 不同壁厚薄壁殼在各溫度下的最大吸收能量

將L和r的具體數(shù)值代入式(1)～(3)，即可得到薄壁殼緩沖器所能吸收的最大能量與壁厚和溫度的關(guān)系，如圖4所示。

反應(yīng)堆中控制棒總質(zhì)量約200 kg，加上鏈條總質(zhì)量約230 kg，最大跌落高度約為8 m，若不計(jì)阻力，控制棒所要吸收的最大能量Etotal=mgh，可得最大沖擊能量為18 kJ，因此選擇壁厚為2 mm的薄壁殼作為吸能體。

2 薄壁殼有限元模型的建立及分析

對(duì)所設(shè)計(jì)的薄壁殼緩沖器建立有限元模型，模擬分析控制棒以13 m/s速度沖擊在薄壁殼上的動(dòng)態(tài)過(guò)程。

真實(shí)的控制棒是由很多節(jié)短棒組成的一個(gè)整體，而每一節(jié)均有復(fù)雜的結(jié)構(gòu)，包括中子吸收體、預(yù)緊彈簧、鉸鏈等，是一類(lèi)似于多節(jié)火車(chē)的結(jié)構(gòu)，如圖5所示?？刂瓢舴譃?0節(jié)，在發(fā)生沖擊前，每?jī)晒?jié)棒之間的可變間隙為2.5 mm，發(fā)生沖擊過(guò)程中，最下面的1節(jié)棒首先沖擊在薄壁殼上，此后9節(jié)棒會(huì)在很短時(shí)間內(nèi)，逐一沖擊在前一節(jié)棒上。理論上，逐節(jié)沖擊的時(shí)間間隔僅為Δt=2.5/13 000=0.2 ms。當(dāng)壓縮波傳至控制棒頂端時(shí)會(huì)反射回拉伸波，拉伸波不能直接穿過(guò)節(jié)與節(jié)之間的間斷面，兩節(jié)控制棒先是在拉伸波的作用下短暫分離，然后通過(guò)鉸鏈傳遞拉伸波。同時(shí)應(yīng)力波在包殼與中子吸收體之間也會(huì)有復(fù)雜的反射與透射現(xiàn)象，因此需對(duì)控制棒建立完整的有限元模型以進(jìn)行計(jì)算。

圖5 控制棒結(jié)構(gòu)示意圖及其中1節(jié)的有限元模型

基于ABAQUS/explicit平臺(tái)建立有限元模型，控制棒整體采用C3D8R單元，網(wǎng)格采用結(jié)構(gòu)化劃分技術(shù)，鉸鏈、石墨孔道和石墨下方的地面設(shè)置為剛體，彈簧簡(jiǎn)化為圓柱形彈性體，薄壁殼采用S4R單元，采用通用接觸，接觸屬性設(shè)置為無(wú)摩擦。

分別計(jì)算從室溫到600 ℃條件下緩沖器的吸能效果，圖6為薄壁殼在常溫下最終的屈曲模態(tài)，圖7為緩沖器受到的最大沖擊力以及緩沖器最終的壓縮比與溫度的變化關(guān)系。由圖7可知，常溫工況下緩沖器受到的沖擊力最大，此時(shí)石墨表面受到的最大壓應(yīng)力為30 MPa，低于石墨的許用應(yīng)力；600 ℃時(shí)緩沖器的壓縮比最大，達(dá)到78%，仍在可接受范圍內(nèi)。

圖6 常溫下薄壁殼受沖擊后的屈曲模態(tài)

圖7 不同溫度下緩沖器的壓縮比和峰值載荷

分析表明，設(shè)計(jì)的薄壁殼緩沖器能有效地吸收控制棒沖擊的能量，保護(hù)石墨的安全。但通過(guò)計(jì)算還發(fā)現(xiàn)，兩節(jié)控制棒之間的沖擊力峰值為1 500 kN，控制棒包殼的最大應(yīng)力為940 MPa，表明此時(shí)控制棒會(huì)發(fā)生一定的破壞。

3 試驗(yàn)設(shè)計(jì)

3.1 試驗(yàn)?zāi)Ｐ偷暮?jiǎn)化

前文分析表明，控制棒在發(fā)生沖擊后，薄壁殼能吸收控制棒的動(dòng)能，避免損壞其底部的石墨結(jié)構(gòu)，但控制棒在沖擊過(guò)程中，存在較大的應(yīng)力，會(huì)出現(xiàn)變形破壞，因此采用實(shí)際結(jié)構(gòu)的控制棒進(jìn)行多次沖擊試驗(yàn)是不可行的。故在試驗(yàn)中采用一等重的單節(jié)實(shí)心棒模型進(jìn)行沖擊試驗(yàn)，考察薄壁殼緩沖器的緩沖效果。

圖8為試驗(yàn)所用的等重鋼棒結(jié)構(gòu)與有限元模型。為研究等重單節(jié)棒的模型是否能有效代表和包絡(luò)多節(jié)棒沖擊的情形，首先對(duì)其建立有限元模型進(jìn)行分析。圖9為單根鋼棒和真實(shí)的多節(jié)控制棒在常溫下以相同的速度(13 m/s)撞擊到緩沖器時(shí)，緩沖器底部的石墨受到的沖擊力。可看出，多節(jié)棒產(chǎn)生的沖擊力較單節(jié)棒多了很多高頻振蕩，反映出多節(jié)棒的沖擊是一不連續(xù)的復(fù)雜過(guò)程；同時(shí)兩者的峰值沖擊力和平均沖擊力大致相等，說(shuō)明試驗(yàn)采用等重的單節(jié)鋼棒模擬真實(shí)的多節(jié)棒是可行的，試驗(yàn)?zāi)茯?yàn)證薄壁殼緩沖器的可靠性。

圖8 試驗(yàn)所用的單節(jié)棒結(jié)構(gòu)和有限元模型

3.2 試驗(yàn)石墨塊的選擇

在真實(shí)的反應(yīng)堆中緩沖器底部是石墨磚結(jié)構(gòu)，試驗(yàn)中同樣采用一石墨材料墊塊模擬反應(yīng)堆中的底部支撐。由于真實(shí)的石墨磚尺寸較大，因此有必要通過(guò)計(jì)算來(lái)分析試驗(yàn)采用何種尺寸的石墨塊最為適合。以一半徑60 mm、高200 mm的石墨圓柱為基礎(chǔ)，通過(guò)不斷加大石墨圓柱塊的尺寸來(lái)觀察控制棒跌落沖擊對(duì)石墨塊的作用，研究石墨表面沖擊力隨石墨尺寸變化的規(guī)律，進(jìn)而推斷采用何種尺寸的石墨塊最適合模擬石墨磚及應(yīng)變信號(hào)采集。

圖9 試驗(yàn)所用的單節(jié)棒和多節(jié)棒的沖擊力

圖10 石墨表面最大壓應(yīng)力隨半徑的變化

由于緩沖器底座半徑為60 mm，以此為最小值，計(jì)算半徑60、80、100、300 mm石墨受到?jīng)_擊時(shí)的受力情況，石墨塊的徑高比保持為3/5。圖10為石墨表面的最大壓應(yīng)力隨半徑的變化，可看出，當(dāng)半徑增加時(shí)最大壓應(yīng)力先是快速下降，當(dāng)半徑超過(guò)100 mm后基本保持不變且均小于石墨的抗壓應(yīng)力極限70 MPa。因此，選擇半徑為60 mm的石墨進(jìn)行試驗(yàn)并以此判斷緩沖器是否滿足要求在工程上是最保守的，同時(shí)選擇更大直徑的石墨塊也是適合的。

試驗(yàn)需通過(guò)測(cè)量石墨側(cè)面的壓應(yīng)變間接測(cè)量沖擊力的大小，圖11為不同半徑的石墨塊側(cè)表面的平均應(yīng)變變化曲線。半徑為60 mm的石墨側(cè)面平均應(yīng)變變化規(guī)律與表面最大壓應(yīng)力變化規(guī)律基本一致，但隨半徑的增加，不但平均應(yīng)變的幅值迅速衰減，而且應(yīng)變變化規(guī)律與沖擊力的偏差也越來(lái)越大，當(dāng)半徑為300 mm時(shí)應(yīng)變幾乎是完全隨機(jī)的，反映到曲線上便是其平均應(yīng)變接近1條直線，即半徑小的石墨塊表面易出現(xiàn)應(yīng)力最大值。綜合這兩方面考慮，試驗(yàn)選擇60 mm半徑的石墨塊是最保守的。

圖11 不同半徑石墨塊側(cè)表面平均應(yīng)變變化

3.3 試驗(yàn)臺(tái)架設(shè)計(jì)

圖12 試驗(yàn)臺(tái)架示意圖

為開(kāi)展8 m控制棒跌落沖擊試驗(yàn)，建立了一高14.4 m的試驗(yàn)臺(tái)架，試驗(yàn)臺(tái)架分為地上和地下兩部分，臺(tái)架結(jié)構(gòu)如圖12所示，地上3層鋼結(jié)構(gòu)臺(tái)架總高7 m，地下2層臺(tái)架總深7.4 m，提升機(jī)、控制棒落棒裝置布置于臺(tái)架的頂層，緩沖器、模擬石墨塊和數(shù)據(jù)采集設(shè)備安裝于臺(tái)架最底層，控制棒的模擬孔道自下而上貫穿臺(tái)架。

試驗(yàn)采用圖8所示的長(zhǎng)2.8 m、外徑110 mm的均質(zhì)鋼棒作為模擬棒，在模擬棒底部加工出1個(gè)倒角與緩沖器頂部的承沖擊環(huán)配合；緩沖器采用全尺寸模型，在緩沖器底部安放1個(gè)圓柱形石墨塊來(lái)驗(yàn)證緩沖效果。石墨底部放有沖擊力傳感器，石墨側(cè)面貼有應(yīng)變片，以此來(lái)衡量沖擊力的大小。

4 試驗(yàn)結(jié)果與分析

沖擊試驗(yàn)后的石墨塊和薄壁殼如圖13所示，試驗(yàn)后石墨塊保持結(jié)構(gòu)完整，證明壁厚2 mm、材料為304不銹鋼的薄壁殼能滿足反應(yīng)堆在冷態(tài)時(shí)的緩沖吸能的作用。同時(shí)屈曲后薄壁殼未與導(dǎo)向筒發(fā)生接觸，薄壁殼發(fā)生了穩(wěn)態(tài)的非軸對(duì)稱屈曲，未發(fā)生整體歐拉失穩(wěn)，與設(shè)計(jì)目標(biāo)一致。

圖13 試驗(yàn)結(jié)束后的薄壁殼和石墨塊

圖14 試驗(yàn)得到的沖擊力與數(shù)值模擬結(jié)果對(duì)比

圖14為石墨所受整體沖擊力試驗(yàn)值與數(shù)值模擬結(jié)果對(duì)比?？煽闯?，因數(shù)值模擬中將地面模擬成剛體忽略了地面的吸能，導(dǎo)致數(shù)值模擬的沖擊時(shí)間偏長(zhǎng)。除此之外，數(shù)值模擬的石墨所受沖擊力的大小與試驗(yàn)結(jié)果符合得很好，石墨所受到的最大沖擊力約為200 kN。

5 結(jié)論

1) 相比HTR-10，提出緩沖器的直徑和長(zhǎng)度要優(yōu)化地選擇合適的值，以防殼體在屈曲過(guò)程中同導(dǎo)向筒接觸，增加額外的沖擊力，同時(shí)防止長(zhǎng)度過(guò)長(zhǎng)導(dǎo)致殼體發(fā)生整體歐拉失穩(wěn)；經(jīng)試驗(yàn)驗(yàn)證，所選擇的半徑和長(zhǎng)度組合是合適的，屈曲后最大半徑小于緩沖器導(dǎo)向外筒，同時(shí)屈曲模態(tài)為穩(wěn)定的非軸對(duì)稱漸進(jìn)屈曲。

2) 模擬分析表明，設(shè)計(jì)的薄壁殼緩沖器能有效地完成吸能緩沖功能，保證石墨基底的完整性，石墨所受到的最大沖擊力為200 kN，試驗(yàn)結(jié)果與模擬結(jié)果符合較好，互相驗(yàn)證了準(zhǔn)確性。

3) 通過(guò)有限元計(jì)算，證明了在試驗(yàn)中采用單節(jié)等重鋼棒代替多節(jié)的控制棒、采用小塊石墨模擬整個(gè)石墨基底的近似是保守的，能包絡(luò)代表實(shí)際工況。

參考文獻(xiàn)：

[1] 杜星文，宋宏偉. 圓柱殼沖擊動(dòng)力學(xué)及耐撞性設(shè)計(jì)[M]. 北京：科學(xué)出版社，2004.

[2] CHING Q S, WAAS A M, HOFFMAN J, et al. The crushing response of braided and CSM glass reinforced composite tubes[J]. Composite Structures, 2001, 52(1): 103-112.

[3] CHEN W, WIERZBICKI T, BREUER O, et al. Torsional crushing of foam-filled thin-walled square columns[J]. International Journal of Mechanical Sciences, 2001, 43(10): 2 297-2 317.

[4] GOLDSMITH W, WANG G T, LI K, et al. Perforation of cellular sandwich plates[J]. International Journal of Impact Engineering, 1997, 19(5): 361-379.

[5] 楊志伸. FFR控制棒驅(qū)動(dòng)機(jī)構(gòu)樣機(jī)的后備減振器[J]. 原子能科學(xué)技術(shù)，1991，25(5)：15-18.

YANG Zhishen. Design of the reserve damper in the FFR control rod drive mechanism prototype[J]. Atomic Energy Science and Technology, 1991, 25(5): 15-18(in Chinese).

[6] 宋威，秦本科，薄涵亮，等. 水力緩沖器三維流場(chǎng)數(shù)值分析[J]. 原子能科學(xué)技術(shù)，2013，47(5)：749-754.

SONG Wei, QIN Benke, BO Hanliang, et al. Numerical analysis of hydraulic absorber flow field[J]. Atomic Energy Science and Technology, 2013, 47(5): 749-754(in Chinese).

[7] 蔣賢國(guó)，胡俊，肖澤軍，等. 反應(yīng)堆控制棒驅(qū)動(dòng)線緩沖試驗(yàn)研究[J]. 原子能科學(xué)技術(shù)，2007，41(增刊)：144-145.

JIANG Xianguo, HU Jun, XIAO Zejun, et al. Research on dashpot structure of rod cluster control assembly of reactor[J]. Atomic Energy Science and Technology, 2007, 41(Suppl.): 144-145(in Chinese).

[8] JONES N. Energy-absorbing effectiveness factor[J]. International Journal of Impact Engineering, 2010, 37(6): 754-765.

[9] ALEXANDER J M. An approximate analysis of the collapse of thin cylindrical shells under axial loading[J]. The Quarterly Journal of Mechanics and Applied Mathematics, 1960, 13(1): 10-15.

[10] ARI-GUR J, WELLER T, SINGER J. Experimental and theoretical studies of columns under axial impact[J]. International Journal of Solids and Structures, 1982, 18(7): 619-641.

[11] 王仁，韓銘寶，黃筑平，等. 受軸向沖擊的圓柱殼塑性動(dòng)力屈曲實(shí)驗(yàn)研究[J]. 力學(xué)學(xué)報(bào)，1983，19(5)：509-514.

WANG Ren, HAN Mingbao, HUANG Zhu-ping, et al. Dynamic plastic buckling of cylindrical shell under axial impact: An experimental study[J]. Acta Mechanica Sinica, 1983, 19(5): 509-514(in Chinese).

[12] 袁碧. HTR-10控制棒減震器設(shè)計(jì)的力學(xué)研究[R]. 北京：清華大學(xué)核能技術(shù)設(shè)計(jì)研究院，2001.

[13] SINGACE A A, ELSOBKY H, REDDY T Y. On the eccentricity factor in the progressive crushing of tubes[J]. International Journal of Solids and Structures, 1995, 32(24): 3 589-3 602.

[14] 魏延鵬. 304及316L不銹鋼激光焊接件力學(xué)性能試驗(yàn)研究[R]. 北京：中國(guó)科學(xué)院力學(xué)研究所，2010.

[15] HSU S S, JONES N. Quasi-static and dynamic axial crushing of thin-walled circular stainless steel, mild steel and aluminum alloy tubes[J]. International Journal of Crashworthiness, 2004, 9(2): 195-217.