基于GPU的快速有限元法求解密度場

李果陽，嚴華，*，張征宇，陳沁梅，祝福順

(1.四川大學 電子信息學院，成都 610065；2.中國空氣動力研究與發(fā)展中心 高速空氣動力研究所，綿陽 621000）

對于風洞試驗，空氣流動包含著最重要的特征信息，對處于空氣流動形成的流場中物理模型的密度場進行定量分析一直都是流場的重要研究方向［1］。傳統(tǒng)的紋影技術由于在定量測量上效果不好而常常作為一種定性分析的方法，干涉技術由于需要嚴苛的要求才能實現(xiàn)，缺乏實用性。因此，背景紋影技術(Background Oriented Schlieren，BOS）就在這種情況下被Meier［2］提出。運用該方法與風洞中的視頻測量技術相結合可以較簡單地實現(xiàn)對處于風洞試驗中的物理模型進行密度場的視頻測量［3-6］，通過測出一束非平行光在流場擾動下的偏折角，進而建立偏折角與擾動流場折射率間的量化關系，得到關于密度場的二階偏微分方程。利用視頻測量法得到關于坐標及坐標偏移量的離散數據，再對于這樣的密度場中二階偏微分方程求解出非平行光經過擾動流場后對應位置的密度投影值。目前有多種常用數值模擬計算方法，經過比較，有限元法因擅長處理各種復雜區(qū)域、精度較高及適于并行計算等優(yōu)點［7-8］，因此本文將采用有限元法求解密度場對應位置的投影值。

從有限元的計算方法可以看出，在用有限元法進行大規(guī)模數值計算時，不僅要考慮多種力學行為的全耦合求解問題，計算過程復雜，而且總剛度方程的求解時間占據了整個計算過程很大的比重［9-10］。由于神經網絡具有強大的映射、實時和并行計算能力［11］，為了解決有限元中的總剛度矩陣求解耗時，簡化有限元法求解的復雜度問題，本文不再采用傳統(tǒng)有限元計算方式，而是引入了神經網絡進行求解。另外，風洞試驗中需要實時分析風洞試驗密度場的分布變化，由此需要快速計算密度場。但是，神經網絡訓練不僅非常耗時且會占據計算機大部分資源，而且有限元計算過程中產生的大規(guī)模稀疏線性方程組數據量也很巨大［12］，計算非常耗時，嚴重影響了密度場的求解效率。對于風洞試驗密度場這樣的時變系統(tǒng)要求實時分析，其核心就在于能夠實時計算，因此，非常有必要探尋求解密度場的快速計算方法。目前主流的高效計算方法主要采用基于GPU的并行計算方式［13-15］，因此本文主要就針對求解過程中耗時的神經網絡、剛度矩陣和載荷向量的求解進行并行加速，以實現(xiàn)密度場的快速求解。

為此，本文提出了一種高效求解密度場的思路，給出了神經網絡與傳統(tǒng)有限元分析融合的有效方法，以及利用GPU快速求解密度場的計算方法和具體步驟，實驗結果證明了所提方法的有效性。

1 有限元法求解密度場

1.1 密度場的二階偏微分方程

利用背景紋影技術的方法測量風洞試驗中的流場［16-17］。坐標系的x軸為逆氣流方向，y軸與x軸垂直，指向試驗段上壁板，z軸與x軸垂直，與光線傳輸方向相反。為了便于分析，不妨假設流動介質的折射率僅與y軸方向的坐標值有關，且隨y軸方向坐標值的增加而增大，沿z軸方向平行射入的光線，在y-z截面內發(fā)生偏折，光線偏折角分析如圖1所示。

圖1 光線的偏折角分析示意圖Fig.1 Schematic diagram of ray deflection angle analysis

折射率n與氣體密度ρ之間的Gladstone-Dale公式為

由光線傳輸方程和式(1），可得y軸方向偏折角εy與氣體密度ρ之間的關系為

同理可得，x軸方向偏折角為

式中：L為風洞試驗段的邊長；KGD為Gladstone-Dale常數，在空氣中標準狀態(tài)下取值為0.226 cm3／g。

對式(4）和式(5）中y、x分別微分再相加即可得到密度場應滿足的二階偏微分方程：

式中：ρm為密度場密度。

將式(6）的等式左端記為Δu，等式右端記為

其求解域為空腔模型，記為Ω。

考慮式(6）的變分問題，對于任意的v∈HΩ，同時乘上式(7）兩邊，使用格林公式，利用邊界條件可得

可證明得到，在一定連續(xù)性要求下，式(6）可等價于求u∈V，使得a(u，v）＝g(u）對任意的v∈V成立。

考慮上述變分問題的有限維逼近，即在有限維子空間Vh?V下，求uh∈Vh，使得a(uh，vh）＝g(uh）對任意的vh∈Vh成立。設在Vh下的一組基函數為，并依次取vh為每個基函數，則可得

分析式(11），其本質是一個線性方程組：

基于上述分析，在初邊值的條件約束下，由Lax-M ilgram定理可知，式(6）即可等價于有限元法下求解一個線性方程組［18-19］：

式中：K為單元力和單元位移關系間的系數矩陣，代表了單元的剛度特性，稱為單元剛度矩陣；F為結點產生單位位移時，在結點上所需要施加的結點力，稱為結點載荷；δ為結點受到擾動流場作用后，結點在其自由度方向偏移的距離，稱為結點位移，且K、F、δ的具體表達式如下：

式中：D為彈性矩陣；t為單元厚度；B為單元應變矩陣；J為雅可比矩陣；α1、α4為剛體位移；α2、α3、α5、α6為單元的常應變；x、y為結點坐標。

1.2 密度場的有限元法求解過程

傳統(tǒng)的有限元分析是利用變分原理將求解域剖分求解的一種數值計算方法，具體步驟如下：

1）連續(xù)體離散化。用有限個離散單元體集合代替原有的連續(xù)體，也稱網格剖分，即進行單元劃分，將全部單元和結點按一定順序編號，每個單元所受荷載均按靜力等效原理作用在結點上，并根據實際情況在位移受約束的結點上設置約束條件。

2）單元分析。建立各個單元的結點位移δe和結點力Fe之間的關系式Fe＝Keδe，Ke表示單元剛度矩陣，即將單元的結點位移作為基本變量，確定一個近似表達式，再應用流體力學理論和虛功原理得到單元中結點力與結點位移的關系式。

3）整體分析。對各個單元組成的整體進行分析，其目的是根據一定法則建立一個線性方程組，得到全部結點載荷F與全部結點位移δ的關系式Kδ＝F，從而求解出結點位移。

本文所實現(xiàn)的快速有限元法采用了神經網絡代替單元分析中求解結點力與結點位移關系的復雜計算，從而使得整體分析計算過程簡化，可大大縮短整體計算時間。

1.2.1 網格剖分

利用有限元法求解泊松方程首要的就是將求解域進行離散，劃分得到有限個互不重疊的單元網格，同時，對這些單元和結點進行編號，這些互不重疊的網格單元通過結點連接相互影響。在劃分單元網格時，對計算結果有著重要影響的是網格單元尺寸大小、網格的疏密程度及網格的類型。三角單元網格能夠適應復雜的求解域，單元大小劃分十分方便，但計算精度相對較低，而四邊形單元劃分不易描述非正交的直邊界和曲邊界，但其計算精度相對較高。

本文采用的是四邊形單元和三角形單元結合的方式，實現(xiàn)了一種四邊形網格單元與三角形網格單元相結合的方式，剖分示意圖如圖2所示。在求解時為便于網格的剖分實現(xiàn)，采用了先對求解域進行等參四邊形單元剖分，再將一個四邊形單元一分為二，形成兩個等參的三角形單元剖分算法，然后對它們進行單元編號和結點編號。

圖2 網格剖分平面示意圖Fig.2 Grid subdivision p lan sketch

1.2.2 神經網絡擬合

利用神經網絡代替?zhèn)鹘y(tǒng)有限元法中的單元分析，關鍵在于要清楚單元分析得到的單元剛度矩陣的特性。

對于本文所采用的三角單元剖分，根據彈性力學二維有限元法［18］，可知有限元法中單元剛度矩陣有以下幾點性質：

1）單元剛度矩陣是對稱方陣。單元剛度矩陣每行元素之和為0，由對稱性質可知，每列元素之和也為0。

2）對于平面單元，單元剛度矩陣不會隨著幾何參數改變。

3）根據單元剛度矩陣的物理意義可以得出其只與單元結點、相對坐標及形函數等相關，與結點絕對坐標無關。

由以上單元剛度矩陣的性質結論可推出，剛度矩陣的計算實質是從單元的尺寸、材料性質到剛度矩陣的映射問題，而神經網絡的優(yōu)勢在于：可以方便地建立輸入和輸出間復雜的映射關系。因此，利用神經網絡可以很方便地從單元結點坐標映射到單元剛度矩陣元素。

本文中，由于密度場二階偏微分方程右端項是關于光線偏折的一階微分式，且實驗數據是使用視頻測量方法在沒有流場擾動時，拍攝一幅包含高密度標記點的背景板圖像作為參考圖像，然后加入擾動流場，再拍攝背景板的時序圖像，根據視頻測量的標定原理測得背景板標記點在坐標系下的(x，y）坐標值及其偏移量Δx和Δy。因此在求解時需要對數據進行預處理，即根據視頻測量法原理和光線偏折理論得到坐標偏移量和偏折角的關系，再對其做歸一化處理，接著采用神經網絡對偏移角擬合得到二階偏微分方程右端項一階微分式，即通過對右端項的一階微分式在點(x，y）處根據泰勒級數展開式可以得到

式中：h為坐標點更新的步長。

神經網絡擬合模型如圖3所示。將經過預處理的標記點坐標數據(其中80%為訓練樣本，20%為驗證樣本）作為網絡輸入，通過一個含有400個神經元的隱含層，激活函數采用tanh型，優(yōu)化函數采用目前普遍使用的Adam進行有監(jiān)督的學習。將網格單元結點通過該訓練好的神經網絡，最終實現(xiàn)網格單元的結點通過神經網絡映射到單元剛度矩陣元素。

圖3 神經網絡擬合模型Fig.3 Neural network fitting model

1.2.3 整體分析

總剛度矩陣K由單元剛度矩陣Ke組裝而成，計算如下：

式中：Ce為單元選擇矩陣。根據單元的結點編號和結構中結點編號的對應關系，將單元剛度矩陣進行式(19）迭加組裝成總體剛度矩陣，得到結點載荷與結點位移間的關系式，如式(13）所示。

然后進行狄利克雷邊界條件約束，采用直接法求解得出結點位移。

2 基于GPU的有限元法并行求解

本文中將有限元法求解密度場的方法采用并行的思想，對求解過程極其耗時且適合并行計算的神經網絡、總剛度矩陣和總載荷向量的求解進行并行加速。融合了神經網絡的單元分析過程，以時間代價簡化了求解復雜度，為進一步加快求解速度采取了數據并行、多線程同時計算的方式。整體分析中，總體剛度矩陣和總載荷向量的求解進行分塊多線程并行計算。

2.1 CUDA架構介紹

GPU具有高性能的浮點計算能力，CPU具有低延遲、邏輯計算能力強的特點，CUDA編程正是利用GPU和CPU這種特點來協(xié)同處理任務。運行流程則是在CPU端準備數據，再傳到GPU端并行計算，最后將計算結果回傳到CPU端。CPU和GPU的這種協(xié)同工作方式大大提高了計算性能［20-21］。圖4為CUDA線程組織結構模型。

為了明白CUDA編程如何進行線程組織，必須了解其結構模型。從圖4中可以看到，CUDA中線程層次分為3層，分別是線程網格(grid）、線程塊(block）、線程(thread）。1個線程索引(thread Idx）可以是1～3維的，進而多個thread組成的block也可以是1～3維的。如同線程組成線程塊，線程塊也能組成1～3維的線程網格。執(zhí)行CUDA應用程序時，整個應用程序的關鍵是kernel函數，它是以grid的形式組織，以block為單位，各block間并行執(zhí)行，不同block間的數據共享只能通過全局顯存，對于同一block內的不同thread間則可以通過共享內存進行通信，即在kernel函數中存在2個層次的并行，即block之間并行計算及thread之間并行計算。

圖4 CUDA線程組織結構模型Fig.4 CUDA thread organization structure model

另外，為了達到更高的效率，在CUDA 編程應格外關注內存的使用。在CUDA編程時要盡量使用寄存器，盡量將數據聲明為局部變量。而當存在著數據的重復利用時，可以把數據存放在共享內存里。而對于全局內存，需要注意用一種合理的方式來進行數據的合并訪問，以盡量減少設備對內存子系統(tǒng)再次發(fā)出訪問操作的次數。

2.2 神經網絡并行加速方法

本文中將有限元法求解密度場的方法采用并行的思想對求解過程極其耗時的部分進行加速處理，其模型如圖5所示。數據擬合的并行化主要是針對神經網絡的學習訓練過程。首先，將在密度場測量得到的離散數據分為訓練樣本、驗證樣本和測試樣本，并將其暫放在內存中，然后對神經網絡進行初始化參數，并且將訓練樣本分成多個批量(batch），加載到GPU顯存中。在神經網絡正向傳播過程中，將原本一整個數據分為多片，調用CUDA的kernel函數為每個batch數據分配處理線程。由于在CUDA中一個W rap線程束有32個線程，在進行線程分配時應當盡量選擇32的整數倍，但是CUDA中一個線程塊里線程數量又不多于1 024，因而可以盡可能提高線程的利用率。從主機端調用kernel函數在GPU設備端對分片好的訓練樣本在大量線程中被同時并行處理，這就極大地加快了正向傳播過程。同理，它的反向傳播過程也是各線程中進行一系列的矩陣運算和向量運算。當每個線程所傳遞的誤差信號都到達網絡的隱含層，就能得到權值偏導和偏置偏導，而每個batch上得到的是誤差信號均值，并且只對其本身參數偏導進行更新，所有線程間互不干擾，所有的參數間也不相關，這就完成了一次并行迭代過程。

圖5 神經網絡并行模型Fig.5 Neural network parallel model

2.3 總剛度矩陣組裝的并行化

2.3.1 總剛度矩陣的并行性分析

對于大規(guī)模的有限元分析中，串行的有限元過程中均不是直接在內存中一次建立總的平衡方程組，因為對于大規(guī)模的有限元問題，其結點總數也很大，導致總剛度矩陣的階數太大，從而計算機計算過程中內存不夠。因此，串行有限元過程中一般有2種處理方法：一是將總剛度矩陣分塊存儲在外存中，二是邊單元分析邊消去，避免總剛度矩陣的生成。但是這2種方法都需要內外存之間頻繁地進行數據交換，因此對于結點數較多的大規(guī)模問題計算一次是非常耗時的。

總剛度矩陣是將所有單元剛度矩陣按照其結點編號對號入座組裝而成，所有單元的單元剛度矩陣不能向總剛度矩陣一次性迭加而成，否則就可能在計算過程中發(fā)生計算錯誤和存取沖突。進一步分析可知，總剛度矩陣組裝最大的并行只可能發(fā)生在每次同時迭加互相沒有共同結點的多個單元的剛度陣所有元素，即實際采用的方案可以是每次同時迭加一個單元剛度矩陣多個位置處的元素，或者每次同時迭加多個單元剛度矩陣在同一位置處的元素。另外，總剛度矩陣的組裝與其存儲格式密切相關，大型有限元問題產生的總剛度矩陣一般是稀疏對稱的。

2.3.2 總剛度矩陣的并行組裝

通過上文分析，本文中的問題則是一個稀疏對稱的總剛度矩陣，對于這樣的總剛度矩陣組裝過程很耗時但卻非常適合并行計算?？倓偠染仃囉擅總€網格單元的單元剛度矩陣組裝而來，每個網格單元的單元剛度矩陣的值只與其形狀大小、對應點次序及在整體坐標系的方位有關，因此本次求解劃分的矩形網格單元內的三角單元具有相同的單元剛度矩陣。圖6表示了劃分的任意一個網格位置結點標號。

圖6 剖分的任一網格位置結點編號Fig.6 Node number of any grid position of subdivision

本文單元剛度矩陣組裝成總剛度矩陣的具體方法是：將單元剛度矩陣Ke的每個子塊Kij送到總剛度矩陣的對應位置上去，然后進行迭加即可得到總剛度矩陣K的子塊，對每個單元都進行如此操作計算即可得到總剛度矩陣。

在執(zhí)行并行組裝總剛度矩陣時，關鍵是如何找出Ke中的子塊在K中的對應位置，然后為每個Ke的組裝分配線程索引，這樣在并行化過程中給每個單元分配的線程求解得到的單元剛度矩陣才能正確地組裝成總剛度矩陣，解決這個問題需要清楚單元中結點編碼(局部碼）和整體結構中的結點編碼(總碼）之間的對應關系。圖7中每個單元的3個結點按逆時針方向編為i、j、m，單元剛度矩陣中的子塊是按結點的局部碼排列的，而總剛度矩陣中的子塊是按總碼排列的，因此，在單元剛度矩陣中把結點的局部碼換成總碼，并把其中的子塊按照總碼次序重新排列，以單元U1為例，局部碼i1、j1、m1對應于總碼m、m＋1、m＋x＋1，即單元剛度矩陣和總剛度矩陣位置對應關系如圖7表示。其中上三角網格單元和下三角網格單元的單元剛度矩陣相同。

圖7 單元剛度矩陣與總剛度矩陣位置對應關系Fig.7 Location of element stiffness matrix corresponding to total stiffness matrix

另外，有限元法產生的總剛度矩陣的階數一般很高，在解算時矩陣的階數很高和存儲容量有限是矛盾的，因此有必要知道總剛度矩陣的特性來解決節(jié)約存儲容量的問題。經分析，總剛度矩陣的特性如下：

1）對稱性。只存儲總剛度矩陣的上三角部分，可以節(jié)約一半的存儲容量。

2）稀疏性。總剛度矩陣的元素絕大部分都是零元素，非零元素只占一小部分。

3）帶形分布?？倓偠染仃嚨姆橇阍胤植荚谝詫蔷€為中心的帶形區(qū)域內，在半個帶形區(qū)域內，包括對角線元素，每行具有的元素個數叫做半帶寬，即可以只存儲該部分元素值來節(jié)約存儲空間。

若存儲時不采取有效的存儲方式，不僅需要大量的存儲空間、時間，還會增大計算量。因此，本文在實現(xiàn)過程中根據總剛度矩陣的特性采用了較為常用的CSR矩陣壓縮存儲格式，用3個一維數組分別存儲稀疏矩陣中上三角的半帶寬非零元素值、對應的列號及行偏移，使得GPU并行更易實現(xiàn)。

在并行實現(xiàn)時，GPU運行kernel核函數，每個CUDA線程執(zhí)行一個網格單元剛度矩陣寫入總剛度矩陣的任務，主要分為如下步驟：

步驟1根據參數計算出網格單元的剛度矩陣。

步驟2根據網格單元的結點編號，判斷網格單元是上三角形單元還是下三角形單元。

步驟3根據判斷結果將值寫入總剛度矩陣對應的位置。

為能夠并行計算剛度矩陣和載荷向量并且不產生線程和存儲沖突，本文根據不同單元網格結點編號進行固定空間的索引存取，即每次同時迭加所有單元的單元剛度矩陣同一位置的元素?？倓偠染仃嚳傒d荷向量和結點編號密切相關，因此在GPU并行求解總剛度矩陣和總載荷向量時，最為關鍵的是如何將各結點寫入到GPU存儲空間。對于本次求解的網格單元，每個結點的載荷在F中最多會迭加6次。為避免同時迭加或者使用CUDA原子操作，CUDA編程時為每個結點分配一個大小為6的空間，即最開始的總載荷向量F大小為N×6。實現(xiàn)中要保證下三角結點的編號順序始終為(m，m ＋x＋1，m ＋x），上三角順序(m＋x＋1，m，m＋1），把所有下三角的結點編號按順序設為索引1、2、3，同理把所有上三角的結點依次設為4、5、6，如圖8所示。相同索引在同一結點處只會出現(xiàn)一次或者不出現(xiàn)，因此在求得結點之后先將載荷向量寫入對應結點對應的索引處，如對于m點，對應的位置為6×m＋Index(索引），這樣即可線程和存儲不沖突，從而正確求解。

圖8 結點存儲位置Fig.8 Node storage location

當所有值全部寫入上述數組之后，則每個結點處真正的載荷值為6個值相加，此時每個CUDA線程只需要將6個數據相加便可得到最后的總剛度矩陣和總載荷向量。

3 實驗驗證

3.1 實驗環(huán)境

運用上述串并行實現(xiàn)原理方法，進行了密度場有限元法求解。實驗平臺操作系統(tǒng)采用Ubuntu 16.04，CPU為Intel i7-6800k，GPU為NVIDIA GTX 1080Ti。實驗數據由某風洞研究所提供。對比串行和并行實現(xiàn)下不同的網格精度、各主要模塊的計算耗時來驗證并行快速有限元法的有效性。

3.2 串行速度

表1列出了不同精度大小的網格單元剖分下主要模塊串行求解時間。圖9展示了有限元法同網格精度下主要模塊占總運行時間的比例。根據圖9結果進行分析，在不同網格規(guī)模下，有限元法求解密度場過程中網絡訓練和網絡預測時間開銷相差較小，但是隨著網格規(guī)模加大，求解總剛度矩陣和總載荷向量的時間開銷卻成倍增加，從求解總時間來看這整個過程極其耗時，完全無法達到風洞試驗分析場密度實時要求。

表1 不同模塊下CPU串行求解時間Table 1 CPU serial solution time under different modules

圖9 不同網格精度下主要模塊占總運行時間的比例Fig.9 Proportion ofmain modules in total running time with different grid precision

3.3 并行速度

為了解決不能實時計算密度場的問題，本文著力于尋求一種基于GPU的快速求解方法，使用現(xiàn)在流行的GPU高效計算硬件平臺，極大提高了求解效率。分別進行了如表2所示幾組實驗。圖10展示了有限元法不同網格精度下串行求解和并行求解密度場的加速比；圖11展示了有限元法不同網格精度下串行和并行求解密度場主要模塊的加速比。本文所指的加速比皆是指串行求解和并行求解運行時間的比值。

圖10 不同網格精度下求解密度場的加速比Fig.10 Acceleration ratio of solved density field with different grid precision

圖11 不同網格精度下求解密度場主要模塊的加速比Fig.11 Acceleration ratio of main modules of solved density field with different grid precision

表2 不同模塊下GPU并行求解時間Table 2 GPU parallel solving time under different modules

從不同的網格劃分疏密程度的實驗結果對比可以看到，剖分的網格較為稀疏時求解密度場總體加速比較小，大概在25.6倍左右，隨著剖分網格越密，總體加速比都在逐漸增加，這是由于網格精度增加，有限元法產生的稀疏矩陣規(guī)模越大，利用GPU并行求解總剛度矩陣和總載荷向量的優(yōu)勢越明顯，而網絡訓練的加速比改變不大是因為導入的實驗數據量并沒改變。另外，可以看到求解過程中的求解總剛度矩陣和總載荷向量占總運行時間的比例逐漸上升，導致這一現(xiàn)象的主要原因是隨著網格數量增加，總剛度矩陣的組裝和總載荷向量的累加操作需要大量的共享存儲器，每個流處理器同時運行的block快速減少，從而GPU的利用率下降，運行時間增加。

3.4 精度比較

本文中采用L2范數和L∞范數來計算誤差，計算如式(20）和式(21）所示，分別表示本文快速有限元法數值計算結果密度xi與精確結果x＾的均方誤差，以及本文快速有限元法數值計算結果密度xi與精確結果絕對值的最大誤差。

式中：Δh為剖分的網格單元面積。

從圖12可以看出，有限元法串行和并行求解密度場的精度隨著網格劃分變密其精度增加，但精度增加緩慢，且串行計算精度高于并行計算，這是因為并行計算加入了一些累加計算，所以截斷誤差影響更顯著。另外，該方法精度在絕大多數應用場景中都處于可接受水平，如在網格精度為496×450時，其均方誤差和最大絕對值誤差均在10－5量級下，在絕大部分場景下，這一誤差可以完全忽略。

圖12 有限元串并行精度比較Fig.12 Serial and parallel precision comparison of finite element method

4 結論

1）給出了基于GPU實現(xiàn)的快速有限元求解風洞試驗密度場的方法。通過對比實驗分析對串行方法中耗時較多的神經網絡擬合、總剛度矩陣組裝及總載荷向量求解實現(xiàn)了GPU并行加速。實驗表明，隨著網格數量的增加，其加速比可達數百倍以上，且其精度達到了工業(yè)要求，提高了有限元法求解密度場的時效性。

2）實現(xiàn)的快速有限元法相對于傳統(tǒng)的有限元簡化了對剛度矩陣的求解，大大減少了計算時間。

3）隨著網格數量增加，加速比相應增加，其精度也逐步提高，但求解更耗時。