龍格庫塔間斷有限元方法在計算爆轟問題中的應(yīng)用

《龍格庫塔間斷有限元方法在計算爆轟問題中的應(yīng)用》由會員分享，可在線閱讀，更多相關(guān)《龍格庫塔間斷有限元方法在計算爆轟問題中的應(yīng)用（13頁珍藏版）》請在裝配圖網(wǎng)上搜索。

1、 ?龍格庫塔間斷有限元方法在計算爆轟問題中的應(yīng)用 張磊1,2[收稿日期] xxxx-xx-xxx；[修改日期] xxxx-xx-xx [基金項目] 國家973計劃（2005CB321703）和國家自然科學基金（10531080,10729101）資助項目 ， 袁禮1 1. 中國科學院數(shù)學與系統(tǒng)科學研究院計算數(shù)學所，科學與工程計算國家重點實驗室，北京100190； 2. 中國礦業(yè)大學(北京)理學院，北京100083 [摘要] 構(gòu)造了求解帶源項守恒律方程組的龍格庫塔間斷有限元（RKDG）方法，并分別結(jié)合源項的Strang分裂法和無分裂法數(shù)值求解模型守恒律方程和反應(yīng)歐

2、拉方程。為了和有限體積型WENO方法進行比較，設(shè)計了計算源項的WENO重構(gòu)格式。對一維帶源項守恒律的計算表明，對于非剛性問題，RKDG方法比有限體積型FVWENO方法的誤差更小，而對于剛性問題，RKDG方法對于間斷面位置的捕捉更為精確。對于一二維爆轟波問題的計算結(jié)果表明，RKDG方法對于爆轟波結(jié)構(gòu)的分辨和爆轟波位置的捕捉能力更強。 [關(guān)鍵詞] 龍格庫塔間斷有限元方法；爆轟波；反應(yīng)Euler方程；剛性源項 [中圖分類號] O241.82 [文獻標識碼] A 引言 在非平衡氣體動力學中，爆轟問題通常被描述成無粘性的化學反應(yīng)流動問題，其物理模型是一個非齊次雙曲守

3、恒律方程組，通常被稱作反應(yīng)歐拉方程組（reactive Euler equations），其中非齊次項（源項）通常被解釋為由于化學反應(yīng)引起的混合組分的質(zhì)量變化率[1]。 最簡單的反應(yīng)歐拉方程組假定混合氣體僅僅由兩種組分構(gòu)成：已燃氣體（burnt gas）和未燃氣體（unburnt gas）[2]。當混合氣體達到點火溫度時，未燃氣體通過一個不可逆的化學反應(yīng)轉(zhuǎn)化為已燃氣體，因此，混合氣體狀態(tài)可以用一個標量，即未燃氣體的質(zhì)量分數(shù)Y來表示，進一步假設(shè)混合氣體各組分具有相同的比熱比和比氣體常數(shù)R，則二維反應(yīng)歐拉方程組可以寫成如下形式 ,

4、 （1） 其中，各個變量和函數(shù)的具體形式為 式中是混合氣體的密度，分別是沿方向的速度分量，分別是壓力，總能量，溫度。稱為化學反應(yīng)率，通常用Arrhenius模型 （2） 其中是化學反應(yīng)率因子，是活化能。此外混合氣體的狀態(tài)方程和溫度分別為 （3） 其中為單位質(zhì)量未燃氣體發(fā)生化學反應(yīng)所釋放的熱量。 在實際爆轟波問題中，化學反應(yīng)的時間尺度遠小于流體流動的時間尺度，所以非齊次方程組具有很強的剛性，給數(shù)值求解帶來很大困難。無論對源項用不用算子分裂方法，采用顯格式還是隱格式，在處理間斷解問題時，都可能得到

5、非物理解。這是因為當源項不能用足夠的空間和時間分辨率進行求解時，激波捕捉方法會得到錯誤的爆轟波傳播速度[3,4,5]。文[5]發(fā)現(xiàn)，當空間分辨率不夠時，一個同時含有已燃氣體和未燃氣體的網(wǎng)格中會發(fā)生虛假反應(yīng)，數(shù)值爆轟波會以每個時間步一個空間步長的速度傳播。文[6]也指出計算的數(shù)值誤差極可能促成溫度敏感化學反應(yīng)的提早發(fā)生。因此，要獲得正確的爆轟波位置，常用的解決方案是采用合適的點火溫度模型來抑制虛假反應(yīng)[3]，或者采用網(wǎng)格自適應(yīng)方法[7, 8, 9]來保證反應(yīng)區(qū)域內(nèi)有足夠多的網(wǎng)格點，或者采用高精度高分辨率方法來保證對反應(yīng)區(qū)有足夠高的分辨率。由于后兩種途徑有更好的普適性而受到更多的關(guān)注。 由Coc

6、kburn和Shu等人發(fā)展的龍格庫塔間斷有限元方法（Runge-Kutta Discontinuous Galerkin Method，RKDG）是一類具有高精度和高分辨率的數(shù)值方法[10]，在解決含有間斷現(xiàn)象的問題中發(fā)揮著越來越重要的作用，它被廣泛地發(fā)展和應(yīng)用于水動力學，氣動力學，波傳播，半導(dǎo)體中的電荷傳輸?shù)葐栴}。該方法既保持了一般有限元方法和有限體積方法的優(yōu)點，又克服了各自的不足。該方法可采用局部高階插值的方法構(gòu)造基函數(shù)，具有靈活處理間斷和邊界條件以及可顯式求解的能力，克服了一般有限元方法不適于間斷問題的缺點，以及一般有限體積方法必須通過擴大模板進行重構(gòu)來提高精度的不足。RKDG方法所具有

7、的單元上連續(xù)分布的高精度高分辨率逼近特性有望能更好地模擬爆轟波問題。 本文應(yīng)用RKDG方法，結(jié)合化學源項的無分裂方法及算子分裂方法，對爆轟波問題進行數(shù)值模擬，并和有限體積型WENO（FVWENO）方法的計算結(jié)果進行比較。由于FVWENO方法需要利用流動變量單元平均值來計算源項的單元積分，本文采用Simpson積分公式，其積分點處的流動變量采用WENO重構(gòu)以提高計算精度。對于帶源項的一維守恒律的數(shù)值測試結(jié)果表明，對于非剛性問題，RKDG方法計算結(jié)果的誤差更小；而對于剛性問題，RKDG方法對于間斷波位置的捕捉能力更強。對典型的一二維爆轟波算例，包括二維不穩(wěn)定爆轟波和爆轟波繞過90°拐角的衍射問題

8、的數(shù)值模擬結(jié)果顯示RKDG法仍有一定優(yōu)勢。 1 數(shù)值方法 1.1 空間離散 本文的空間離散采用RKDG方法[10]。設(shè)是區(qū)域的一個有限剖分。單元表示多邊形單元的一條邊界，表示單元邊界的外法向。是單元上的局部有限元空間，取作次多項式集合。，在間斷有限元空間 中尋找近似解，其中。 首先在單元上用試探函數(shù)乘以方程（1）的兩端，并用近似解代替方程（1）的精確解，用代替試探函數(shù)，將（1）寫成變分形式，并由分部積分和Green公式，可得 （4） 對于（4）式的后三項，采用數(shù)值積分計算 （5） （6） （7） 流通量采

9、用與之相容的數(shù)值流通量代替，數(shù)值流通量定義為為某種形式的Riemann解算器。分別代表邊積分和單元積分的Gauss積分點個數(shù)，和分別代表邊積分和單元積分的Gauss積分點位置[10]。這樣得到半離散的數(shù)值格式 （8） 為了計算方便，在單元中取正交基函數(shù)（如勒讓德多項式），則質(zhì)量矩陣成為對角矩陣，故有限元解表示為 （9） 在式（8）中，取試探函數(shù)為所有基函數(shù)，得到半離散的常微分方程組 （10） 其中，為單元的質(zhì)量矩陣，和分別為 （11） 和

10、 （12） 高階RKDG方法會產(chǎn)生數(shù)值振蕩，需要在每一個時間步后對數(shù)值解應(yīng)用斜率限制器，本文采用文獻[10]的TVB限制器。當基函數(shù)的階數(shù)增加時，流通量的選擇對于數(shù)值結(jié)果的影響較小[10]，所以本文采用局部的Lax-Friedrichs流通量。 1.2 無分裂方法和分裂方法 對于空間離散后得到的常微分方程組（10），如果源項的剛性不太強，可以直接采用常微分方程組的求解方法進行時間推進。本文的無分裂方法是采用三階強穩(wěn)定保持的龍格庫塔（SSPRK）方法[11]直接求解方程組（10），該方法在空間離散取階基函數(shù)時，可以保證光滑解的三階時間和空間精度。但當源項的剛性比較強

11、時，通常要對流動項和化學源項分裂求解，即對于方程組（10）分別求解流動項的半離散方程 （13） 和化學源項的半離散方程 （14） 本文采用的分裂法是具有二階時間精度的Strang分裂方法[12] （15） 其中和分別表示以和為時間步長求解（14）和（13）的算子。 在分裂方法中，守恒律的半離散方程組（13）仍采用顯式的三階的SSPRK方法求解，而源項的方程組（14）采用常微分方程求解器VODE求解[13]。 1.3 有限體積型WENO方法 本文將比

12、較RKDG方法和高精度的FVWENO方法[14]。在§1.1節(jié)介紹的空間離散中，如果基函數(shù)只取為1，則可得到守恒律方程組（1）的有限體積方法。在二維矩形網(wǎng)格上，半離散有限體積法的方程為 （16） 上式中對流項的離散采用維數(shù)分裂的基于特征變量的五階WENO格式[14]。在應(yīng)用無分裂方法計算時，由于FVWENO方法的待求變量是流動變量的網(wǎng)格單元平均值，如果直接用它計算（16）式中源項的積分，則只有二階空間精度。為盡量提高精度，本文對一維問題采用四階的Simpson公式來計算積分 （17） 其中積分點的值和采用基于特征變量的五階WENO重構(gòu)得到。

13、對于的重構(gòu)會產(chǎn)生負權(quán)，故需要對重構(gòu)結(jié)果進行處理[15]。 對于二維問題，可以采用（17）式的張量積形式，這時積分具有四階精度，但涉及重構(gòu)積分點處值，計算公式相對復(fù)雜。為計算方便，本文采用一種具有三階精度的二維Simpson積分 （18） 其中由方向的單元平均值重構(gòu)得到，由方向的單元平均值重構(gòu)得到，采用兩個方向重構(gòu)值的平均值。 2 數(shù)值算例 2.1 精度測試 第一個例子是一個帶有源項的Burgers方程 （19） 這個方程的解是連續(xù)的，其解析解為 （20） 式中控制方程的剛性，即為初值。為了測試算法

14、的精度，取，此時問題的剛性不強。計算區(qū)間取為[-15,25]，RKDG方法中TVB限制器的常數(shù)取為5，計算到0，CFL數(shù)取為0.18。表1是數(shù)值解的誤差及其數(shù)值精度階。 表1 帶有源項的Burgers方程 Table 1 Burgers equation with source term N Strang split method Unsplit method DG2 R WENO5 R DG2 R WENO5 R WENOM R 32 7.66′10-5 － 4.79′10-3 － 7.55′10-5 － 4.78′10-3 －

15、 5.25′10-3 － 64 5.29′10-6 3.86 9.62′10-4 2.31 4.77′10-6 3.98 9.64′10-4 2.31 5.04′10-4 3.38 128 4.58′10-7 3.53 2.42′10-4 1.99 1.91′10-7 4.64 2.43′10-4 1.99 3.19′10-5 3.98 256 8.37′10-8 2.45 5.76′10-5 2.07 1.36′10-8 3.81 5.76′10-5 2.07 2.41′10-6 3.72 512 1.84′

16、10-8 2.18 1.43′10-5 2.01 1.42′10-9 3.26 1.43′10-5 2.01 2.80′10-7 3.11 1024 4.33′10-9 2.09 3.56′10-6 2.00 1.69′10-10 3.07 3.56′10-6 2.00 3.53′10-8 2.99 2048 1.05′10-9 2.04 8.90′10-7 2.00 2.09′10-11 3.02 8.90′10-7 2.00 4.46′10-9 2.99 注：DG2表示P2基函數(shù)的RKDG方法，WENO5表示五階F

17、VWENO方法， WENOM表示五階FVWENO方法+源項四階精度積分公式（17） 從表1中可以看出，用Strang分裂方法，由于在時間上只有二階精度，在時間步長與空間步長同階的情況下，即使采用高精度的基函數(shù)DG方法也只能達到二階精度。而無分裂的元DG方法和修正FVWENOM方法可以達到和SSPRK方法同階的三階精度，但沒有采用高精度源項積分的FVWENO只達到二階精度，無論是分裂方法還是無分裂方法，RKDG方法在網(wǎng)格規(guī)模相同的情況下，其數(shù)值解的誤差明顯要比FVWENO方法的小得多。 2.2 間斷計算能力測試 第二個例子是一個帶有源項的對流方程

18、 （21） 取初值為 （22） 是間斷的初始位置，這個例子的解析解是一個以速度1向右傳播的間斷解。 文[5]指出，當網(wǎng)格數(shù)不夠多時，計算的間斷波速度變慢甚至是不動的。為了衡量方法的間斷捕捉能力，定義數(shù)值平均速度 其中為一維空間步長，為計算中止時間，和分別是時刻初值和數(shù)值解的求和。 計算中取，CFL數(shù)取0.1，計算區(qū)間為[0,2]，初始間斷位置為。 表2 帶有剛性源項的對流方程 Table 2 Linear advection equation with stiff source term N Strang spli

19、t method Unsplit method DG2 WENO5 DG2 WENO5 WENOM 256 5.3′10-4 4.7′10-4 0.7106 4.8′10-4 4.8′10-4 512 0.7249 4.0′10-2 0.9029 7.8′10-4 0.4032 1024 0.9332 0.8145 0.9689 0.8094 0.8327 2048 0.9877 0.9398 0.9956 0.9387 0.9421 從表2中可以看出，在網(wǎng)格規(guī)模相同的情況下，RKDG方法的數(shù)值平均速度比

20、FVWENO的更接近于1，說明前者對于間斷波的捕捉更為準確。而對于同一數(shù)值格式，當網(wǎng)格數(shù)較少時，計算的數(shù)值平均速度是不正確的。只有當網(wǎng)格數(shù)足夠多時，各種格式得到的數(shù)值速度才趨近于1。無分裂的DG2和WENOM比對應(yīng)的分裂方法的結(jié)果稍微準確些，這與無分裂方法比相應(yīng)的分裂方法的精度高一階有關(guān)。 2.3 一維ZND爆轟波 第三個算例是一維穩(wěn)定爆轟波問題，初始時刻，假定爆轟波的von Neumann點位于。首先考慮一個CJ-ZND爆轟波情形。未燃氣體狀態(tài)為 狀態(tài)方程及Arrhenius模型中各參數(shù)分別為。這些參數(shù)下理論解半反應(yīng)長度為1[3]。計算區(qū)間取為[-40,40]，網(wǎng)格數(shù)N=400，

21、該網(wǎng)格數(shù)可以分辨反應(yīng)區(qū)，計算的CFL取為0.18。該問題的理論爆轟波速度為5.4419。計算到時間，此時爆轟波傳播到27.2附近。圖1是計算的密度分布的局部放大圖，可以看出RKDG方法的解比FVWENO的更接近于精確解，且間斷面過渡網(wǎng)格點要少。但是由于RKDG方法使用了TVB型的限制器（限制器常數(shù)取為1.0），算出的von Neumann點比FVWENO的略低，而分裂方法對于解的影響幾乎看不出。 圖1 CJ-ZND爆轟波結(jié)果局部放大圖 Fig.1 Detail view of CJ-ZND detonation wave 圖2 強ZND爆轟波結(jié)果局部放大圖 Fig.2

22、 Detail view of strong ZND detonation wave 其次，我們考慮一個強爆轟情形，其未燃氣體狀態(tài)和前面的CJ-ZND爆轟波波前的相同，狀態(tài)方程及Arrhenius模型中各參數(shù)分別為 。這時半反應(yīng)長度為0.01[3]，計算區(qū)間為[-10,110]，網(wǎng)格數(shù)N=800。該網(wǎng)格數(shù)對這個問題是“under-resolved”，即無法分辨反應(yīng)過程中組分的空間變化，只能捕捉爆轟波的位置。該情形下理論強爆轟波的速度為4.9093，計算到時間，此時爆轟波傳播到98.186附近。圖2是局部放大圖，圖中實線是用4000個網(wǎng)格求得的理論解?？梢钥闯鯮KDG方法的間斷面過渡網(wǎng)格

23、點數(shù)比FVWENO的更少，而分裂方法對于兩種數(shù)值方法的影響不大。 2.4 二維不穩(wěn)定爆轟波 第四個算例考慮一個初始擾動了的CJ-ZND爆轟波經(jīng)過一段時間發(fā)展之后，它的化學反應(yīng)區(qū)里具有胞格這樣一種規(guī)則的結(jié)構(gòu)。初始時刻，該ZND爆轟波的von Neumann點位于處，右側(cè)的未燃物狀態(tài)為 狀態(tài)方程及Arrhenius模型中各參數(shù)分。初值的周期擾動在方向，初始流動狀態(tài)其中為在上述參數(shù)下的精確CJ-ZND解。物理區(qū)域為。網(wǎng)格數(shù)為，CFL數(shù)取0.18，上下邊界采用反射邊界條件，左右邊界分別采用入流和出流邊界條件。圖3和圖4分別是元RKDG方法和五階WENO方法結(jié)合Strang分裂法的結(jié)果。由于

24、所用的網(wǎng)格足夠密，二者的結(jié)果差別看不出。在兩種方法得到的圖中，都可以清楚的看到三波點沿著爆轟波陣面的橫向移動[17]。 圖3 P2-RKDG結(jié)合Strang分裂模擬二維不穩(wěn)定爆轟波結(jié)果，從左到右、從上到下依次為以Dt=1/30為時間間隔，從t=3/30到8/30的密度等值線圖 Fig.3 Numerical results of 2d unstable detonation wave by P2-RKDG with Strang split. Sequences from left to right and top to bottom correspond t

25、o those from t=3/30 to 8/30 with timestep Dt=1/30 圖4 五階FVWENO結(jié)合Strang分裂模擬二維不穩(wěn)定爆轟波結(jié)果， 從左到右、從上到下依次為以Dt=1/30為時間間隔，從t=3/30到8/30的密度等值線圖 Fig.4 Numerical results of 2d unstable detonation wave by fifth order FVWENO with Strang split. Sequences from left to right and top to bottom correspon

26、d to those from t=3/30 to 8/30 with timestep Dt=1/30 2.5 CJ-ZND爆轟波繞過拐角的衍射現(xiàn)象 第五個算例為一個CJ-ZND爆轟波繞過一個的拐角，計算區(qū)域為，其中是一固壁，初始時刻一向右傳播的CJ-ZND爆轟波位于固壁上面處。右側(cè)的未燃氣體狀態(tài)為狀態(tài)方程及Arrhenius模型中的參數(shù)分別為 計算取到時間，CFL數(shù)取為0.18，網(wǎng)格數(shù)為。上下邊界和壁面采用反射邊界條件，左右邊界分別采用入流和出流邊界條件。 當爆轟波傳播到拐角后，支撐爆轟波燃燒的激波在壁面附近的強度減弱，溫度降低，沿著繞射激波，激波后面的溫度將會明顯地降低到比相應(yīng)Z

27、ND波的溫度低，這時化學反應(yīng)將會停止。隨著未反應(yīng)激波垂直于豎直墻的表面繼續(xù)往前傳播，緊跟其后的流體溫度將會上升，這時依賴于活化能的強度，如果波后溫度相對于活化能足夠高，燃料可被重新點燃，爆燃波被重新激發(fā)，并在遠處與水平向右傳播的爆轟波相連接。圖5(a)和5(b)分別是用階RKDG和五階WENO結(jié)合Strang分裂得到的結(jié)果，可以清楚的看到爆轟波繞過拐角后爆燃波和前導(dǎo)激波分離的現(xiàn)象。我們將FVWENO在網(wǎng)格上的計算結(jié)果作為更準確的參考解，并取的截面進行了比較，發(fā)現(xiàn)在網(wǎng)格上RKDG的結(jié)果和參考解的符合程度要略好于FVWENO的結(jié)果。 (a) RKDG方法的密度等值線圖 (a) Dens

28、ity contour of P2 RKDG method (b) FVWENO方法的密度等值線圖 (b) Density contour of 5th order FVWENO 圖5 CJ-ZND爆轟波繞過90°拐角 Fig. 5 CJ-ZND detonation wave passing a 90° corner 圖6 x=25截面的密度 Fig. 6 Density of slice plane x=25 3 結(jié)論 本文采用龍格庫塔間斷有限元方法對帶有源項的守恒律方程組進行了數(shù)值求解并和有限體積型WENO方法進行了比較。數(shù)值結(jié)果表明，間斷有限元方法

29、不僅在光滑的區(qū)域具有較高精度，而且具有非常好的間斷捕捉能力，因而可以更好的模擬爆轟波的ZND結(jié)構(gòu)。即使對于強剛性問題，由于RKDG法在單元內(nèi)通過選取適當?shù)幕瘮?shù)提高數(shù)值精度，能夠在較少的網(wǎng)格數(shù)下得到正確的爆轟波位置，克服了一般有限體積方法在求解強剛性問題時由于網(wǎng)格分辨率不夠而容易產(chǎn)生虛假反應(yīng)的現(xiàn)象，而且由于RKDG方法不要求近似解的全局連續(xù)性，因此非常易于自適應(yīng)和并行計算，有望成為模擬復(fù)雜化學反應(yīng)流動的一種有效方法。 參考文獻 [1] Tosatto L, Vigevano L. Numerical solution of under-resolved detonations [J]. J

30、 Comput Phys, 2008, 227: 2317-2343. [2] Helzel C, LeVeque R, Warneke G. A modified fractional step method for the accurate approximation of detonation waves [J], SIAM J Sci Comput, 1999, 22: 1489-1510. [3] Berkenbosch A C. Capturing detonation waves for reactive Euler equations: [Ph.D. thesis][D].

31、 Edndhoven, The Netherlands. Technische Universiteit Eindhoven. 1995. [4] Colella P, Majda A, Roytburd V. Theoretical and numerical structure for reacting shock waves [J]. SIAM J Sci Stat Comput, 1986, 7: 1059-1080. [5] LeVeque R J, Yee H C. A study of numerical methods for hyperbolic conservation

32、 laws with stiff source term [J]. J Comput Phys, 1990, 86: 187-210. [6] Oran E S, Boris J P. Numerical simulation of reactive flow [M]. Elsevier, New York, 1987, Chap.13. [7] Yuan L, Tang T. Resolving the shock-induced combustion by an adaptive mesh redistribution method [J]. J Comput Phys, 2007,

33、224: 587-600. [8] Azarenok B, Tang T. Second-order Godunov scheme for reactive flow calculations on moving meshes [J]. J Comput Phys, 2005, 206: 48-80. [9] Deiterding R. Parallel adaptive simulation of multi-dimensional detonation structures. PhD thesis [D], Brandenburgische Technische Universit?t

34、 Cottbus, Sep. 2003. [10] Cockburn B. Discontinuous Galerkin methods for convection-dominated problem [A]. //Barth T J, Deconinck H. High-Order Methods for Computational Physics, Vol. 9, LNCSE[C], Berlin: Springer, 1999: 69-224. [11] Gottlieb S, Shu C W, Tadmor E. Strong stability-preserving high-

35、order time discretization methods [J]. SIAM REVIEW, 2001, 43: 89-112. [12] LeVeque R J. Finite volume methods for hyperbolic problems [M], Cambridge: Cambridge University Press, 2002. [13] Brown P N, Byrne G D, Hindmarsh A C. VODE: A variable coefficient ODE solver [J]. SIAM J. Sci. Stat. Comput,

36、1989, 10: 1038-1051. [14] Shu C W. High order ENO and WENO schemes [A]. //Barth T J, Deconinck H. High-order Methods for Computational Physics, Vol. 9, LNCSE[C], Berlin: Springer, 1999: 439-582. [15] Shi J, Hu C Q, Shu C W. A technique of treating negative weights in WENO schemes [J]. J Comput Phy

37、s, 2002, 175: 108-127. [16] Bao W, Jin S. The random projection method for stiff detonation capturing [J]. SIAM J Sci Comput, 2001, 22: 1000-1026. [17] Lian Y S, Xu K. A gas-kinetic scheme for multimaterial flows and its application in chemical reactions [J]. J Comput Phys, 2000, 163: 349-375. 作者簡介： 1. 張磊(1977-)，男，北京，講師，在讀博士生，主要從事間斷有限元方法及在流動計算中應(yīng)用的研究。通信地址：中國礦業(yè)大學(北京)理學院，北京100083 第一作者照片： 其他： ------------------------------------------ 身份證號，130603197704170319 13