馮希金, 朱光苗, 危銀濤, 賀炅皓

（1.清華大學(xué) 汽車工程系, 北京 100084; 2.青島大學(xué) 應(yīng)用技術(shù)學(xué)院, 山東 青島 266061; 3. 大冢材料科技（上海）有限公司, 上海 200233）

全鋼子午胎滾動阻力的仿真計(jì)算和實(shí)驗(yàn)對比分析

馮希金1, 朱光苗2, 危銀濤1, 賀炅皓3

（1.清華大學(xué) 汽車工程系, 北京 100084; 2.青島大學(xué) 應(yīng)用技術(shù)學(xué)院, 山東 青島 266061; 3. 大冢材料科技（上海）有限公司, 上海 200233）

簡述了橡膠材料的滯后特性和輪胎滾動阻力產(chǎn)生的機(jī)理，推導(dǎo)了輪胎滾動阻力與輪胎單周滾動損失之間的關(guān)系，從而確定了輪胎滾動阻力的仿真方法。在此基礎(chǔ)上，探討了某規(guī)格全鋼子午胎的滾動阻力受氣壓、實(shí)驗(yàn)速度和載荷的影響規(guī)律，并將仿真結(jié)果與實(shí)驗(yàn)結(jié)果進(jìn)行了詳細(xì)對比。從對比結(jié)果看本文所述的仿真方法是可靠的，同時(shí)也指出了滾動阻力仿真需要進(jìn)一步優(yōu)化的方向。

輪胎; 滾動阻力; 有限元仿真

0 前 言

輪胎的滾動阻力是汽車行駛阻力的重要組成部分，由于輪胎的滾動阻力而導(dǎo)致的汽車燃油消耗占總消耗的20%左右[1]。在進(jìn)入能源緊缺時(shí)代后，降低輪胎的滾動阻力以降低汽車燃油消耗，并隨之降低二氧化碳等溫室氣體排放成為輪胎企業(yè)的重要任務(wù)。尤其是近些年來，汽車制造企業(yè)、政府機(jī)構(gòu)、社會公眾對輪胎的滾動阻力的要求越來越高，歐盟的輪胎標(biāo)簽法案也正式把滾動阻力作為一個重要指標(biāo)[2]，降低輪胎滾動阻力，提高產(chǎn)品的滾動阻力等級顯得更為緊迫。而開展輪胎的滾動阻力仿真和實(shí)驗(yàn)研究，以實(shí)驗(yàn)結(jié)果來驗(yàn)證和完善仿真方法，從而為低滾動阻力輪胎的開發(fā)提供可靠的仿真工具是一個非常有意義的工作。

在輪胎的滾動阻力實(shí)驗(yàn)方面已經(jīng)有了規(guī)范的實(shí)驗(yàn)標(biāo)準(zhǔn)[3,4]，在此不再贅述。輪胎的滾動阻力仿真方面，多年來國內(nèi)外有許多人進(jìn)行了大量的研究[5-10]。這些研究中有些是基于自主編程的方法進(jìn)行輪胎的有限元仿真[11,12]，然后在此基礎(chǔ)上進(jìn)行輪胎的滾動阻力計(jì)算；有些是基于Marc仿真軟件進(jìn)行輪胎滾動阻力的二次開發(fā)[13]；有些是基于Ansys和ls-dyna仿真軟件進(jìn)行二次開發(fā)[14]。通過這些研究工作，使輪胎滾動阻力的仿真趨于成熟。

本文基于Abaqus軟件進(jìn)行輪胎的穩(wěn)態(tài)自由滾動計(jì)算，得到輪胎在穩(wěn)態(tài)滾動下的應(yīng)力和應(yīng)變場，然后在此基礎(chǔ)上編制了滾動阻力仿真計(jì)算軟件，用于子午胎滾動阻力仿真分析。為了驗(yàn)證仿真結(jié)果的可靠性，進(jìn)行了輪胎在多個狀態(tài)下的滾動阻力實(shí)驗(yàn)，將兩者進(jìn)行了對比驗(yàn)證分析，重點(diǎn)剖析輪胎滾動阻力仿真中存在的應(yīng)該進(jìn)一步改進(jìn)的問題和研究方向。

1 橡膠材料的滯后特性和輪胎滾動阻力的機(jī)理

1.1 橡膠材料的滯后特性

橡膠材料是一種黏彈性材料，其突出的特點(diǎn)是滯后性。滯后現(xiàn)象被定義為：橡膠試樣在交變應(yīng)力作用下，應(yīng)力的變化落后于應(yīng)變的變化的現(xiàn)象。

如果把橡膠產(chǎn)生的動態(tài)應(yīng)變ε看作是角頻率為ω、振幅為ε0的正弦波，則動態(tài)應(yīng)變可用下式來描述：

如果橡膠是純粹的彈性體，其應(yīng)力應(yīng)變將同相位變化。但是由于橡膠材料黏彈性的存在，應(yīng)力是具有相位滯后角δ的正弦波，可用下式來表示：

將上式展開，可以得到如下公式：

式（4）中的右邊第一項(xiàng)表示與應(yīng)變ε同相位的部分，右邊第二項(xiàng)表示與應(yīng)變ε正交的部分。σst和σl則分別是整個應(yīng)力中存儲部分和損耗部分所作的貢獻(xiàn)。很容易給出E'/E''與滯后相位δ的關(guān)系：

tanδ稱為損耗正切。

1.2 輪胎滾動阻力的機(jī)理和表現(xiàn)

根據(jù)車輛動力學(xué)的基本理論和馮希金的相關(guān)分析[15]，車輛在行駛中其從動輪的受力如圖1所示。Fv2是輪胎所受車輛的垂直載荷，F(xiàn)p2是輪胎所受車軸的拉力，F(xiàn)x2是地面作用到輪胎上的切向力。Fz2是輪胎所受地面垂直反力，G2是輪胎的重力，Tj2是輪胎旋轉(zhuǎn)慣性阻力矩，Tf2是由于輪胎滾動阻力引起的滾動阻力矩，Ω2為旋轉(zhuǎn)角速度。

圖1 行駛過程中從動輪的受力分析示意圖

如果輪胎是完全剛性體，那么在旋轉(zhuǎn)過程中其所受體面垂直反力Fz2應(yīng)該與重力和Fv2的合力在同一條垂線上，但是由于輪胎的黏彈性，它們并不在同一條直線上，而是偏移一個距離e2，從而產(chǎn)生了滾動阻力距Tf2。輪胎水平方向的力平衡方程式：

也就是說輪胎處于自由滾動狀態(tài)時(shí)地面切向力就等于滾動阻力。將公式（8）右側(cè)分子分母同乘以2π，可以得到公式：

其中Eloss=2πTf2為輪胎滾動一周所損失的能量。

從公式（9）看到，輪胎的滾動阻力就等于輪胎滾動單位長度所損失的能量。

簡諧應(yīng)變情況下單位體積的橡膠材料在每一個周期中的能量損耗為：

將式（1）和式（2）代入上式可得：

將上式展開，并積分可得：

則一個周期內(nèi)輪胎的滾動能量損耗為：

其中：?Wi為一個周期內(nèi)輪胎各部分的能量損耗，Vi為輪胎各部分的體積。

綜合以上各式即可進(jìn)行輪胎滾動阻力的計(jì)算，最終輪胎滾動阻力可表示為：

由公式（14）可以看到，只要求出了輪胎滾動一周過程中的能量損失，就可以求出其滾動阻力。而根據(jù)公式（12）和（13），知道了輪胎滾動一周過程中的應(yīng)力和應(yīng)變，在此基礎(chǔ)上對應(yīng)力和應(yīng)變進(jìn)行傅立葉變換，將周期函數(shù)變成若干個簡諧函數(shù)的和，就可以根據(jù)公式（12）求出?W，進(jìn)而求出ELoss，所以問題轉(zhuǎn)化成求輪胎滾動過程中的應(yīng)力和應(yīng)變場。輪胎的應(yīng)力應(yīng)變場是采用有限元方法求解的。

2 輪胎滾動狀態(tài)應(yīng)力和應(yīng)變場的求解

2.1 輪胎的有限元模型

根據(jù)該規(guī)格輪胎的材料分布圖，可以作出其二維有限元網(wǎng)格，如圖2所示。本二維有限元模型共有13種橡膠單元選擇集，8種骨架材料單元選擇集，共有1903個節(jié)點(diǎn)、1749個單元。橡膠的二維單元類型是CGAX4H（四邊形單元）和CGAX3H（三角形單元），骨架材料rebar單元的類型是SFMGAX1。

圖2 二維有限元網(wǎng)格圖

橡膠材料模型采用Yeoh材料模型，此為超彈性材料。橡膠的滯后特性用tanδ表示，由DMA動態(tài)黏彈譜儀測量得到。

三維有限元模型（如圖3所示）周向等角度劃分120份，保證了輪胎應(yīng)力和應(yīng)變場的計(jì)算精度，從而保證可以提供良好的滾動阻力計(jì)算結(jié)果。

2.2 滾動狀態(tài)應(yīng)力和應(yīng)變場的求解

首先進(jìn)行輪胎的二維安裝輪輞和充氣分析，然后采用Abaqus的*SYMMETRIC MODEL GENERATION，REVOLVE關(guān)鍵詞，將二維的充氣狀態(tài)空間旋轉(zhuǎn)成為三維的充氣狀態(tài)。然后給輪胎施加標(biāo)準(zhǔn)載荷，得到輪胎的靜態(tài)加載狀態(tài)。在此基礎(chǔ)上，給輪胎施加旋轉(zhuǎn)速度，得到輪胎的自由滾動狀態(tài)，自由滾動狀態(tài)下的應(yīng)力和應(yīng)變場就是我們用來進(jìn)行滾動阻力計(jì)算的應(yīng)力和應(yīng)變場。

圖3 輪胎三維有限元模型圖

3 滾動阻力的分析計(jì)算及其與實(shí)驗(yàn)結(jié)果的對比

得到輪胎的穩(wěn)態(tài)自由滾動狀態(tài)后，可以提取出單元的應(yīng)力和應(yīng)變，然后對所有單元的應(yīng)力和應(yīng)變進(jìn)行處理，得到輪胎的滾動阻力和滾動阻力系數(shù)。下面即為在各種工況條件下12R22.5輪胎的滾動阻力計(jì)算結(jié)果，并同時(shí)給出了實(shí)驗(yàn)結(jié)果，便于兩者對比分析。

3.1 氣壓對輪胎滾動阻力的影響

給輪胎施加36 138N的載荷，輪胎自由滾動速度為110km/h，輪胎的氣壓分別為860kPa、990kPa和1120kPa三種情況。之所以選擇這三個工況，是因?yàn)閷?shí)驗(yàn)報(bào)告中只有在110km/h時(shí)的情況下才給出了具體的氣壓值，而氣壓決定著輪胎的變形和剛度，所以我們選擇這個速度下的氣壓，以便更好地模擬實(shí)驗(yàn)條件。

計(jì)算得到的滾動阻力和滾動阻力系數(shù)與實(shí)驗(yàn)結(jié)果的對比，如表1所示。

從表1看出，從仿真結(jié)果來說，滾動阻力隨著氣壓的增大而減小，這與理論是相符的，因?yàn)闅鈮涸酱?，輪胎剛度越大，各部分的變形減小，尤其是輪胎滾動一周的過程中，其應(yīng)力應(yīng)變狀態(tài)變化減小，那么導(dǎo)致的滾動損失就會減小，這也和人們的用胎常識一致。

表1 載荷36 138N、速度110km/h條件下氣壓對滾動阻力的影響

從實(shí)驗(yàn)結(jié)果來看，隨著氣壓的增加，滾動阻力也是減小的。但在110km這個速度條件下，氣壓990kPa和1120kPa時(shí)滾動阻力相差很小，認(rèn)為這是由于實(shí)驗(yàn)誤差所導(dǎo)致的。進(jìn)一步對實(shí)驗(yàn)過程中其他速度條件下的滾動阻力隨氣壓的變化進(jìn)行統(tǒng)計(jì)分析，得到如表2所示的結(jié)果。

表2 實(shí)驗(yàn)條件下氣壓對滾動阻力的影響

表2中給出了輪胎在不同的初始?xì)鈮?、不同的載荷下以40km/h跑30min后的滾動阻力測試結(jié)果。第一列是氣壓，第二列是載荷為36 138N的結(jié)果，第三列和第四列分別是載荷為32 852N和26 282N的結(jié)果。之所以選擇40km/h跑30min后的結(jié)果來對比，是因?yàn)檩喬L動的過程中氣壓會升高，其他速度是在40km/h的速度基礎(chǔ)上得到的，因此它們的氣壓是無法確知的，只有這個初始的速度條件下才最接近初始?xì)鈮?。從?看到，隨著氣壓的升高，滾動阻力是下降的。

從表1的滾動阻力系數(shù)來看，氣壓越高，滾動阻力系數(shù)越小。從實(shí)驗(yàn)結(jié)果和仿真結(jié)果的對比來看，滾動阻力系數(shù)的誤差在10%以下，還是比較可靠的。

3.2 速度對輪胎滾動阻力的影響

在討論之前，我們首先說明一下實(shí)驗(yàn)測量滾動阻力的方法和仿真的條件設(shè)置。在實(shí)驗(yàn)報(bào)告里，輪胎滾動阻力實(shí)驗(yàn)時(shí)首先讓輪胎預(yù)熱，然后冷卻6h，再進(jìn)行正常的實(shí)驗(yàn)。首先以40km/h的速度跑90min，測量滾動阻力，然后不經(jīng)冷卻，速度設(shè)置成60km/h跑30min，測量此時(shí)的滾動阻力。然后速度再設(shè)置成80km/h跑30min，測量滾動阻力。再接著以100km/h跑30min，測量滾動阻力；最后以110km/h跑30min，測量滾動阻力。初始的氣壓是830kPa。但這么連續(xù)實(shí)驗(yàn)下來，到110km/h結(jié)束時(shí)胎內(nèi)氣壓達(dá)到了980kPa，整個過程氣壓升高了整整150kPa，可見輪胎在滾動過程中產(chǎn)生了大量熱量。溫度升高，導(dǎo)致胎內(nèi)的氣壓升高。胎內(nèi)氣壓升高，必然使輪胎的剛度增大，也會使橡膠材料的tanδ發(fā)生變化。而實(shí)驗(yàn)報(bào)告只給出了初始的充氣壓力830kPa和實(shí)驗(yàn)結(jié)束時(shí)的充氣壓力980KPa，也就是我們不知道40km/h、60km/h、80km/h、100km/h結(jié)束時(shí)的輪胎內(nèi)部氣壓。

依據(jù)目前的仿真計(jì)算流程，我們還無法實(shí)現(xiàn)完全的熱-機(jī)耦合仿真，也就是在仿真過程中還無法實(shí)現(xiàn)如下圖所示的閉環(huán)仿真：

圖4 熱-機(jī)耦合閉合仿真

當(dāng)前我們的仿真技術(shù)是先給定氣壓、載荷、速度，然后求解輪胎的穩(wěn)態(tài)自由滾動狀態(tài)下的應(yīng)力應(yīng)變，根據(jù)這個應(yīng)力應(yīng)變，采用能量法求解滾動阻力。沒有考慮滾動中的生熱造成的氣壓和物理性能的改變。而實(shí)驗(yàn)過程中速度是從40km/h逐步遞增的，中間的速度沒有測量其氣壓值。所以為了使仿真的結(jié)果更好地反映實(shí)際情況，我們在進(jìn)行仿真時(shí)認(rèn)為40、60、80km/h這三個速度實(shí)驗(yàn)完畢時(shí)氣壓仍然為830kPa，110km/h的條件下我們采用實(shí)測的氣壓980kPa，110km/h下的氣壓我們采用950kPa。用這些條件去模擬實(shí)驗(yàn)條件，得到的實(shí)驗(yàn)滾動阻力和仿真滾動阻力隨速度的關(guān)系如表3所示。

從表3看到，隨著速度的增加，仿真得到的滾動阻力稍微增加一些，而實(shí)驗(yàn)得到的滾動阻力也是稍微增加一些，在100km/h和110km/h的速度下，仿真結(jié)果大幅下降，原因是氣壓比上述3個速度增加了100kPa多。實(shí)驗(yàn)得到的滾動阻力也有這種現(xiàn)象，只是下降幅度沒有仿真多。

從表3中還可以看到，不同的速度下，仿真得到的滾動阻力系數(shù)和實(shí)驗(yàn)得到滾動阻力系數(shù)基本一致，誤差的絕對值在10%以下，說明仿真結(jié)果還是比較可靠的。

3.3 載荷對滾動阻力的影響

在實(shí)驗(yàn)報(bào)告中，分別給出了在26 282N、32 852N、36 138N三個載荷在速度為110km/h的氣壓測量值，我們進(jìn)行仿真分析時(shí)也以這個載荷、氣壓和速度為基礎(chǔ)進(jìn)行仿真。結(jié)果對比如表4所示。

表3 不同速度下滾動阻力實(shí)驗(yàn)結(jié)果和仿真結(jié)果的對比

表4 速度為110km/h的情況下，載荷對滾動阻力的影響

從表4中看到，仿真計(jì)算得到的滾動阻力比較小，而實(shí)驗(yàn)得到的滾動阻力比較大。其中氣壓為26 282N時(shí)誤差較大，估計(jì)原因是在很大的氣壓、較小的載荷下，輪胎剛度比較大。由于仿真所采用的是只有縱溝的光面輪胎，本身剛度就比實(shí)際輪胎大，再加上大氣壓，造成其仿真計(jì)算的輪胎的剛度進(jìn)一步增加，從而導(dǎo)致滾動阻力減小。而實(shí)際的輪胎，由于有橫向花紋溝，因此其胎面剛度要比只有縱溝的光面胎要小，變形大，因此實(shí)驗(yàn)得到滾動阻力比仿真的大。

這就說明，建模時(shí)由于沒有考慮花紋溝，只把輪胎當(dāng)作帶有縱溝的光面胎，會使輪胎剛度比實(shí)際剛度大，從而得到較小的滾動阻力。因此，為了進(jìn)一步提高仿真精度，應(yīng)該考慮帶花紋輪胎的滾動阻力仿真。

4 結(jié) 論

通過上面的分析和對比，我們可以得出如下的結(jié)論：

（1）目前開發(fā)的輪胎滾動阻力仿真技術(shù)可以實(shí)現(xiàn)輪胎各種條件下的滾動阻力計(jì)算，在正常的載荷和氣壓范圍內(nèi)，計(jì)算結(jié)果和仿真結(jié)果比較接近，誤差在10%以下。

（2）如果是大氣壓、小載荷的受力條件，由于仿真沒有考慮花紋的影響，輪胎剛性大為增加，輪胎變形減小，雖然輪胎的受力有可能增大，但是由于輪胎轉(zhuǎn)動一周的過程中應(yīng)力和應(yīng)變的變化幅值比較小，因此其滾動阻力要減小。而實(shí)驗(yàn)結(jié)果，由于真實(shí)輪胎是帶有橫向花紋溝的，所以其胎面剛度受氣壓的影響小，測量得到的滾動阻力仍會較大，這就會導(dǎo)致仿真和實(shí)驗(yàn)之間的差別變大。

（3）當(dāng)前的仿真，大部分是基于光面輪胎，或者是帶有縱溝的光面輪胎，但是實(shí)際的輪胎都是有復(fù)雜花紋的，花紋的溝槽導(dǎo)致胎面橡膠塊剛度的變化，相同的條件下，胎面橡膠塊的剛度要比光面輪胎的剛度小，剛度是要影響胎面橡膠的變形和受力的，因此這種剛度的差別會造成仿真的滾動阻力和實(shí)測的滾動阻力有差別。

（4）從計(jì)算原理上說，我們提取的應(yīng)力和應(yīng)變是基于輪胎的穩(wěn)態(tài)自由滾動狀態(tài)的，從數(shù)學(xué)上說，穩(wěn)態(tài)自由滾動狀態(tài)的輪胎沒有考慮動力學(xué)的影響。因?yàn)檩喬L動起來后，胎面橡膠是不斷經(jīng)歷著從圓弧狀態(tài)到平直狀態(tài)這個變化的，也就是路面會給胎面施加沖擊作用，這個沖擊作用導(dǎo)致輪胎的振動，如果不考慮動力學(xué)的因素，那么這個振動的激勵就考慮不到，而振動激勵也會導(dǎo)致輪胎的額外的變形。也就是說，仿真所用的穩(wěn)態(tài)自由滾動輪胎沒有考慮到真實(shí)輪胎的動力學(xué)效應(yīng)，因此從這點(diǎn)上將也會有誤差。

（5）輪胎的滾動阻力測量是在一個準(zhǔn)穩(wěn)態(tài)的條件下進(jìn)行的，正如圖4所示的那樣，在實(shí)驗(yàn)的過程中輪胎溫度升高，導(dǎo)致氣壓和物理性能發(fā)生變化，而這些變化又反過來影響輪胎的受力狀態(tài)，而滾動阻力是基于受力狀態(tài)的。仿真時(shí)沒有能夠完成這個閉環(huán)的仿真，所以有可能使仿真結(jié)果和實(shí)驗(yàn)結(jié)果有一定的差距。不過從我們實(shí)際情況看，誤差并不大。

總之，本文基于穩(wěn)態(tài)自由滾動的受力狀態(tài)所得的滾動阻力與實(shí)測的滾動阻力誤差在10%以下，可以提供比較好的可靠度，可以用來進(jìn)行滾動阻力的預(yù)報(bào)研究。

[1] Michelin. Tire Rolling Resistance, 2001.15-23.

[2] On the labelling of tyres with respect to fuel efficiency and other essential parameters, Regulation (EC) No 1222/2009 ofthe European Parliament and of the Council of 25 November 2009.

[3] GB/T 18861-2012. 汽車輪胎和摩托車輪胎滾動阻力試驗(yàn)方法 多點(diǎn)試驗(yàn)[S].

[4] ISO 28580.Passenger car, truck and bus tyres -- Methods of measuring rolling resistance -- Single point test and correlation of measurement results[S].

[5] Walter J D , Conant F S. Energy losses in tires[J]. Tire Science and Technology, 1974, 2(04)

[6] Shida Z, Koishi M, Kogure T, Kabe K. A rolling resistance simulation of tires using static finite element analysis[J]. Tire Science and Technology, 1999, 27(02).

[7] Clark S K. Rolling resistance of pneumatic tires[J]. Tire Science and Technology, 1978, 6(03).

[8] 葉進(jìn)雄. 全鋼載重子午胎滾動阻力有限元仿真與試驗(yàn)分析[D]. 清華大學(xué)碩士學(xué)位論文, 2007.

[9] 丁劍平, 賈德民, 周小舒, 等.子午線輪胎滾動阻力的數(shù)值模擬[J]. 汽車技術(shù), 2005(03)：26-28.

[10] 王國林, 伍建軍. 載重子午線輪胎滾動阻力有限元仿真分析[J]. 計(jì)算機(jī)仿真, 2001(11)：350-352.

[11] 危銀濤. 輪胎熱力學(xué)分析及耐久性評價(jià)[D]. 哈爾濱工業(yè)大學(xué)博士學(xué)位論文, 1997.

[12] 危銀濤, 劉宇艷, 杜星文, 吳寶國. 子午線輪胎滾動阻力與溫度場非線性有限元分析[J]. 輪胎工業(yè),1998(06)：10-15.

[13] 丁劍平, 賈德民, 俞 淇. 用有限元法分析輪胎結(jié)構(gòu)與滾動阻力的關(guān)系[J]. 廣西大學(xué)學(xué)報(bào)(自然科學(xué)版), 2007,32(4)：81-85.

[14] 韓秀枝. 子午線輪胎穩(wěn)態(tài)滾動阻力及水滑特性的研究[D]. 北京化工大學(xué),2009.

[15] 馮希金, 危銀濤, 鄭小剛, 等，論輪胎的自由滾動力學(xué)[J], 橡膠工業(yè), 2013,(9)：6-15.

[責(zé)任編輯：朱 胤]

FEA Simulation and Testing Comparison for the Rolling Resistance of All-steel Radial Tire

Feng Xijin1, Zhu Guangmiao2, Wei Yintao1, He Jionghao3
(1.Tsinghua University, Beijing 100084; 2. TTA (Qingdao) Tire Technology Co., Ltd, Qingdao 266061; 3.Otsuka Material Science and Technology(Shanghai) Co., Ltd, Shanghai 200233)

In this article, the hysteresis of rubber compounds and mechanism of tire rolling resistance are described. And the relationship of tire rolling resistance and single-turn rolling hysteresis has been derived to simulate the tire rolling resistance. Based on this, the infl uence of pressure, testing speed and load on the rolling resistance of a certain size of truck tire was studied. We found the FEA method to calculate tire rolling resistance is reliable when compared with tire testing results. Also it is suggested for the future improvement of the FEA simulation on rolling resistance.

Tire; Rolling Resistance; FEA Simulation

TQ 336.1

B

1671-8232(2014)04-0020-07

2014-02-07