面向高階精度CFD的JFNK算法及其并行計算*

程 彬，李大力，徐傳福，劉 巍，王光學，鄧小剛

1.國防科學技術大學 計算機學院，長沙 410073

2.中國空氣動力研究與發展中心 空氣動力學國家重點實驗室，四川 綿陽 621000

3.國防科學技術大學，長沙 410073

面向高階精度CFD的JFNK算法及其并行計算*

程 彬1+，李大力1，徐傳福1，劉 巍1，王光學2，鄧小剛3

1.國防科學技術大學 計算機學院，長沙 410073

2.中國空氣動力研究與發展中心 空氣動力學國家重點實驗室，四川 綿陽 621000

3.國防科學技術大學，長沙 410073

CHENG Bin,LI Dali,XU Chuanfu,et al.Research on Jacobian-Free Newton-Krylov method for high-order CFD applications and its parallel computing.Journal of Frontiers of Computer Science and Technology, 2017,11(1)：61-69.

目前計算效率低是限制計算流體力學（computational fluid dynamics，CFD）高階精度格式方法的重要因素之一。由于高階精度格式計算模板相對復雜，很難精確計算其Jacobian矩陣，從而影響傳統LU-SGS（lowerupper symmetric Gauss-Seidel）等算法的收斂效率。JFNK（Jacobian-free Newton-Krylov）算法是Krylov子空間方法與非精確牛頓方法的結合，擁有較好的迭代收斂效率，采用無矩陣思想，只計算Jacobian矩陣與矢量的乘積，從而有效避免Jacobian矩陣的計算和存儲。在真實高精度結構網格CFD應用程序中，設計并實現了JFNK時間求解算法。在有粘低速圓柱繞流的算例測試中，和傳統LU-SGS算法相比，JFNK算法擁有更好的計算穩定性，同時可使迭代收斂效率提高2倍以上。以天河2號超級計算機為并行計算平臺，對JFNK算法和傳統的LU-SGS算法的并行強可擴展性進行了測試，二者均表現出良好的并行效率。

高階精度格式；JFNK算法；計算效率；可擴展性

1 引言

計算流體力學（computationalfluid dynamics，CFD）是一門采用數值方法求解流動控制方程，從而揭示流動本質和規律的學科[1]。隨著CFD的發展以及高階精度格式的應用，可模擬的流動問題越來越復雜，時空尺度越來越精細，迫切需要提高模擬效率以加快復雜問題的模擬過程，這就對時間求解算法和并行計算提出了更高的要求。

高精度格式具有計算精度高，流場分辨率高的優點，能夠很好地適用于激波、湍流、邊界層等復雜流體問題的模擬。Harten等人提出了本質無震蕩格式[2]（essentially non-oscillatory scheme，ENO）；Liu等人基于ENO格式提出了加權本質無震蕩格式[3]weighted essentially non-oscillatory scheme，WENO）；鄧小剛等人提出了加權緊致非線性格式[4-5]，包括WCNS（weighted compact nonlinear scheme）、HWCNS（hybrid cell-edge and cell-node weighted compact nonlinear scheme）以及混合型耗散緊致格式[6]hybrid cell-edge and cell-node dissipative compact scheme，HDCS）；Reed和Hill提出了間斷有限元方法[7]discontinuous Galerkin method，DG）。但是相對于二階精度格式，高階精度格式單時間步計算量較大，計算穩定性較差，導致高階精度格式的模擬效率較低，甚至不容易收斂。解決這一問題則需要從時間求解算法和并行計算兩個方面來綜合考慮。

在算法方面，采用高效的迭代求解算法能夠有效減少迭代步數，從而提高CFD模擬的效率，經典的時間求解算法有點松弛SGS（point relaxation symmetric Gauss-Seidel）、線松弛SGS（line relaxation symmetric Gauss-Seidel）、LU-SGS（lower-upper symmetric Gauss-Seidel）、LU-ADI（lower-upper alternating direction implicit）和Jacobi迭代法。經典迭代法計算穩定性高，應用廣泛，然而隨著計算規模持續膨脹，經典迭代法計算效率受到限制，無法滿足工程實際問題需求。Krylov子空間方法作為一種高效的線性系統求解方法在近年來得到了廣泛的研究和應用。Krylov子空間方法是對基于Krylov子空間的一類算法的總稱，包含很多種不同的解法，常用的有CG（conjugate gra-dient）、GMRES（generatized minimal residual）、JFNK（Jacobian-free Newton-Krylov）等。其中，JFNK算法是采用非線性近似計算Jacobian矩陣和矢量乘積的Newton-Krylov子空間方法。研究表明預條件能夠有效提高JFNK算法的收斂速度[8-9]。預條件技術將方程組變換為另一個更容易迭代求解的同解方程組。其實質就是在保持方程組解不變的前提下，縮小矩陣特征值范圍，以改善矩陣的性質，從而達到加速收斂的目的。McHugh等人[10-12]研究了低階空間離散格式CFD應用中各種預條件子在Krylov子空間方法中的效率；Knoll等人[13]研究了多重網格法作為預條件子，將JFNK算法應用于求解不可壓Navier-Stokes方程。上述方法都得到了較好的加速收斂效果。不過關于Krylov子空間方法及其預條件在高階精度格式中的收斂性及穩定性研究還不夠深入全面。

CFD并行計算的基本思想是采用區域分解的數據并行方法，把流場區域劃分為不同的子區域，將流場信息、幾何信息載入相應區域的CPU內存中，每個子區域分配到節點上計算，在子區域內進行串行迭代，一步計算完成后子區域之間交換邊界，然后進行下一步計算，直到收斂。并行算法的并行擴展性對CFD模擬效率同樣有著重要的影響。

本文基于自主開發的真實復雜結構網格高階精度CFD程序HOSTA（high-order simulator for aerodynamics）研究了JFNK求解方法及其并行計算技術。HOSTA采用有限差分法離散控制方程，空間離散采用WCNS[14]高階精度格式，在多塊結構網格上求解Navier-Stokes方程，已應用于航空航天領域等復雜外形復雜流動問題的機理研究和工程應用中。HOSTA中實現了一系列的定常、非定常時間求解方法，包括三步Runge-Kutta、LU-SGS、標量點松弛、矩陣點松弛等。本文在原有時間求解方法上實現了預條件JFNK求解方法，采用圓柱繞流算例深入研究分析了JFNK算法在高精度格式WCNS中的迭代收斂效率和算法穩定性，同時在天河2號超級計算機上實現了JFNK求解算法的并行計算，分析了HOSTA程序的計算效率和并行可擴展性。

2 高階精度格式常用時間求解算法

計算坐標系下，不考慮源項，無量綱Navier-Stokes方程表示為[15]：

式中，為無量綱形式的守恒變量矢量；是l(l=ξ,η,ζ)方向的無粘通量矢量；是l(l=ξ,η,ζ)方向的粘性通量矢量。由式（1）可得：

式中，J為坐標變換Jacobian行列式。

時間求解的收斂速度由左手項決定，一般采用Jacobi、LU-ADI、Gauss-Seidel等線性方程組方法或者Krylov子空間方法求解。其中LU-SGS隱式格式是Navier-Stokes方程定常問題中隱式時間推進方法應用最廣泛的[16]。LU-SGS方法構造簡單，在二階格式中具有穩定性好和收斂效率高的特點；但在高階格式中，由于右端項采用高階精度格式離散，左端項和右端項離散精度不匹配，影響其收斂性和收斂效率[17]。

3 JFNK算法及預條件

JFNK算法是求解大型非線性方程組的算法[10,18]，是非精確Newton法與Krylov子空間法的結合。JFNK算法是在GMRES算法的基礎上，采用無矩陣思想，利用有限差商計算來近似Jacobian矩陣和矢量的乘積計算，因此不需要計算和存儲Jacobian矩陣。在時間步迭代之內，JFNK算法至少包含兩層迭代，外層是非精確Newton迭代，內層是Krylov子空間線性迭代。本文著眼于JFNK算法的預條件等參數對于計算效率的影響，同時對比分析JFNK算法和LU-SGS算法的并行擴展性。

3.1 JFNK算法

JFNK算法的計算步驟如下：

式中，m代表此算法最大內層迭代數；ε是最小范數收斂判據；r0=b-Ax0和wj=Avj采用無矩陣思想，利用有限差商計算即為JFNK算法。

3.2 Jacobian矩陣和矢量乘積的計算

由于Krylov子空間方法只用到Jacobian矩陣與向量的乘積，從而可以不計算和存儲Jacobian矩陣。Krylov空間記為：

令v=Ai-1r(k)，則可采用有限差商近似計算Jacobian矩陣與向量的乘積：

式中，ε對于方程求解的精度和穩定性影響很大，一種計算方法[19]是：

3.3 JFNK算法的預條件

對于復雜數值計算，JFNK算法迭代中需要占用大量計算資源，為了加快迭代收斂速度，考慮到特征值分布越集中收斂越快的特點，在求解過程中，可以采用預條件來改善線性系統。

例如，引入線性左預條件，寫成典型線性方程組形式，然后左乘M-1得：

式中，M是矩陣A的近似矩陣，M-1位于線性層次，迭代過程中不變，為M(Q(k))-1。預條件后的方程殘差為：

預條件后的Krylov子空間為：

Km(M-1A,r(k))≡span{r(k),M-1Ar(k),…,(M-1A)m-1r(k)}采用JFNK算法時，令v=(M-1A)i-1r(k)，則：

本文采用了標量和矩陣的LU-SGS作為JFNK算法的預條件，二者在求解預條件方程時對Jacobian近似矩陣的分裂上存在差異。矩陣預條件采用的是按特征值分裂，而標量預條件采用的是按譜半徑分裂。關于標量和矩陣方法的介紹見文獻[15]。

4 JFNK算法及預條件的并行實現

并行計算方面，HOSTA程序采用了區域分解的思想結合MPI（message passing interface）并行編程實現多進程的并行計算，每次迭代計算過程中僅對邊界信息進行通信交互。

當引入JFNK算法后，需要采用基向量構造Krylov子空間，因此在JFNK算法之前要將原來的多維數據結構映射到一維向量中。采用各塊按序分段映射的方法，一個網格塊對應向量中的一段，保證向量的分段與區域分解的分塊對應，同時又方便塊內的內積運算。

如3.1節的算法流程中所述，JFNK算法中存在向量內積運算，因此需要引入進程歸約通信。首先各個進程對所擁有的向量部分做內積運算；然后通過MPI歸約操作，將各個進程上的局部內積累加到根進程；最后再通過MPI廣播操作將最終的內積結果廣播給所有進程。

另外如圖1和圖2所示，預條件的引入還將涉及到預條件方程右端項的計算（包括通信）以及預條件方程的求解等過程。

5 算例與測試平臺介紹

本文所采用的測試算例是定常低速圓柱繞流問題，分析對比預條件JFNK算法和LU-SGS算法的收斂效率及其影響因素。采用的是單塊網格，網格規模為96×104×5（徑向×周向×展向），網格截面圖如圖3所示。圓柱的無量綱直徑和高度等于1，第一層無量綱網格高度等于0.05。計算狀態為馬赫數Ma=0.1，雷諾數Re=5，T=287 K。使用3種隱式時間推進算法，包括標量LU-SGS算法以及標量預條件和矩陣預條件的JFNK算法。圖4為收斂流場無量綱速度U的分布圖。本文還對算法的強并行可擴展性進行了測試分析，算例網格規模為512×512×5（徑向×周向×展向），利用區域分解技術將網格均勻分割。

Fig.1 Flow chart of HOSTAprogram圖1 HOSTA程序流程圖

Fig.2 Flow chart of HOSTAprogram using JFNK algorithm圖2 使用JFNK算法的HOSTA程序流程圖

Fig.3 Computational grid of flow around a circular cylinder圖3 圓柱繞流算例計算網格

Fig.4 Dimensionless velocity contour distribution of flow around a circular cylinder圖4 圓柱繞流算例無量綱速度分布

本文的測試平臺是天河2號超級計算機系統。天河2系統由16 000個異構體系結構的計算節點組成，總處理器核數為312萬個，雙精度浮點運算的總峰值性能為54 PFlop/s，2013年6月至今一直位居全球超級計算機500強排名榜首。每個天河2號計算節點由2塊Intel Xeon E5-2692處理器和3塊Intel Xeon Phi 31S1P（MIC）協處理器組成。本文測試僅在CPU同構平臺進行，每個計算節點有24個CPU核，共享64GB內存空間，雙精度峰值浮點性能高達422 GFlop/s。

6 收斂效率測試分析

圖5和圖6分別給出了平均殘差相對于時間迭代步和墻鐘時間的平均殘差收斂曲線，結果表明不管是在迭代步數還是墻鐘時間上JFNK算法都要優于LU-SGS算法。JFNK算法殘差達到收斂的時間步數僅為LU-SGS算法的1/4，墻鐘時間也縮短了超過一半，其中矩陣預條件JFNK算法迭代步和墻鐘時間略多于標量預條件JFNK算法。

Fig.5 Average residual curves vs.time step number of LU-SGS and JFNK with nonlinear scalar and matrix type preconditioner圖5 標量、矩陣預條件JFNK和LU-SGS關于時間步的平均殘差曲線

Fig.6 Average residual curves vs.wall time of LU-SGS and JFNK with nonlinear scalar and matrix type preconditioner圖6 標量、矩陣預條件JFNK和LU-SGS關于墻鐘時間的平均殘差曲線

圖7和圖8給出了CFL數對標量預條件和矩陣預條件JFNK算法收斂效率的影響，總體上CFL數越大算法的收斂性越好。當CFL數較小時，迭代收斂速度隨CFL數變化較快，當CFL數較大時，迭代收斂速度隨CFL數變化較慢，最終，收斂速度基本保持不變。對標量預條件和矩陣預條件JFNK算法測試都體現了類似的規律。同時還發現JFNK算法可以取到比LU-SGS算法更大的CFL數，即JFNK算法體現出更優的計算穩定性。綜上分析，JFNK算法擁有相對LU-SGS算法更好的收斂效率和計算穩定性。

Fig.7 Influence of CFL number on convergence rate of JFNK with nonlinear scalar type preconditioner圖7 CFL數對標量預條件JFNK收斂速度的影響

Fig.8 Influence of CFL number on convergence rate of JFNK with nonlinear matrix type preconditioner圖8 CFL數對矩陣預條件JFNK收斂速度的影響

7 并行擴展性分析

圖9給出了單節點內不同進程數的程序強可擴展性結果，其中網格分為24塊，在單節點上采用如橫坐標所示的進程數對其進行測試。從圖中可以看出，在進程數小于4時，并行效率可以保持在80%以上，進程數大于8時，受訪存帶寬限制，程序的并行效率下降明顯。JFNK算法和LU-SGS算法具有類似的問題，因此在進行多節點測試時，每個節點的進程數最多取到4。

Fig.9 Strong program scalability on a single node圖9 單個節點強可擴展性測試

圖10給出了多個節點不同進程的程序強可擴展性結果，其中網格分為64塊，在16個節點上對進程擴展性進行測試，最大進程數為64。從圖中可以看出，隨著進程數增加，程序的并行效率能保持在85%以上，JFNK算法和LU-SGS算法的并行效率非常接近，都具有較好的強可擴展性。

Fig.10 Strong program scalability on multiple nodes圖10 多個節點強可擴展性測試

8 結束語

高階精度格式CFD有限差分模型的計算是非常復雜的，由于難以得到精確的Jacobian矩陣，采用近似Jacobian矩陣的LU-SGS算法收斂效率不高，也無法與高階空間離散格式匹配。JFNK算法避免了計算和存儲Jacobian矩陣，不僅收斂效率更高，同時表現出了更好的穩定性，尤其適合于采用高階空間離散的CFD應用。本文測試算例表明，JFNK算法可以取到比LU-SGS算法更大的CFL數，使用標量預條件、矩陣預條件的JFNK算法都要好于LU-SGS算法，墻鐘時間節省超過一半。

在天河2號超級計算機上進行了并行測試，結果表明JFNK算法和LU-SGS算法都具有良好的強可擴展性，并行計算可以有效地提高HOSTA程序的計算效率，HOSTA程序在高性能應用方面具有優勢。

后續工作除了進一步分析JFNK算法的其他影響因素，包括內層迭代數、不同類型的預處理條件、預處理迭代數、松弛因子等對高階精度CFD應用收斂效率和穩定性的影響，同時將在天河2號超級計算機上對JFNK算法實現CPU+MIC異構協同并行計算和深度性能優化。

[1]Yan Chao.Methodology and applications of computational fluid fynamics[M].Beijing:Beijing University of Aeronautics andAeronautics Press,2006.

[2]Harten A,Engquist B,Osher S,et al.Uniformly high order accurate essentially non-oscillatory schemes[J].Journal of Computational Physics,1987,71(2):231-303.

[3]Liu Xudong,Osher S,Chan T.Weighted essentially nonoscillatory schemes[J].Journal of Computational Physics, 1994,115(1):200-212.

[4]Deng Xiaogang,Mao Meiliang.Weighted compact high-order nonlinear schemes for the Euler equations[C]//Proceedings of the 13th AIAA Computational Fluid Dynamics Conference,Snowmass Village,USA,1997.Reston,USA:AIAA, 1997.

[5]Deng Xiaogang,Mao Meiliang,Jiang Yi,et al.New highorder hybrid cell-edge and cell-node weighted compact nonlinear schemes[C]//Proceedings of the 20th AIAA Computational Fluid Dynamics Conference,Honolulu,USA,Jun 27-30,2011.Reston,USA:AIAA,2011.

[6]Deng Xiaogang,Jiang Yi,Mao Meiliang,et al.Developing hybrid cell-edge and cell-node dissipative compact scheme for complex geometry flows[J].Science China Technological Sciences,2013,56(10):2361-2369.

[7]Zhang Laiping,Liu Wei,He Lixin,et al.A class of hybrid dg/fv methods for conservation laws I:basic formulation and one-dimensional systems[J].Journal of Computational Physics,2012,231(4):1081-1103.

[8]Wu Jianping,Wang Zhenghua,Li Xiaomei.Efficient solving and parallel computing of sparse linear systems[M].Changsha:Hunan Science and Technology Press,2004.

[9]Cai Xiaochuan,Keyes D E.Nonlinearly preconditioned inexact Newton algorithms[J].SIAM Journal on Scientific Computing,2001,24(1):183-200.

[10]Johnson R W,McHugh P R,Knoll D A.High-order scheme implementation using Newton-Krylov solution methods[J]. Numerical Heat Transfer:Part B Fundamentals,1997,31 (3):295-312.

[11]McHugh P R,Knoll D A.Fully implicit finite volume solutions of the incompressible Navier-Stokes and energy equations using inexact Newton's method[J].International Journal for Numerical Methods in Fluids,1994,19(5):439-455. [12]McHugh Knoll P R,D A.Inexact Newton's method solution to the incompressible Navier-Stokes and energy equations using standard and matrix-free implementations[C]// Proceedings of the 11th Computational Fluid Dynamics Conference,Orlando,USA,1994.Reston,USA:AIAA,1994.

[13]Knoll D A,Rider W J.A multigrid preconditioned Newton-Krylov method[J].SIAM Journal on Scientific Computing, 1999,21(2):691-710.

[14]Che Yonggang,Zhang Lilun,Wang Yongxian,et al.OpenMP performance analysis of CFD application on Intel multicore and manycore architectures[J].Journal of Frontiers of Computer Science and Technology,2015,9(10):1153-1162.

[15]Liu Wei,Zhang Lilun,Wang Yongxian,et al.Foundations of computational serodynamics parallel programming[M]. Beijing:National Defense Industry Press,2013.

[16]Jameson A,Yoon S.Lower-upper implicit schemes with multiple grids for the Euler equations[J].AIAA Journal, 1987,25(7):929-935.

[17]Klopfer G H,Hung C M,van der Wijgaart R F,et al.A diagonalized diagonal dominant alternating direction implicit (D3ADI)scheme and subiteration correction[C]//Proceedings of the 29th AIAA Fluid Dynamics Conference,Albuquerque,USA,1998.Reston,USA:AIAA,1998.

[18]Brown P N,Saad Y.Hybrid Krylov methods for nonlinear systems of equations[J].SIAM Journal on Scientific and Statistical Computing,1990,11(3):450-481.

[19]Park H K,Nourgaliev R R,Martineau R C,et al.On physicsbased preconditioning of the Navier-Stokes equations[J]. Journal of Computational Physics,2009,228(24):9131-9146.

附中文參考文獻：

[1]閻超.計算流體力學方法及應用[M].北京:北京航空航天大學出版社,2006.

[8]吳建平,王正華,李曉梅.稀疏線性方程組的高效求解與并行計算[M].長沙:湖南科學技術出版社,2004.

[14]車永剛,張理論,王勇獻,等.Intel多核與集成眾核上CFD程序的OpenMP性能分析[J].計算機科學與探索,2015,9 (10):1153-1162.

[15]劉巍,張理論,王勇獻,等.計算空氣動力學并行編程基礎[M].北京:國防工業出版社,2013.

CHENG Bin was born in 1991.He received the M.S.degree in computer science from National University of Defense Technology.His research interest is applications of high performance computing.

程彬（1991—），男，湖南澧縣人，國防科技大學計算機科學與技術專業碩士，主要研究領域為高性能計算應用。

LI Dali was born in 1987.He is a Ph.D.candidate in computer science at National University of Defense Technology. His research interest is applications of high performance computing.

李大力（1987—），男，湖北荊州人，國防科技大學計算機科學與技術專業博士研究生，主要研究領域為高性能計算應用。

XU Chuanfu was born in 1980.He received the Ph.D.degree in computer science from National University of Defense Technology.Now he is an assistant professor at College of Computer,National University of Defense Technology.His research interests include parallel algorithm and program optimization.

徐傳福（1980—），男，安徽六安人，國防科技大學計算機科學與技術專業博士，現為國防科技大學計算機學院助理研究員，主要研究領域為并行算法，程序優化。

LIU Wei was born in 1980.He received the Ph.D.degree in aeronautical and astronautical science and technology from National University of Defense Technology.Now he is an assistant professor at College of Computer,National University of Defense Technology.His research interest is computational fluid dynamics.

劉巍（1980—），男，吉林遼源人，國防科技大學航空宇航科學與技術專業博士，現為國防科技大學計算機學院助理研究員，主要研究領域為計算流體力學。

WANG Guangxue was born in 1976.He is an associate professor at State Key Laboratory of Aerodynamics,China Aerodynamics Research and Development Center.His research interest is computational fluid dynamics.

王光學（1976—），男，重慶人，中國空氣動力研究與發展中心空氣動力學國家重點實驗室副研究員，主要研究領域為計算流體力學。

DENG Xiaogang was born in 1960.He is a professor and Ph.D.supervisor at National University of Defense Technology.His research interest is computational fluid dynamics.

鄧小剛（1960—），男，四川綿陽人，國防科技大學教授、博士生導師，主要研究領域為計算流體力學。

Research on Jacobian-Free Newton-Krylov Method for High-Order CFD Applications and Its Parallel Computing*

CHENG Bin1+,LI Dali1,XU Chuanfu1,LIU Wei1,WANG Guangxue2,DENG Xiaogang3
1.College of Computer,National University of Defense Technology,Changsha 410073,China
2.State Key Laboratory of Aerodynamics,China Aerodynamics Research and Development Center,Mianyang,Sichuan 621000,China
3.National University of Defense Technology,Changsha 410073,China
+Corresponding author:E-mail:13297499104@163.com

Computational efficiency is a significant factor of high-order accurate scheme for computational fluid dynamics(CFD)finite difference method.It is hard to get Jacobian matrix for lower-upper symmetric Gauss-Seidel(LUSGS)method because of its complicated computing stencil,which affects the computational efficiency.Jacobian-free Newton-Krylov(JFNK)method is a combination of inexact Newton method and Krylov subspace method.With the adoption of matrix free technique,JFNK method only calculates the approximate product of Jacobian matrix and vec-tor using the finite difference quotient of the nonlinear function,thus successfully avoids the calculation and store of approximate Jacobian matrix.Based on the domestic high-order simulator for aerodynamics program,this paper designs and realizes JFNK method.In the viscid and low speed test case of steady flow around a circular cylinder,compared with LU-SGS time stepping method,JFNK method is more stable and efficient.Some numerical experiments of JFNK method and LU-SGS method are performed on Tianhe-2 supercomputer system and good scalability is observed from the test results.

high-order accurate scheme;JFNK method;computational efficiency;scalability

A

：TP301

10.3778/j.issn.1673-9418.1512054

*The National Natural Science Foundation of China under Grant Nos.11272352,61561146395(國家自然科學基金);the Basic Research Program of National University of Defense Technology under Grant No.ZDYYJCYJ20140101(國防科技大學科研計劃重大應用基礎研究項目);the Open Research Program of State Key Laboratory of Aerodynamics of China under Grant No.SKLA20140104 (空氣動力學國家重點實驗室開放課題).

Received 2015-11,Accepted 2016-01.

CNKI網絡優先出版:2016-01-27,http://www.cnki.net/kcms/detail/11.5602.TP.20160127.1635.004.html