于 強，王天舒

(清華大學(xué)航天航空學(xué)院，北京 100084)

0 引 言

一直以來，充液航天器貯箱內(nèi)液體燃料及氧化劑的晃動作用對于航天器的姿軌控制都有著不可忽視的影響，甚至可能導(dǎo)致任務(wù)失敗，對液體晃動問題進行準確分析具有重要的意義。許多學(xué)者都對該問題進行了研究。

等效力學(xué)模型方法是目前廣泛應(yīng)用于航天領(lǐng)域液體小幅晃動分析的研究方法[1-3]，目前也拓展到了液面非線性不太強的大幅晃動問題中[4-5]。然而，當航天器在進行變軌或軟著陸等較大幅度機動的時候，可能會激發(fā)貯箱內(nèi)液體的強非線性晃動，甚至產(chǎn)生液面翻卷破碎的現(xiàn)象。此時等效力學(xué)模型不再適用，合理快速的數(shù)值仿真方法研究有其必要性。此外，航天器在進行大幅度機動的時候往往伴隨著燃料的消耗，貯箱內(nèi)的液體質(zhì)量在不斷減少，所以對液體晃動的數(shù)值仿真方法還要考慮變質(zhì)量的因素。

基于網(wǎng)格的計算流體動力學(xué)方法(CFD)是處理液體晃動問題的一類常用方法，如任意拉格朗日-歐拉(ALE)方法[6]。但網(wǎng)格類方法處理自由液面非線性很強的情況，需要很精細的網(wǎng)格劃分才能比較好地刻畫自由液面，會降低計算效率，在實際工程應(yīng)用時有一定的局限性。光滑粒子流體動力學(xué)(SPH)方法是一種基于拉格朗日描述的無網(wǎng)格粒子類方法，適合自由表面流動問題的仿真[7]。其特點是將流體離散成一系列任意分布的粒子，并利用空間核函數(shù)近似以及粒子求和近似估算場函數(shù)及其空間導(dǎo)數(shù)。粒子之間不存在物理上的連接，在計算大變形自由液面問題時存在優(yōu)勢。

變質(zhì)量晃動問題反映在SPH方法中是流出流量可控的開口邊界問題。邊界條件處理一直都是SPH方法的重要研究方向，大致可分為固壁和開口兩類。目前為止，SPH處理固壁邊界的方法主要可以分為兩大類，第一類是設(shè)置虛粒子的方法，主要有三種實現(xiàn)方式：一是邊界虛粒子對流體粒子直接建立斥力模型[8]，阻止粒子穿透邊界；二是認為相對靜止的邊界虛粒子具有流體粒子的性質(zhì)[9-10]，通過計算合理的邊界虛粒子壓力阻止流體粒子流出；三是通過設(shè)置與流體粒子相對于邊界對稱的鏡像虛粒子，使得流體粒子在固壁邊界處形成近似反彈的效果以達到固壁邊界條件的實現(xiàn)[11]。第二類方法是一致化半解析邊界條件(Unified semi-analytical wall boundary conditions, USAW)處理方法[12]。在該方法中，邊界是按照網(wǎng)格進行劃分，核函數(shù)近似解在邊界被截斷的部分離散到邊界網(wǎng)格節(jié)點上，進而推導(dǎo)出邊界上的控制方程。相比于虛粒子設(shè)置方法，USAW方法具備一定的數(shù)學(xué)理論基礎(chǔ)，可以比較準確地反映固壁邊界附近的流動，并描述任意紐曼(Neumann)或狄利克雷(Dirichlet)邊界條件，但該方法在計算效率方面相對于虛粒子方法有明顯不足。Mayrhofer等[13]將USAW方法應(yīng)用到了三維情況，拓展了其應(yīng)用范圍。駱釗等[14]在USAW方法基礎(chǔ)上提出了無質(zhì)量邊界粒子概念，避免了邊界零粒子層現(xiàn)象。

對應(yīng)固壁邊界條件，SPH方法對于開口邊界條件的處理方法也主要分為虛粒子類方法和半解析類方法兩類。前者主要采用設(shè)置緩沖層并填滿虛粒子的處理方式[15]，借鑒網(wǎng)格類處理開口邊界的思想，為緩沖層虛粒子設(shè)置滿足邊界物理條件的物理量，緩沖區(qū)粒子參與SPH計算。該方法實現(xiàn)簡單，但由于粒子的離散特性，難以準確描述邊界處的流量。后者則利用USAW方法的處理思想控制開口邊界處的流量以實現(xiàn)開口邊界條件。Leroy等[16]結(jié)合ISPH方法和USAW條件，實現(xiàn)了二維及三維情況下定流速開口邊界工況的計算，并得到了很好的效果。但是，半解析的開口邊界處理方法在實際應(yīng)用時，尤其對于三維情況，龐大的計算量限制了其在實際工程領(lǐng)域上的應(yīng)用。

針對充液航天器燃料消耗過程中大機動可能引起的液體晃動問題，本文出于對計算效率以及與航天器動力學(xué)程序聯(lián)合仿真的需求，在非慣性系下SPH(Non-inertial SPH, NI-SPH)方法[17]的基礎(chǔ)上，對虛粒子類開口邊界處理方法進行了改進，使其能夠準確控制出口邊界處的流量，從而實現(xiàn)變質(zhì)量液體晃動的SPH仿真。此外，本文還利用變質(zhì)量質(zhì)點系的動量和動量矩定理計算了液體晃動的作用力及作用力矩，與航天器系統(tǒng)動力學(xué)仿真程序形成了閉環(huán)。最后，通過將仿真結(jié)果與CFD商業(yè)軟件Flow-3d結(jié)果進行對比，驗證了算法的有效性。

1 控制方程

1.1 SPH方法基本理論

標準SPH方法是將液體部分離散成一系列具有質(zhì)量的拉格朗日粒子，不存在碰撞體積，然后利用核函數(shù)近似和粒子求和近似的方法用粒子所攜帶的物理量對宏觀物理量f(xi)(如密度、壓強等)及其空間導(dǎo)數(shù)進行近似，形式如下：

(1)

(2)

其中，下標代表粒子編號，Wij=W(xi-xj,h)為核函數(shù)；x，ρ，m分別表示相應(yīng)粒子的位置向量、密度和質(zhì)量；h表示核函數(shù)的光滑長度，根據(jù)核函數(shù)的不同，其作用范圍，即支持域的半徑為κh(κ取2或3)，h一般可以取粒子初始間距r0的倍數(shù)，即h=kr0；N為粒子i支持域內(nèi)的粒子數(shù)目；此外，有

(3)

本文的核函數(shù)采用的是目前最為廣泛的B-樣條函數(shù)(取κ=2,k=1.3)，如下所示：

(4)

其中,R=r/h，r為粒子間距；對于一維、二維和三維問題分別有

(5)

1.2 N-S方程的SPH離散形式

利用SPH近似方法對拉格朗日描述下Navier-Stokes方程進行離散，可以得到N-S方程的SPH離散形式[18]：

(6)

Γi+fout,i

(7)

其中,v和fout分別代表速度矢量和重力產(chǎn)生的加速度，物理量的下標ij代表粒子i和j的相應(yīng)物理量之差。Γi為簡化的物理黏性耗散項，采用由Lo和Shao[19]提出的公式：

(8)

式中：ν是液體的運動黏性系數(shù)。

1.3 壓力計算方法

SPH方法在計算不可壓縮流體問題時，通過引進人工壓縮率將不可壓縮流近似看作弱可壓流體。對于一般液體的晃動問題，本文采用下面的狀態(tài)方程由密度變化顯式計算壓力[20]：

(9)

式中:ρ0是液體初始密度；γ為常數(shù)，一般γ=7；pb為背景壓力，一般可置為0，式(9)計算的便是流體的相對壓力場；cs為人工聲速，本文取cs=15。γ和cs的選取決定了液體的剛性，限制液體粒子的密度變化率在1%之內(nèi)，能夠保證流體的弱可壓性質(zhì)。

2 邊界處理方法

由于SPH離散粒子的拉格朗日性質(zhì)，SPH方法在表面張力作用不明顯的情況下可以自然滿足自由邊界條件，而對于固壁邊界條件和流量可控的出口邊界條件，本文均基于虛粒子類方法進行了實現(xiàn)。

2.1 固壁邊界條件

基于邊界虛粒子的設(shè)置，本文采用的固壁邊界處理方法詳見文獻[10]?？紤]無滑移邊界條件，給邊界虛粒子設(shè)置了與鏡像虛粒子邊界處理方法意義類似的鏡像速度：

(10)

其中，w代表邊界虛粒子。該速度用于方程(6)和(8)的計算，但不用于虛粒子相對位置更新，邊界虛粒子在非慣性系下的位置始終保持不變。

邊界虛粒子的壓力則需要較好地描述固壁邊界附近的液體壓力梯度，綜合考慮非慣性系下的慣性力和重力，邊界虛粒子的壓力由SPH方法離散得到：

(11)

邊界虛粒子的密度根據(jù)式(11)求出的壓力，利用式(9)的反函數(shù)計算得到。邊界虛粒子參與流體域?qū)嵙Ｗ拥腟PH計算，可以在沒有強沖擊的情況下有效防止流體實粒子的非物理穿透，這種處理方法對于充液航天器貯箱內(nèi)液體晃動仿真非常適用。

2.2 流量可控的出口邊界條件

由于SPH粒子的離散性質(zhì)，SPH方法中的物理量分布是離散的，不是連續(xù)變化的，難以準確刻畫一個截面上物理量的變化率，因此基于虛粒子設(shè)置的緩沖層開口邊界處理方法在處理控制邊界流量的開口邊界問題有著天然的劣勢，無法像網(wǎng)格類方法簡單地在邊界結(jié)點處設(shè)置固定速度。

針對這一問題，本文將貯箱及緩沖層區(qū)域(可認為是貯箱出口管)統(tǒng)一用固壁邊界虛粒子封閉，正常情況下均按照SPH控制方程進行計算。在每一個航天器動力學(xué)時間步Δtd(遠大于SPH仿真時間步)中，根據(jù)所給的當前時刻充液比σ，計算得到該動力學(xué)時間步要流出流體區(qū)域的粒子數(shù)目：

(12)

其中，N和N0分別是上一時刻剩余的流體粒子數(shù)和初始流體粒子數(shù)；σ0是初始的充液比；方括號代表高斯取整函數(shù)。根據(jù)需要流出粒子的數(shù)目ΔN，在緩沖層流出邊界的鄰近網(wǎng)格(鏈表搜索法搜索鄰近粒子對設(shè)置的網(wǎng)格)中搜索最接近邊界的ΔN個粒子設(shè)置為邊緣流出粒子，如圖1所示，記錄此時這些粒子與邊界的垂直距離分別為si0,i=1,2,…,ΔN。這些粒子具有如下性質(zhì)：

圖1 流出邊界附近的粒子分布示意圖Fig.1 Particles distribution near the outlet boundary

1)解除邊界虛粒子對邊緣流出粒子的影響，該類粒子具備垂直于流出邊界向外的速度，其大小為

vi=si0/Δtd

(13)

2)邊緣流出粒子的壓力和密度分別按照固壁邊界條件的方法進行計算。

3)邊緣流出粒子對于其他流體粒子的作用均乘以一個衰減系數(shù)ξi，該系數(shù)由粒子與流出邊界的實時距離si計算得到，其表達式為：

ξi=(esi/si0-1)/(e-1)

(14)

第一條性質(zhì)使得邊緣流出粒子能夠在該動力學(xué)時間步的計算后運動到邊界處，進而被設(shè)置成場外粒子，不再參與計算。而后兩條性質(zhì)則是在保持邊緣流出粒子對于緩沖層粒子作用的同時，保證邊界虛粒子對緩沖層其他流體實粒子的阻擋作用依然合理。該方法可以有效控制流出邊界的流量，并讓出口邊界附近的流場比較合理，減少數(shù)值震蕩。這種離散的邊界處理方法是基于SPH粒子進行人工操作的，數(shù)學(xué)上不能嚴格對應(yīng)流體在邊界處要滿足的方程，但由于緩沖層的存在，這種人工處理方法對于貯箱內(nèi)流體實粒子的擾動很小，可以適用于充液航天器內(nèi)的變質(zhì)量液體晃動問題。此外，方法中衰減系數(shù)的表達式也可以采用其他類型函數(shù)，滿足初始為1，粒子到達邊界時為0即可。

3 變質(zhì)量液體晃動作用力及力矩的求解方法

在實際的應(yīng)用當中，航天器的姿態(tài)信息變化是比較復(fù)雜的，動力學(xué)閉環(huán)仿真中的晃動仿真模塊需要通過航天器的姿態(tài)信息和充液比，反饋該時刻液體晃動對航天器產(chǎn)生的作用力和作用力矩。本文采用非慣性系下的SPH方法[16]，根據(jù)航天器姿態(tài)信息對每一個流體實粒子施加慣性力Finertial,i：

Finertial,i=-mi·

(15)

式中:V0是航天器質(zhì)心的平動速度，dV0/dt為航天器質(zhì)心的絕對加速度，ω和β分別為航天器質(zhì)心本體坐標系的角速度和角加速度的矢量形式，ri為流體實粒子在航天器質(zhì)心本體坐標系下的位置矢量，vi為流體實粒子相對于航天器質(zhì)心本體坐標系的相對速度矢量。

至于液體晃動作用力及作用力矩的求解，本文采用了變質(zhì)量質(zhì)點系的動量定理及動量矩定理進行計算，不計算緩沖層流體實粒子的動量和動量矩，在每一個航天器動力學(xué)時間步Δtd中對累積流出貯箱進入緩沖層的粒子動量和角動量求和，分別得到ΔPout和ΔLout，則該時間步內(nèi)液體晃動對于航天器產(chǎn)生的作用力Fslosh及作用力矩Mslosh分別為：

(16)

(17)

式中：N為該動力學(xué)時間步計算之后的貯箱內(nèi)流體實粒子數(shù)目，Δ后加括號表示括號內(nèi)物理量在該時間步的變化量，finertial,i和fout,i分別為第i個流體實粒子所受慣性力和重力導(dǎo)致的加速度。

SPH仿真液體晃動的計算模塊可以作為一個反饋模塊反映液體晃動產(chǎn)生的影響，與航天器系統(tǒng)動力學(xué)仿真形成一個閉環(huán)，從而實現(xiàn)兩者的聯(lián)合仿真。

4 算例分析和討論

本文所有算例測試均在一臺配置4核8線程4.20 GHz主頻的Intel i7處理器和16 GB內(nèi)存的個人電腦上完成。根據(jù)文獻[21]中的驗證，非慣性系下SPH方法可以在粒子初始間距較大，即粒子數(shù)較少的情況下得到合理的結(jié)果。為了體現(xiàn)算法的效率，本文算例采取了粒子初始間距較大的情況(將粒子初始間距選為貯箱半徑的1/10左右)。此外，算例以CFD商業(yè)軟件Flow-3d的計算結(jié)果作為參照，該軟件采用了VOF方法計算自由表面流，在仿真貯箱內(nèi)液體晃動問題上有穩(wěn)定良好的表現(xiàn)。

4.1 常重環(huán)境球形貯箱短時間仿真

球形貯箱的半徑為0.62 m，貯箱中心位置為(0.62，0.62，0.62)m，液體為水，其密度為1000 kg/m3，動力黏性系數(shù)為9.29×10-4kg/(m·s)。粒子初始間距為0.05 m，初始充液比為15%，設(shè)置了1322個實粒子(包括緩沖區(qū))和2284個邊界虛粒子。初始狀態(tài)下的液體分布及相應(yīng)SPH粒子離散如圖2所示。

圖2 球形貯箱初始液體及SPH粒子分布圖Fig.2 Initial liquid and SPH particles distribution in the sphere tank

仿真時間為10 s，充液比從15%勻速減少到5%?？紤]環(huán)境為常重環(huán)境，重力加速度g=9.8 m/s2，垂直初始水平液面向下(設(shè)為z軸負方向)，然后在水平x和y方向(與z方向構(gòu)成了貯箱的隨體坐標系)分別施加強制的簡諧振動，振動位移形式分別如下所示：

x=0.15sinπt

(18)

y=0.1sin0.25πt

(19)

預(yù)平衡仿真時間2 s(實際仿真前僅施加重力進行的初始狀態(tài)計算)，動力學(xué)時間步長為0.01 s，SPH仿真時間步長Δt由CFL(Courant-Friedrichs-Lewy)穩(wěn)定性條件[22]決定：

(20)

式中:h為光滑長度，cs為人工聲速,vi為粒子i的速度大小。

仿真實際用時52 s，效率上能夠達到嵌入航天器動力學(xué)仿真的要求。表1展示了部分時間點， SPH計算得到的液體粒子分布與Flow-3d計算結(jié)果的對比圖吻合程度較好，SPH方法可以比較好地刻畫自由液面的非線性現(xiàn)象?；蝿幼饔昧傲卦谫A箱隨體坐標系下投影的對比曲線如圖3和圖4所示。

表1 常重情況下球形貯箱內(nèi)變質(zhì)量液體晃動：SPH方法與Flow-3d液面對比Table 1 Variable-mass liquid sloshing in the sphere tank under normal gravity: SPH results of flow patterns in comparison with those from Flow-3d

由圖3～圖4可知，此時SPH方法求解得到的晃動作用力及力矩曲線與Flow-3d仿真結(jié)果吻合良好，計算精度可以得到保證。相比于Flow-3d的計算結(jié)果，SPH求解得到的三軸晃動作用力峰值的相對偏差分別為0.94%，15.74%，9.72%；三軸晃動作用力矩峰值的相對偏差分別為8.64%，9.05%，2.12%。其中相對偏差最大的是水平y(tǒng)方向上的作用力，這主要是由于y方向上的激勵加速度峰值是事實上，包括SPH方法在內(nèi)的粒子類方法會有比網(wǎng)格類方法更強的數(shù)值振蕩，即便采用數(shù)值耗散的方法，也會使計算得到的晃動作用力及力矩存在小幅高頻的非物理振蕩。但是，該振蕩沒有改變作用力變化趨勢，且不會累積沖量，在晃動幅度較大方向上衰減明顯，對整體計算結(jié)果不會產(chǎn)生顯著影響?？傮w而言，對于常重環(huán)境下的變質(zhì)量液體晃動問題，本文提出的方法能夠得到合理的計算結(jié)果。

圖3 常重環(huán)境SPH方法與Flow-3d仿真球形貯箱內(nèi)變質(zhì)量液體晃動的作用力對比曲線Fig.3 Comparison of variable-mass liquid sloshing force histories predicted by the SPH and those simulated by Flow-3d in the sphere tank under normal gravity

圖4 常重環(huán)境SPH方法與Flow-3d仿真球形貯箱內(nèi)變質(zhì)量液體晃動的作用力矩對比曲線Fig.4 Comparison of variable-mass liquid sloshing moment histories predicted by the SPH and those simulated by Flow-3d in the sphere tank under normal gravity

x方向上的1/24，主要的晃動現(xiàn)象會呈現(xiàn)在x這一個方向上(見表1)。從曲線可以看出實際計算的y方向上晃動作用力量級大約是x方向上的1/10，其大小受離散粒子產(chǎn)生的數(shù)值振蕩影響相對較大。

4.2 低重環(huán)境卡西尼貯箱長時間仿真

卡西尼(Cassini)貯箱的中間段為圓柱，兩端為半球，在航天器中有廣泛的應(yīng)用。本節(jié)針對卡西尼貯箱在低重環(huán)境下的復(fù)雜激勵情況，進行了長時間的變質(zhì)量液體晃動仿真。

設(shè)置卡西尼貯箱的半球及圓柱段半徑為0.492 m，圓柱段的高度為0.1 m，貯箱幾何中心在航天器質(zhì)心本體坐標系中的位置坐標為(0.0，-0.875，-0.681)m。液體為航天器中常用的燃燒劑甲基肼(MMH)，密度為874.4 kg/m3，動力黏性系數(shù)為8.50×10-4kg/(m·s)。初始充液比為60%，仿真時間為300 s，勻速減少到35%。粒子初始間距為0.045 m，初始共設(shè)置了3324個流體實粒子(含緩沖層)和1982個虛粒子，初始狀態(tài)下的液體分布及相應(yīng)SPH粒子離散如圖5所示。

圖5 卡西尼貯箱初始液體及粒子分布圖Fig.5 Initial liquid and particles distribution in the Cassini tank

實際航天器在進行變軌或軟著陸時除了用推進的方式，還會使用動量輪調(diào)整姿態(tài)。算例采用了復(fù)雜的運動形式作為激勵輸入，航天器的質(zhì)心平動加速度以及姿態(tài)變化規(guī)律在貯箱隨體坐標系下的投影如圖6和圖7所示(垂直于初始液面向上為z軸正方向)。

圖6 航天器質(zhì)心的平動加速度Fig.6 Translational acceleration of the spacecraft mass center

圖7 航天器角速度Fig.7 Angular velocity of the spacecraft

SPH仿真實際用時2795 s，約為仿真時間的9倍，效率能夠滿足閉環(huán)仿真。計算得到的晃動作用力及力矩在貯箱隨體坐標系下的投影與Flow-3d計算結(jié)果對比曲線如圖8和圖9所示。

圖8 低重環(huán)境SPH方法與Flow-3d仿真卡西尼貯箱內(nèi)變質(zhì)量液體晃動的作用力對比曲線Fig.8 Comparison of variable-mass liquid sloshing force histories predicted by the SPH and those simulated by Flow-3d in the Cassini tank under low gravity

圖9 低重環(huán)境SPH方法與Flow-3d仿真卡西尼貯箱內(nèi)變質(zhì)量液體晃動的作用力矩的對比曲線Fig.9 Comparison of variable-mass liquid sloshing moment histories predicted by the SPH and those simulated by Flow-3d in the Cassini tank under low gravity

從圖8～9可以看出，整體上SPH計算得到的結(jié)果與Flow-3d計算結(jié)果吻合較好；相比于Flow-3d計算結(jié)果，SPH方法求解得到的三軸晃動作用力峰值相對偏差分別為1.64%，10.50%，5.80%；三軸晃動作用力矩峰值的相對偏差分別為6.93%，2.12%，1.35%。

為了更清晰地區(qū)分兩種仿真結(jié)果，對圖8進行了局部放大，選取100～150 s區(qū)間內(nèi)的結(jié)果如圖10所示。

圖10 低重環(huán)境SPH方法與Flow-3d仿真卡西尼貯箱內(nèi)變質(zhì)量液體晃動的作用力局部對比曲線Fig.10 Comparison of variable-mass liquid sloshing local force histories predicted by the SPH and those simulated by Flow-3d in the Cassini tank under low gravity

兩種仿真結(jié)果變化趨勢吻合良好，但通過進一步觀察，可以發(fā)現(xiàn)隨著時間的推移，兩種仿真結(jié)果差異有明顯的增大趨勢，尤其是z方向上的晃動作用力，在仿真最后已經(jīng)可以觀察出有明顯的偏移。經(jīng)過檢查，發(fā)現(xiàn)Flow-3d在計算長時間自由表面流動的時候，由于自由液面處網(wǎng)格之間的流量計算存在誤差，會出現(xiàn)比較明顯的體積計算誤差。針對本算例，F(xiàn)low-3d中液體體積在300 s的仿真過程中大約增加了7%，后期的實際液體體積大于既定值，最后的晃動作用力及力矩也因此相對偏大。當仿真時間進一步加長時，F(xiàn)low-3d的體積累積誤差會進一步增大。此外，根據(jù)實際測試，當自由液面面積更大，變化更加劇烈的時候，液體體積的計算誤差累積速度也會變得更快。因此，考慮到SPH粒子系的質(zhì)量守恒性，SPH方法在長時間仿真時具有更好的穩(wěn)定性。

綜上所述，本文提出的算法能夠解決低重環(huán)境復(fù)雜激勵下長時間變質(zhì)量液體晃動的動力學(xué)問題，對充液航天器貯箱內(nèi)變質(zhì)量液體晃動的仿真分析十分有效。

5 結(jié) 論

本文基于非慣性系下SPH方法的基本理論，針對存在燃料消耗的充液航天器大幅機動情況，提出了一種可以控制流量的虛粒子類出口邊界處理方法。該方法能夠根據(jù)輸入的航天器姿態(tài)信息和充液比，完成單動力學(xué)時間步的SPH更新，并采用變質(zhì)量質(zhì)點系動量定理和動量矩定理計算并輸出晃動作用力及力矩，與航天器系統(tǒng)動力學(xué)仿真形成閉環(huán)。

最后通過與商業(yè)軟件Flow-3d計算結(jié)果的對比，該算法對于解決充液航天器貯箱內(nèi)變質(zhì)量液體大幅晃動問題的適用性和有效性得到了驗證。

SPH方法能夠較好地解決常重及低重環(huán)境下的液體晃動問題，但對于微重環(huán)境表面張力作為主要驅(qū)動力的情況還未有簡潔有效的處理方法，這也是今后需要研究的重要方向。