車輛-軌道非線性耦合動(dòng)力學(xué)的精細(xì)積分法及其應(yīng)用

張 健，譚述君，2，吳昌華

(1.大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室，大連 116024;2.大連理工大學(xué) 航空航天學(xué)院，大連 116024)

移動(dòng)車輛與軌道相互作用的數(shù)值模擬可高效準(zhǔn)確地獲得二者的動(dòng)態(tài)響應(yīng)，為車輛和軌道設(shè)計(jì)提供理論依據(jù)。車輛與軌道相互作用時(shí)，車輛系統(tǒng)、軌道系統(tǒng)及輪軌關(guān)系存在非線性特性，這時(shí)需要在時(shí)域下通過(guò)逐步積分求解車輛-軌道耦合系統(tǒng)的運(yùn)動(dòng)方程組。

目前計(jì)算車輛-軌道耦合系統(tǒng)的運(yùn)動(dòng)方程的逐步積分方法主要有:Newmark 法［1-2］，Runge-Kutta 法［3］，Admas法［4］，翟婉明的顯式積分法［5-6］和精細(xì)積分法［7］。文獻(xiàn)［1-6］采用逐步積分方法求解，移動(dòng)車輛與軌道相互作用時(shí)，在每個(gè)積分步內(nèi)載荷位置與大小是固定不變的，經(jīng)歷一個(gè)積分步長(zhǎng)后，突然移到下一個(gè)積分步的作用位置。這種計(jì)算方法未考慮移動(dòng)載荷在結(jié)構(gòu)上的連續(xù)移動(dòng)，載荷中的高頻分量很難準(zhǔn)確獲得［8］。然而利用鐘萬(wàn)勰提出的精細(xì)積分方法可求解移動(dòng)載荷在每個(gè)積分步內(nèi)連續(xù)變化情況［8-12］。但文獻(xiàn)［8-9，11-12］忽略質(zhì)量慣性作用對(duì)結(jié)構(gòu)的影響。全玉云［7］采用精細(xì)積分方法分析車輛-軌道耦合系統(tǒng)動(dòng)態(tài)響應(yīng);將車輛和軌道的運(yùn)動(dòng)方程耦合到整體方程中，這導(dǎo)致車輛在軌道運(yùn)行時(shí)系統(tǒng)矩陣有時(shí)變性。當(dāng)系統(tǒng)矩陣發(fā)生改變時(shí)，相應(yīng)的指數(shù)矩陣也要重新計(jì)算。張志超［10］提出利用精細(xì)積分方法計(jì)算移動(dòng)質(zhì)量與橋梁相互作用的迭代算法。而其假設(shè)移動(dòng)質(zhì)量與橋梁始終保持接觸而不發(fā)生分開(kāi)，無(wú)法模擬質(zhì)量與橋梁瞬間脫開(kāi)現(xiàn)象。文獻(xiàn)［7-12］在求非齊次項(xiàng)的Duhamel積分時(shí)需要求系統(tǒng)矩陣的逆陣。譚述君提出非齊次項(xiàng)的Duhamel積分的精細(xì)積分算法［13］，避免計(jì)算系統(tǒng)矩陣求逆。

本文吸收他人的研究工作優(yōu)點(diǎn)，提出利用精細(xì)積分求解移動(dòng)質(zhì)量與結(jié)構(gòu)相互作用的新的方法。采用多項(xiàng)式與指數(shù)矩陣乘法的積分精細(xì)積分算法計(jì)算非齊次項(xiàng)的Duhamel積分，避免系統(tǒng)矩陣求逆;利用Lagrange插值多項(xiàng)式獲得載荷項(xiàng)外插公式，避免迭代求解移動(dòng)質(zhì)量與結(jié)構(gòu)線性/非線性相互作用;分別建立車輛和軌道各個(gè)部件的運(yùn)動(dòng)方程，通過(guò)部件之間相互作用關(guān)系把各個(gè)部件連接起來(lái)，避免當(dāng)車輛在軌道移動(dòng)時(shí)重復(fù)計(jì)算各個(gè)部件的指數(shù)矩陣。通過(guò)有限元形函數(shù)分解移動(dòng)載荷項(xiàng)，獲得移動(dòng)載荷在積分步內(nèi)連續(xù)變化的情況。本文建立移動(dòng)車輛-軌道耦合垂向動(dòng)力學(xué)模型。模型中考慮車輛與軌道非線性輪軌關(guān)系和軌枕懸空后與道床再次接觸情況。分析鋼軌存在焊接接頭不平順時(shí)，軌枕正常支撐狀態(tài)和軌枕懸空對(duì)輪軌力和軌道部件之間相互作用力的影響。

1 精細(xì)積分

系統(tǒng)非線性振動(dòng)方程為:

把式(1)轉(zhuǎn)換到狀態(tài)空間下描述:

其中:

I為單位陣。

利用 Duhamel積分求解式(2)得［16］:

其中T=eHη，η為積分步長(zhǎng)。

根據(jù)鐘萬(wàn)勰［17］提出指數(shù)矩陣的求解方法可以得到T的高精度數(shù)值解。

1.1 非移動(dòng)載荷項(xiàng)處理方法

由于式(1)中的外力項(xiàng)與系統(tǒng)位移有關(guān)，為避免精細(xì)積分方法迭代求解式(1)，利用Lagrange插值多項(xiàng)式近似獲得外力項(xiàng)。

利用的Lagrange插值多項(xiàng)式獲得在第k時(shí)刻到k+1時(shí)刻之間的外力項(xiàng)外插公式:

其中:

fk，fk-1和fk-2是第k，k-1和k-2時(shí)刻的外力項(xiàng)值。

把式(4)代入式(3)可得:

其中:

利用四階Runge-Kutta法解決精細(xì)積分法起步問(wèn)題。

1.2 移動(dòng)載荷項(xiàng)處理方法

當(dāng)載荷的作用位置隨時(shí)間變化時(shí)，利用形函數(shù)獲得移動(dòng)載荷在每個(gè)積分步內(nèi)連續(xù)變化情況［8］。

圖1中，單元長(zhǎng)度為L(zhǎng)。在tk時(shí)刻，載荷作用位置在A處，經(jīng)過(guò)積分步長(zhǎng)η后，作用點(diǎn)移到B處，A和B在同一單元內(nèi)。當(dāng)載荷從A運(yùn)行到C處經(jīng)歷τ時(shí)間，τ∈［0，η］。A和B距節(jié)點(diǎn)I的距離分別為l1和l2。

圖1 移動(dòng)載荷作用位置Fig.1 Location of moving load

移動(dòng)載荷P通過(guò)形函數(shù)分解獲得節(jié)點(diǎn)載荷向量:

其中:E為節(jié)點(diǎn)載荷指示矩陣，E隨載荷作用位置而改變。

n為形函數(shù)組成的向量:

把式(7)代入式(9)得:

其中:bij(i=1～4，j=0～3)的具體表達(dá)式見(jiàn)附錄A。

將式(4)和式(10)代入式(8)可得:

其中:ei=E{d1i，d2i，d3i，d4i}T(i=0 ～ 5)

dij(i=1～4，j=0～5)的具體表達(dá)式見(jiàn)附錄B。

把式(11)代入式(3)可得:

其中:T0，T1和T2具體表達(dá)式見(jiàn)式(6)。

如果共有N個(gè)移動(dòng)載荷同時(shí)作用在結(jié)構(gòu)上，把每個(gè)移動(dòng)載荷項(xiàng)按上述方法計(jì)算非齊次項(xiàng)的Duhamel積分，最后進(jìn)行累加。

當(dāng)積分步長(zhǎng)跨越一個(gè)單元或多個(gè)單元時(shí)，計(jì)算方法與積分步長(zhǎng)不跨越一個(gè)單元時(shí)類似。假設(shè)積分步起始點(diǎn)和終點(diǎn)構(gòu)成虛擬單元，利用形函數(shù)把移動(dòng)載荷分配到虛擬單元的節(jié)點(diǎn)上，然后再次利用形函數(shù)把虛擬節(jié)點(diǎn)上的載荷分配到實(shí)際節(jié)點(diǎn)上。具體推導(dǎo)公式見(jiàn)文獻(xiàn)［12］。

1.3 非齊次項(xiàng)的Duhamel積分的算法

鐘萬(wàn)勰［17］給出式(6)的表達(dá)式:

但是當(dāng)剛度矩陣為奇異矩陣時(shí)，無(wú)法求解系統(tǒng)矩陣H逆陣。然而根據(jù)譚述君［13］提出的多項(xiàng)式與指數(shù)矩陣相乘的積分的精細(xì)積分算法可以精確獲得T0、T1和T2的數(shù)值解，這樣就避免系統(tǒng)矩陣H的求逆。

構(gòu)造非齊次項(xiàng)的Duhamel積分的加法定理:

利用二項(xiàng)式展開(kāi)定理式(16)變?yōu)?

將積分步長(zhǎng)分為2N份:

這時(shí)H(ξ-τ)是非常小的量，可以利用Taylor展開(kāi)求 eH(ξ-τ)，并取展開(kāi)項(xiàng)的前 5 項(xiàng):

通過(guò)執(zhí)行N次加法定理可以得到多項(xiàng)式與指數(shù)矩陣乘積的積分的精確數(shù)值解。

2 車輛－軌道非線性耦合動(dòng)力學(xué)模型

車輛-軌道非線性耦合動(dòng)力學(xué)模型主要包括車輛模型、軌道模型和輪軌關(guān)系。僅考慮車輛和軌道的垂向振動(dòng)，假設(shè)車輛和軌道結(jié)構(gòu)關(guān)于軌道縱向中心線對(duì)稱，因此只建立單側(cè)垂向半節(jié)車輛-軌道耦合動(dòng)力學(xué)模型，如圖2所示。

圖2 垂向車輛-軌道耦合動(dòng)力學(xué)模型Fig.2 Vertical vehicle-track nonlinear coupling dynamic model

鋼軌模型采用離散支撐的兩端簡(jiǎn)支Euler梁模擬，其余部件簡(jiǎn)化成質(zhì)量塊。車輛模型共包括5個(gè)自由度:車體沉浮、轉(zhuǎn)向架沉浮與點(diǎn)頭和輪對(duì)沉浮。鋼軌考慮橫向和彎曲運(yùn)動(dòng)。軌枕和道床只考慮垂向運(yùn)動(dòng)。利用并聯(lián)的彈簧阻尼單元模擬一系和二系懸掛裝置的剛度和阻尼特性、鋼軌與軌枕、軌枕與道床塊、相鄰道床塊之間和道床塊與地基的相互作用。忽略地基的質(zhì)量。

由于沒(méi)有考慮鋼軌的阻尼，因此從方程(1)變成鋼軌的運(yùn)動(dòng)方程為:

其中:Mr和Kr分別為鋼軌質(zhì)量陣和剛度陣;xr和分別為鋼軌的位移和加速度向量;fwr為移動(dòng)輪軌載荷向量;fr為軌枕對(duì)鋼軌作用力向量。

除了鋼軌采用彈性體模擬外，車輛-軌道耦合系統(tǒng)中其他部件采用剛體模擬。

不考慮剛體的阻尼和剛度對(duì)剛體運(yùn)動(dòng)影響，從方程(1)轉(zhuǎn)變成單自由度剛體運(yùn)動(dòng)方程:

例如轉(zhuǎn)向架的運(yùn)動(dòng)方程為:

其中:mb為轉(zhuǎn)向架質(zhì)量，為轉(zhuǎn)向架加速度，F(xiàn)c為車體對(duì)轉(zhuǎn)向架的作用力，F(xiàn)w為輪對(duì)對(duì)轉(zhuǎn)向架的作用力。

輪軌法向力由赫茲非線性彈性接觸理論確定［14］:

其中:Fwr為輪軌力，G為輪軌接觸常數(shù)，xw，xr和xirr分別為車輪位移、鋼軌位移和軌道不平順。當(dāng)輪軌法向力為零時(shí)，表示車輪與鋼軌脫開(kāi)。

懸空軌枕與道床采用非線性彈簧單元模擬［15］，具體表達(dá)式為:

其中:Fsb為軌枕與道床相互作用力，kb為道床剛度，xs、xb和dgap分別為軌枕位移、道床位移和軌枕與道床之間間隙。當(dāng)軌枕與道床之間有相互作用力，這時(shí)反映軌枕與道床再次接觸情況。

采用余弦函數(shù)模擬焊接接頭不平順［14］:

其中:V為車輛運(yùn)行速度，a為不平順波深，L為不平順波長(zhǎng)。

鋼軌運(yùn)動(dòng)方程的載荷向量包括移動(dòng)載荷(輪軌力)和固定載荷(軌枕對(duì)鋼軌的作用載荷)。移動(dòng)載荷和固定載荷都是由兩個(gè)物體位移決定。例如輪軌力是由車輪與鋼軌位移決定的。然而車輪和鋼軌的位移是由輪軌力決定的。所以在每個(gè)時(shí)間積分步內(nèi)，輪軌力和車輪與鋼軌的位移是相互耦合的。利用1.1和1.2節(jié)的方法分別可以解除固定載荷和移動(dòng)載荷與位移的耦合關(guān)系。因?yàn)槔们皫讉€(gè)積分步的載荷信息預(yù)測(cè)當(dāng)前步的載荷信息，再利用當(dāng)前步的載荷求解車輪和鋼軌的位移。解除載荷與位移之間耦合關(guān)系的優(yōu)點(diǎn)是不需迭代求解鋼軌的運(yùn)動(dòng)方程。其它部件的運(yùn)動(dòng)方程處理方法與鋼軌類似。

把方程(22)轉(zhuǎn)換成狀態(tài)空間下的方程(2)形式，其中H矩陣為:

這時(shí)無(wú)法求解H-1，即無(wú)法利用傳統(tǒng)精細(xì)積分求解式(14)中的T0、T1和T2的值。這時(shí)需要采用1.3節(jié)提供的方法獲得T0、T1和T2的值。

因此采用顯示精細(xì)積分方法可以無(wú)需迭代和求解H-1而計(jì)算運(yùn)動(dòng)方程(21)和(22)。

3 焊接接頭和軌枕懸空對(duì)車輛－軌道耦合系統(tǒng)的影響

車輛和軌道具體參數(shù)見(jiàn)文獻(xiàn)［5］。軌道模型長(zhǎng)度24 m。焊接接頭屬于上凸型，位于軌道中央的軌枕正上方，波長(zhǎng)=0.1 m，波深=0.000 3 m。設(shè)焊接接頭下方軌枕發(fā)生懸空，軌枕與道床之間的間隙為0.001 m。車輛運(yùn)行速度為160 km/h。

圖3為軌枕正常支撐和無(wú)軌道不平順時(shí)采用精細(xì)積分與Runge-Kutta法的輪軌力計(jì)算結(jié)果。兩種結(jié)果非常吻合，這表明本文精細(xì)積分方法的正確性。

圖4是軌枕正常支撐和軌枕發(fā)生懸空時(shí)，圖2中右側(cè)輪軌力隨時(shí)間變化情況。在車輪未進(jìn)入焊接接頭不平順和軌枕懸空的影響區(qū)內(nèi)時(shí)，軌枕正常支撐和軌枕懸空時(shí)輪軌力相同(圖4未具體給出)。在車輪未進(jìn)入焊接接頭不平順影響區(qū)，但進(jìn)入軌枕懸空的影響區(qū)時(shí)，兩種情況的輪軌力略有差別，主要是軌枕懸空后，軌道的支撐剛度變低。由于只有一根軌枕發(fā)生懸空，它對(duì)軌道支撐剛度影響范圍很小，所以這時(shí)兩種輪軌力差別不大。當(dāng)車輪進(jìn)入焊接接頭不平順的影響區(qū)內(nèi)時(shí)，軌枕懸空時(shí)輪軌力的第一個(gè)波峰略大于軌枕正常支撐時(shí)的輪軌力。這是因?yàn)榈谝粋€(gè)輪軌力波峰是主要由焊接接頭不平順造成的，懸空軌枕對(duì)輪軌力影響較小，所以兩種情況的輪軌力相差不大。輪軌力的最大值發(fā)生在焊接接頭不平順深度逐漸增加時(shí)不平順的最大斜率附近，這時(shí)不平順深度隨時(shí)間變化非常快。經(jīng)過(guò)斜率最大值后，不平順深度增加速度漸漸降低，輪軌力也隨著降低。當(dāng)不平順深度達(dá)到最高值時(shí)，斜率為零，而輪軌力基本恢復(fù)到無(wú)焊接接頭不平順時(shí)情況。不平順深度漸漸降低和較大輪軌力，造成車輪與軌道發(fā)生脫離現(xiàn)象。兩種情況下車輪和鋼軌脫離時(shí)間幾乎相同，這是因?yàn)閮煞N情況的第一個(gè)輪軌力波峰相差不大。由于車輛重力影響，車輪再次與鋼軌發(fā)生接觸，輪軌力產(chǎn)生第二個(gè)波峰。這個(gè)輪軌力波峰值小于第一個(gè)輪軌力波峰值，軌枕懸空時(shí)的輪軌力明顯大于軌枕正常支撐時(shí)的輪軌力。因?yàn)檫@時(shí)車輪已經(jīng)過(guò)焊接接頭不平順影響區(qū)域，但是還在軌枕懸空區(qū)影響區(qū)域內(nèi)，輪軌力差異由軌道的支持剛度影響。由于輪軌力第二個(gè)波峰持續(xù)時(shí)間較長(zhǎng)，即沖擊力從車輪與鋼軌接觸界面向車輛和軌道傳播時(shí)間長(zhǎng)，會(huì)對(duì)遠(yuǎn)離車輪與鋼軌接觸界面的車輛和軌道的構(gòu)件造成一定損傷。車輪經(jīng)過(guò)軌枕懸空影響區(qū)后，由于車輛和軌道阻尼的作用，輪軌力漸恢復(fù)到未經(jīng)歷軌道不平順情況。圖2中左側(cè)輪軌力與右側(cè)的相似，這里就不再介紹。

圖3 精細(xì)積分與Runge-Kutta法的輪軌力計(jì)算結(jié)果對(duì)比Fig.3 Comparsion of wheel/rail force calculated by precise integration method and Runge-Kutta

圖4 輪軌力時(shí)間歷程Fig.4 Wheel/rail force time history

圖5 輪軌力、鋼軌與軌枕相互作用力和軌枕與道床相互作用力的最大值Fig.5 The maximum wheel/rail force，interaction force of the rail and the sleeper and interaction force of the sleeper and the ballast

在右側(cè)車輪作用下，輪軌力、鋼軌與軌枕相互作用力和軌枕與道床相互作用力的最大值發(fā)生的運(yùn)行時(shí)間，如圖5所示。在軌枕正常支撐時(shí)，當(dāng)輪軌力從車輪與鋼軌接觸界面向下傳時(shí)，由于阻尼吸收能量，它們數(shù)值依次減小，而且有一定時(shí)間差。當(dāng)軌枕與道床發(fā)生再次接觸，它們之間的相互作用力比較小。這是因?yàn)閼铱哲壵淼膬蓚?cè)軌枕承擔(dān)大部分來(lái)自輪軌的作用力。

圖6是右側(cè)車輪通過(guò)焊接接頭不平順時(shí)，焊接接頭下方軌枕及其左右各一根軌枕與道床相互作用力比例因子圖。比例因子是軌枕正常支撐時(shí)或軌枕懸空時(shí)的軌枕與道床相互作用力與軌道無(wú)任何缺陷時(shí)的比值。軌道存在接頭不平順時(shí)，焊接接頭正下方的軌枕與道床相互作用力為軌道正常時(shí)的1.6倍。這就容易造成道床下沉，最終導(dǎo)致道床失去對(duì)軌枕的支撐能力，出現(xiàn)軌枕懸空情況。當(dāng)軌枕懸空出現(xiàn)后，鄰近懸空區(qū)的軌枕與道床作用力明顯增大，這就容易造成軌枕懸空區(qū)進(jìn)一步擴(kuò)大。

圖6 軌枕與道床相互作用力Fig.6 Interaction force of the sleeper and the ballast

4 結(jié)論

利用本文方法求解移動(dòng)質(zhì)量與結(jié)構(gòu)相互作用的特點(diǎn)是不需要求解系統(tǒng)矩陣逆陣，不需要迭代求解考慮移動(dòng)質(zhì)量和結(jié)構(gòu)相互作用的非線性和結(jié)構(gòu)本身的非線性的運(yùn)動(dòng)方程，系統(tǒng)矩陣具有時(shí)不變性，而且可以考慮移動(dòng)質(zhì)量在每個(gè)積分步內(nèi)連續(xù)變化情況。

分析車輛在軌道運(yùn)行時(shí)，輪軌力和軌道部件之間的相互作用力的時(shí)間歷程。計(jì)算結(jié)果表明，在焊接接頭不平順的影響區(qū)內(nèi)，一根懸空軌枕對(duì)輪軌力影響不大。軌枕懸空主要影響第二個(gè)輪軌力峰值。焊接接頭不平順引起較大的輪軌力易造成道床下沉，出現(xiàn)軌枕懸空現(xiàn)象。出現(xiàn)一根軌枕懸空后，易造成鄰近軌枕也發(fā)生懸空現(xiàn)象。

本文方法可以進(jìn)一步應(yīng)用到車輛與無(wú)碴軌道相互作用、車輛與橋梁相互作用和汽車與路面相互作用等移動(dòng)質(zhì)量與結(jié)構(gòu)相互作用的領(lǐng)域。

［1］Lei X，Noda N A.Analyses of dynamic response of vehicle and track coupling system with random irregularity of track vertical profile［J］.Journal of Sound and Vibration，2002，258(1):147-165.

［2］Zhang Y W，Lin J H，Zhao Y，et al.Symplectic random vibration analysis of a vehicle moving on an infinitely long periodic track［J］.Journal of Sound and Vibration，2010，329(21):4440-4454.

［3］ Li Z G，Wu T X.On vehicle/track impact at connection between a floating slab and ballasted track and floating slab track's effectiveness of force isolation［J］.Vehicle System Dynamics，2009，47(5):513-531.

［4］Nielsen J C O，Igeland A.Vertical dynamic interaction between train and track-influence ofwheeland track imperfections［J］.Journal of Sound and Vibration，1995，187(5):825-839.

［5］ Jin X S，Wen Z F，Wang K Y，et al.Three-dimensional train-track model for study of rail corrugation［J］.Journal of Sound and Vibration，2006，293(3-5):830-855.

［6］Zhai W M，Wang K Y，Cai C B.Fundamentals of vehicletrack coupled dynamics［J］.Vehicle System Dynamics，2009，47(11):1349-1376.

［7］全玉云.機(jī)車車輛/軌道系統(tǒng)垂向耦合動(dòng)力學(xué)有限元分析的研究［D］.北京:鐵道部科學(xué)研究院，2000.

［8］林家浩，張守云，呂 峰，等.移動(dòng)簡(jiǎn)諧載荷作用下橋梁響應(yīng)的高效計(jì)算［J］.計(jì)算力學(xué)學(xué)報(bào)，2006，23(4):385-390.

［9］余 華，吳定俊，項(xiàng)海帆.移動(dòng)荷載過(guò)橋的精細(xì)計(jì)算［J］.振動(dòng)與沖擊，2009，28(5):17-21.

［10］張志超，林家浩.移動(dòng)質(zhì)量作用下橋梁響應(yīng)的精細(xì)積分［J］.大連理工大學(xué)學(xué)報(bào)，2006，46(S1):26-32.

［11］ Lu F，Kennedy D，Williams F W，et al.Non-stationary random vibration of fe structures subjected to moving loads［J］.Shock and Vibration，2009，16(3):291-305.

［12］呂 峰，林家浩，張亞輝.移動(dòng)隨機(jī)荷載作用下橋梁振動(dòng)分析［J］.振動(dòng)與沖擊，2008，27(12):73-78.

［13］譚述君，鐘萬(wàn)勰.非齊次動(dòng)力方程duhamel項(xiàng)的精細(xì)積分［J］.力學(xué)學(xué)報(bào)，2007，39(3):374-381.

［14］翟婉明.車輛-軌道耦合動(dòng)力學(xué)(第三版)［M］.北京:科學(xué)出版社，2007.

［15］Bezin Y，Iwnicki S D，Cavalletti M，et al.An investigation of sleeper voids using a flexible track model integrated with railway multi-body dynamics［J］. Proceedings ofthe Institution of Mechanical Engineers Part F-Journal of Rail and Rapid Transit，2009，223(6):597-607.

［16］鐘萬(wàn)勰.應(yīng)用力學(xué)的辛數(shù)學(xué)方法［M］.北京:高等教育出版社，2006.

［17］Zhong W X，Williams F W.A precise time-step integration method［J］.Proceedings of the Institution of Mechanical Engineers PartC-JournalofMechanicalEngineering Science，1994，208(6):427-430.

附錄A

bij系數(shù)表達(dá)式:

附錄B

dij系數(shù)表達(dá)式

同理可得其它值。