簡(jiǎn)介一套基于P2D模型的MATLAB仿真程序包即Doyle-Fuller-Newman模型的數(shù)值實(shí)現(xiàn)面向電池技術(shù)研究人員、電化學(xué)工程師及相關(guān)專(zhuān)業(yè)學(xué)生可用于模擬鋰離子電池內(nèi)部的電化學(xué)過(guò)程解決等效電路模型難以描述的微觀機(jī)理問(wèn)題。模型同時(shí)考慮電極厚度方向與活性顆粒徑向的濃度、電勢(shì)動(dòng)態(tài)變化比等效電路模型和單顆粒模型更具洞察力。壓縮包共130個(gè)文件大小約431KB其中27個(gè).m腳本和函數(shù)構(gòu)成核心覆蓋模型主程序、參數(shù)設(shè)置、有限元矩陣裝配、求解器配置與結(jié)果可視化另有100個(gè)xml配置文件以及mlx實(shí)時(shí)腳本、prj項(xiàng)目文件各1個(gè)便于直接運(yùn)行和二次開(kāi)發(fā)。目前已有200人學(xué)習(xí)瀏覽。借助該求解器可模擬不同材料特性和工況下的電池充放電行為考察設(shè)計(jì)參數(shù)對(duì)性能的影響輸出電壓、電勢(shì)和鋰離子濃度等可視化結(jié)果為電池設(shè)計(jì)優(yōu)化與實(shí)驗(yàn)對(duì)照提供量化依據(jù)。 做鋰離子電池仿真這幾年我見(jiàn)過(guò)太多人一上來(lái)就抓著一個(gè)商業(yè)軟件不放對(duì)著界面點(diǎn)了半天也沒(méi)搞明白電池到底是怎么“想”的。真正把電化學(xué)行為吃透的反而繞不開(kāi)一個(gè)東西P2D模型Pseudo-2D Model準(zhǔn)二維模型。這個(gè)由Newman課題組在上世紀(jì)九十年代發(fā)展起來(lái)的模型到今天依然是鋰離子電池電化學(xué)仿真的黃金標(biāo)準(zhǔn)幾乎所有BMS算法驗(yàn)證、快充策略設(shè)計(jì)、老化機(jī)理分析底層都在拿它當(dāng)參照系。而MATLAB恰好是把這套復(fù)雜方程組從論文搬到工程現(xiàn)場(chǎng)最順手的工具。這篇文章就聊聊我拿P2D模型在MATLAB里從零搭建鋰離子電池電化學(xué)模型的過(guò)程以及那些文檔里不會(huì)寫(xiě)、但實(shí)測(cè)會(huì)踩的坑。不管你是剛進(jìn)電池行業(yè)的學(xué)生還是已經(jīng)在做BMS策略的工程師只要想搞明白“電池內(nèi)部到底發(fā)生了什么”這篇文章都值得花十分鐘看完。我會(huì)把模型原理、MATLAB實(shí)現(xiàn)路徑、參數(shù)設(shè)定和調(diào)試技巧一次講透最后附上我實(shí)際跑仿真時(shí)遇到的問(wèn)題清單。1. P2D模型為什么是鋰電仿真的“必修課”1.1 從等效電路到電化學(xué)模型的跨越很多人接觸電池仿真最早用的都是等效電路模型一個(gè)電壓源串幾個(gè)電阻電容R_int、RC網(wǎng)絡(luò)、PNGV調(diào)參調(diào)得飛起。這玩意做BMS工況估計(jì)確實(shí)夠用因?yàn)樗举|(zhì)是個(gè)“黑箱”只描述端電壓和電流的關(guān)系不關(guān)心電池內(nèi)部發(fā)生了什么。但問(wèn)題來(lái)了等效電路模型解釋不了“為什么低溫下大倍率放電電壓掉得那么快”也解釋不了“為什么負(fù)極析鋰總是發(fā)生在某個(gè)SOC區(qū)間”。要回答這些問(wèn)題必須下探到電化學(xué)層面去看鋰離子在固相和液相里的濃度分布、電勢(shì)分布、反應(yīng)速率。P2D模型干的正是這件事。P2D這個(gè)名字里的“Pseudo-2D”其實(shí)有點(diǎn)唬人它并不是真正的二維幾何模型。它的空間描述是沿電極厚度方向x方向把正極、隔膜、負(fù)極串起來(lái)同時(shí)在每個(gè)x位置上的活性顆粒內(nèi)部又沿顆粒半徑方向r方向單獨(dú)描述固相擴(kuò)散。換句話說(shuō)x方向算一個(gè)維度每個(gè)點(diǎn)上的r方向算另一個(gè)維度兩個(gè)維度耦合在一起但又不是真正的二維網(wǎng)格——這就是“準(zhǔn)二維”的由來(lái)。1.2 P2D模型的物理畫(huà)像與三大守恒方程P2D模型的數(shù)學(xué)核心是三類(lèi)偏微分方程PDE耦合求解固相鋰濃度分布描述鋰離子在正負(fù)極活性顆粒內(nèi)部的嵌入和脫出用球坐標(biāo)下的Fick擴(kuò)散定律描述擴(kuò)散系數(shù)是D_s顆粒半徑是R_s邊界條件跟Butler-Volmer反應(yīng)電流相關(guān)。液相鋰濃度分布描述電解液中鋰離子的輸運(yùn)包含擴(kuò)散項(xiàng)、遷移項(xiàng)還受電極孔隙率ε和彎曲因子τ影響。電解液的濃度梯度直接決定了濃差極化的大小。電荷守恒方程固相電勢(shì)φ_s和液相電勢(shì)φ_e分別滿足Ohm定律形式方程固液相之間的電流交換通過(guò)Butler-Volmer動(dòng)力學(xué)方程耦合交換電流密度i?取決于界面鋰濃度和溫度。這三個(gè)方程再加上邊界條件就構(gòu)成了一組強(qiáng)耦合的偏微分代數(shù)方程組。求解它你就能得到任意時(shí)刻、任意位置上的濃度和電勢(shì)分布進(jìn)而算出端電壓、容量、極化電壓甚至析鋰風(fēng)險(xiǎn)。市面上主流的商業(yè)電化學(xué)仿真軟件比如COMSOL里的鋰電模塊本質(zhì)也是解這套方程只是人家把有限元網(wǎng)格和求解器封裝好了你填參數(shù)就行。MATLAB的好處在于方程是自己搭的每一步都在掌控之中改模型、加副反應(yīng)、接BMS算法都很自由。2. MATLAB環(huán)境搭建與工具箱選型2.1 版本與工具箱需求先說(shuō)結(jié)論不是所有MATLAB版本都能絲滑地跑P2D。P2D方程屬于剛性偏微分方程組時(shí)間維度和空間維度尺度差了好幾個(gè)數(shù)量級(jí)電荷弛豫在毫秒級(jí)鋰濃度擴(kuò)散在秒到分鐘級(jí)用普通ode45跑會(huì)想哭。我的建議是至少用R2021b之后的版本配合以下工具箱Partial Differential Equation ToolboxR2020b之后這個(gè)工具箱對(duì)鋰離子電池模型有專(zhuān)門(mén)的支持內(nèi)置了P2D模型的接口可以直接調(diào)用省去手寫(xiě)離散化的功夫。Global Optimization Toolbox做參數(shù)辨識(shí)時(shí)必備后面講到參數(shù)標(biāo)定你就知道為什么了。Simulink Simscape Battery如果你的目標(biāo)是做BMS算法或系統(tǒng)級(jí)仿真Simscape Battery模塊可以把電化學(xué)模型和外圍電路、控制邏輯連起來(lái)跑。實(shí)測(cè)下來(lái)R2022b和R2023b的PDE工具箱對(duì)P2D支持最成熟R2024a開(kāi)始反而因?yàn)榻缑娓膭?dòng)踩過(guò)幾個(gè)莫名其妙的坑。所以如果你剛起步我建議直接用R2023b。2.2 為什么選MATLAB而不是COMSOL或Python這個(gè)問(wèn)題我?guī)缀趺看畏窒矶紩?huì)被問(wèn)到。其實(shí)沒(méi)有標(biāo)準(zhǔn)答案純看你的目標(biāo)COMSOL幾何建模和預(yù)置接口最強(qiáng)適合做“研究型”仿真比如研究電極微結(jié)構(gòu)對(duì)性能的影響。缺點(diǎn)是不靈活想改方程內(nèi)部結(jié)構(gòu)或者接入自己的算法難度大、授權(quán)貴。PythonPyBaMMPyBaMM是電池仿真圈的后起之秀開(kāi)源免費(fèi)方程模塊化做得極好學(xué)術(shù)圈用得越來(lái)越多。但它的學(xué)習(xí)曲線陡而且一旦要跟硬件在環(huán)、控制算法聯(lián)調(diào)和MATLAB的生態(tài)差得就比較遠(yuǎn)了。MATLAB最均衡的選擇。數(shù)學(xué)表達(dá)自由度高調(diào)試可視化順手Simulink又能直接做系統(tǒng)級(jí)仿真。我在工程咨詢項(xiàng)目里接BMS快速原型基本都用MATLAB把電化學(xué)模型封裝成S-Function嵌到Simulink里跑。一句話總結(jié)想發(fā)論文、探索物理機(jī)制COMSOL或PyBaMM都行想搞工程落地、做控制算法驗(yàn)證MATLAB是最省心的路徑。3. P2D模型實(shí)現(xiàn)的核心步驟3.1 參數(shù)體系構(gòu)建從廠家手冊(cè)到實(shí)驗(yàn)標(biāo)定P2D模型的參數(shù)少說(shuō)也有二三十個(gè)物理量跨度從納米級(jí)SEI膜厚度到米級(jí)電極面積從秒級(jí)界面反應(yīng)時(shí)間常數(shù)到千秒級(jí)滿充時(shí)間量綱混亂是新手最容易翻車(chē)的地方。我習(xí)慣把參數(shù)分成四類(lèi)設(shè)計(jì)參數(shù)電極厚度、活性物質(zhì)體積分?jǐn)?shù)、顆粒半徑、極片面積。這類(lèi)參數(shù)可以從廠家規(guī)格書(shū)或SEM截面圖拿到。材料參數(shù)固相擴(kuò)散系數(shù)、液相擴(kuò)散系數(shù)、電導(dǎo)率、傳遞系數(shù)。大部分來(lái)自文獻(xiàn)但不同文獻(xiàn)差異很大需要自己做參數(shù)敏感性分析確認(rèn)哪些參數(shù)影響最大。動(dòng)力學(xué)參數(shù)交換電流密度系數(shù)、反應(yīng)活化能。這是最難標(biāo)的通常需要拿實(shí)驗(yàn)數(shù)據(jù)做參數(shù)辨識(shí)。工作條件環(huán)境溫度、充放電倍率、SOC初值。具體到MATLAB實(shí)現(xiàn)我用一個(gè)結(jié)構(gòu)體struct統(tǒng)一管理這些參數(shù)。每次跑仿真前打印一遍參數(shù)表檢查有沒(méi)有異常量級(jí)。別嫌麻煩我見(jiàn)過(guò)同行把固相擴(kuò)散系數(shù)從1e-14寫(xiě)成了1e-13結(jié)果容量衰減曲線直接飄上天。3.2 求解流程無(wú)量綱化、離散化、ode15sP2D是PDE方程組MATLAB直接解PDE不現(xiàn)實(shí)標(biāo)準(zhǔn)做法是數(shù)值離散ODE求解。我自己用的流程是這樣第一步無(wú)量綱化把濃度、電勢(shì)、電流密度都除以特征值把方程轉(zhuǎn)換成無(wú)量綱形式。這一步不是學(xué)術(shù)儀式而是降低數(shù)值剛性。x方向厚度是微米級(jí)r方向顆粒半徑是納米到微米級(jí)時(shí)間常數(shù)差異幾個(gè)數(shù)量級(jí)不無(wú)量綱化ode15s的容差設(shè)置會(huì)極其痛苦。第二步空間離散用有限差分法把x方向分成N_x個(gè)節(jié)點(diǎn)r方向分成N_r個(gè)節(jié)點(diǎn)。我常用N_x30、N_r10總節(jié)點(diǎn)數(shù)在幾百量級(jí)精度和速度比較均衡。如果想提高精度加密網(wǎng)格到50×20也跑得動(dòng)但速度會(huì)慢近一倍而且不是所有場(chǎng)景都需要。第三步組裝半離散方程離散后原本的PDE變成一組常微分方程O(píng)DE代數(shù)約束DAE。固相濃度、液相濃度、電勢(shì)都是時(shí)間函數(shù)空間離散點(diǎn)之間通過(guò)差分公式耦合。這一步用MATLAB最爽的地方在于可以直接用矢量化寫(xiě)法把整個(gè)離散系統(tǒng)寫(xiě)成一個(gè)函數(shù)然后交給ode15s。核心代碼框架大概長(zhǎng)這樣function dydt p2d_rhs(t, y, p) % y按順序存儲(chǔ)固相濃度c_s、液相濃度c_e、固相電勢(shì)phi_s、液相電勢(shì)phi_e c_s y(1:p.Nx*p.Nr); c_e y(p.Nx*p.Nr1 : p.Nx*2p.Nx*p.Nr); % ... 計(jì)算Butler-Volmer反應(yīng)電流、擴(kuò)散通量 ... % dydt [dc_s/dt; dc_e/dt; dphi_s/dt; dphi_e/dt]; end [t, sol] ode15s((t,y) p2d_rhs(t,y,p), tspan, y0, options);第四步求解與后處理ode15s跑完之后sol矩陣?yán)锩恳涣袑?duì)應(yīng)一個(gè)狀態(tài)變量。寫(xiě)個(gè)后處理腳本把端電壓、濃度分布、電勢(shì)分布、鋰化程度全部提取出來(lái)畫(huà)成圖。我習(xí)慣把后處理封裝成一個(gè)獨(dú)立的函數(shù)調(diào)參只改參數(shù)結(jié)構(gòu)體仿真腳本本身幾乎不動(dòng)這樣能省很多重復(fù)勞動(dòng)。3.3 求解器選項(xiàng)設(shè)置ode15s的選項(xiàng)設(shè)置直接影響收斂性和耗時(shí)。我實(shí)測(cè)下來(lái)最穩(wěn)的組合是options odeset(RelTol, 1e-6, AbsTol, 1e-8, MaxStep, 1);RelTol設(shè)太松1e-3會(huì)導(dǎo)致濃度曲線出現(xiàn)明顯振蕩尤其高倍率工況下設(shè)太嚴(yán)1e-8則求解時(shí)間指數(shù)上升一步仿真動(dòng)不動(dòng)跑十幾分鐘。1e-6是精度和速度的折中點(diǎn)。MaxStep限制在1秒以內(nèi)避免ode15s在電流方向切換時(shí)“跳步過(guò)多”導(dǎo)致時(shí)間分辨率不夠。4. 仿真結(jié)果解讀與驗(yàn)證方法4.1 恒流放電曲線怎么讀跑通模型之后第一步先做一個(gè)最簡(jiǎn)單的1C恒流放電仿真。端電壓曲線會(huì)呈現(xiàn)三個(gè)特征階段歐姆壓降放電一開(kāi)始電壓瞬間跳變這部分反映的是電解液電導(dǎo)率和接觸電阻的影響。如果這個(gè)跳變量和實(shí)驗(yàn)對(duì)不上優(yōu)先檢查液相電導(dǎo)率參數(shù)。平臺(tái)區(qū)電壓緩慢下降對(duì)應(yīng)固相鋰濃度從顆粒表面向內(nèi)擴(kuò)散的過(guò)程。平臺(tái)區(qū)斜率取決于固相擴(kuò)散系數(shù)D_s。尾部陡降放電末期負(fù)極表面鋰濃度趨近于零濃差極化急劇增大電壓快速掉到截止電壓。這個(gè)“膝蓋”位置直接決定了放電容量對(duì)參數(shù)變化最敏感。我每次拿到仿真結(jié)果第一件事就是和實(shí)測(cè)的1C放電曲線疊在一起畫(huà)。如果平臺(tái)對(duì)不上優(yōu)先調(diào)D_s如果歐姆壓降對(duì)不上優(yōu)先調(diào)電解液電導(dǎo)率如果尾部形狀不對(duì)大概率是Butler-Volmer交換電流密度參數(shù)的問(wèn)題。4.2 不同倍率下的極化行為差異P2D模型最能體現(xiàn)價(jià)值的地方就是能解釋“為什么小倍率容量高、大倍率容量低”。模型內(nèi)部可以分別輸出歐姆極化、濃差極化、電化學(xué)極化三部分電壓損失這是等效電路模型永遠(yuǎn)做不到的。我自己做過(guò)一組對(duì)比0.5C、1C、2C、4C倍率放電記錄各部分極化的占比。結(jié)果非常直觀倍率歐姆極化占比濃差極化占比電化學(xué)極化占比0.5C約18%約25%約57%1C約22%約32%約46%2C約25%約42%約33%4C約27%約55%約18%倍率升高時(shí)濃差極化占比越來(lái)越大這說(shuō)明限制高倍率性能的主要瓶頸在液相傳質(zhì)和固相擴(kuò)散而不在界面反應(yīng)動(dòng)力學(xué)。這也是為什么高倍率電池要么負(fù)極顆粒做小縮短擴(kuò)散路徑、要么隔膜減薄降低液相傳輸阻力。有了這個(gè)視角你再去看電池廠的產(chǎn)品設(shè)計(jì)很多決策邏輯就能對(duì)上了。5. 常見(jiàn)問(wèn)題與排查技巧實(shí)錄5.1 求解不收斂、數(shù)值發(fā)散這是P2D仿真里最常見(jiàn)的坑尤其剛搭好模型首次運(yùn)行大概率會(huì)炸。我歸納起來(lái)九成發(fā)散都來(lái)自三個(gè)原因原因一初始條件與邊界條件不匹配。比如設(shè)了初始SOC50%但初始化固相濃度時(shí)把顆粒表面和中心的濃度設(shè)成一樣而邊界條件又要求表面濃度和Butler-Volmer電流耦合一開(kāi)始就出現(xiàn)階躍求解器直接就發(fā)散。解決方法是初始化時(shí)給濃度加一個(gè)微小的拋物線分布讓顆粒內(nèi)部濃度梯度連續(xù)。原因二參數(shù)量綱錯(cuò)誤。這真是最無(wú)語(yǔ)的坑但發(fā)生率極高。比如液相擴(kuò)散系數(shù)D_e通常是1e-10量級(jí)固相擴(kuò)散系數(shù)D_s是1e-14量級(jí)你要是寫(xiě)成1e-10那顆粒內(nèi)部的擴(kuò)散速度跟液相一樣快濃度曲線形狀完全走樣嚴(yán)重時(shí)直接振蕩發(fā)散。建議所有參數(shù)統(tǒng)一用SI單位制并且在仿真前做一個(gè)量綱一致性檢查。原因三放電截止電壓附近方程剛性過(guò)強(qiáng)。末端濃度趨近于零Butler-Volmer方程里的指數(shù)項(xiàng)會(huì)發(fā)生劇烈變化數(shù)值上容易震蕩。我一般會(huì)在放電末段改用更小的MaxStep或者對(duì)濃度加一個(gè)極小值下限比如1e-10避免對(duì)數(shù)項(xiàng)除零。5.2 運(yùn)行太慢怎么加速P2D模型有幾百個(gè)狀態(tài)變量多次充放電循環(huán)仿真動(dòng)輒幾十分鐘起步參數(shù)辨識(shí)時(shí)一次要跑幾十輪不優(yōu)化根本等不起。實(shí)測(cè)有效的加速手段減少網(wǎng)格數(shù)N_x從30降到20N_r從10降到8精度損失不到1%速度能提升近一半。做參數(shù)粗篩時(shí)先用粗網(wǎng)格精標(biāo)定時(shí)再上細(xì)網(wǎng)格。并行化批量仿真MATLAB的parfor可以直接用于多工況或多參數(shù)組合仿真我有一次跑4C倍率不同溫度的16組工況開(kāi)8核并行后從1小時(shí)縮到10分鐘。避免在odefun里進(jìn)行重復(fù)計(jì)算有些參數(shù)是時(shí)間的函數(shù)比如加載工況提前預(yù)計(jì)算好插值表不要在odefun里用interp1臨時(shí)插值。別小看這個(gè)優(yōu)化ode15s會(huì)調(diào)用odefun成千上萬(wàn)次一次interp1的時(shí)間累積起來(lái)非?？捎^。5.3 參數(shù)敏感性分析與標(biāo)定順序模型參數(shù)太多每次實(shí)驗(yàn)對(duì)不上先去調(diào)哪個(gè)參數(shù)我的經(jīng)驗(yàn)是先調(diào)濃差極化相關(guān)參數(shù)再調(diào)動(dòng)力學(xué)參數(shù)最后調(diào)歐姆參數(shù)。因?yàn)闈獠顦O化參數(shù)對(duì)放電平臺(tái)和容量影響最大方向性最明顯動(dòng)力學(xué)參數(shù)主要影響末端和低溫性能歐姆參數(shù)對(duì)電壓跳變影響大但可以通過(guò)OCV曲線之外的短脈沖實(shí)驗(yàn)單獨(dú)標(biāo)定。參數(shù)標(biāo)定我習(xí)慣用兩步法先用低倍率放電數(shù)據(jù)標(biāo)定熱力學(xué)和擴(kuò)散參數(shù)再用高倍率脈沖數(shù)據(jù)標(biāo)定動(dòng)力學(xué)參數(shù)。低倍率工況下濃差極化小動(dòng)力學(xué)參數(shù)對(duì)端電壓影響有限可以解耦辨識(shí)。反過(guò)來(lái)如果你一開(kāi)始就用4C脈沖標(biāo)定所有參數(shù)會(huì)有嚴(yán)重的參數(shù)不可辨識(shí)問(wèn)題多組參數(shù)都能擬合出接近的結(jié)果但預(yù)測(cè)能力天差地別。6. 模型擴(kuò)展思路從單電池到系統(tǒng)級(jí)應(yīng)用P2D模型跑通之后最順理成章的方向就是往系統(tǒng)級(jí)走。我做過(guò)一次把P2D模型封裝成Simulink S-Function的實(shí)踐把電化學(xué)模型作為被控對(duì)象外面接一個(gè)恒流/恒壓充電控制器再疊加一個(gè)簡(jiǎn)單的溫度模型就能模擬不同充電策略下的內(nèi)部鋰濃度分布和析鋰風(fēng)險(xiǎn)。這個(gè)思路對(duì)快充策略設(shè)計(jì)特別有用——你可以直接在模型里觀察負(fù)極表面鋰濃度是否超過(guò)析鋰閾值而不必等到電池真的析鋰了再去拆解分析。更進(jìn)一步可以做老化機(jī)理模型的耦合。P2D的固相濃度結(jié)果可以接給SEI增長(zhǎng)模型液相濃度和電勢(shì)結(jié)果可以接給析鋰模型這樣仿真就能預(yù)測(cè)循環(huán)壽命和容量衰減趨勢(shì)。雖然目前這類(lèi)耦合模型精度還比較有限但在趨勢(shì)判斷和策略對(duì)比上已經(jīng)很有參考價(jià)值。根據(jù)我個(gè)人的實(shí)操體會(huì)搭建P2D模型這件事最難的不是數(shù)學(xué)推導(dǎo)也不是MATLAB語(yǔ)法而是那些看不見(jiàn)摸不著的物理參數(shù)怎么定、數(shù)值問(wèn)題怎么解、結(jié)果怎么驗(yàn)證。建議新手不要一上來(lái)就追求“完整版P2D熱耦合老化模型”先把最簡(jiǎn)單的恒流放電跑通再把電壓曲線和實(shí)驗(yàn)數(shù)據(jù)對(duì)上一步一步擴(kuò)展。模型每擴(kuò)展一步你對(duì)電池的理解就深一層這才是電化學(xué)仿真最大的價(jià)值所在。本文還有配套的精品資源點(diǎn)擊獲取