JavaScript is not enabled on browser.
更新履歴
Ver.0.13.3 ・spec更新
Ver.0.12.3 ・spec追記
・spec更新
Ver.0.11.3 ・spec追記
Ver.0.10.3 ・時間刻みdtを可変（低温で収束不能の対策）
・spec追記
Ver.0.9.3 ・spec追記
Ver.0.8.3 ・ce時間発展式の致命的な不具合を修正
　事象：ce濃度分布が温度に依らず、ほぼ緩和しない（限界電流密度の温度依存性が不明確）
　対策：fluxの発散項を修正（両辺の次元を統一）
　結果：ce濃度分布の緩和が正常に進行（高温で緩和速度大を確認）
・spec追記
Ver.0.8.2 ・spec追記
Ver.0.7.2 ・spec追記
Ver.0.6.2 ・休止中の電流値を設定（It=0で収束不能の対策）
　It=0　⇒　It_rest=1e-8（極小値）
Ver.0.5.2 ・内部発熱と外部放熱の温度変化を弱連成
Ver.0.4.2 ・spec追記
Ver.0.3.2 ・Rfの温度因子ratio_factorの演算子を訂正（高温で抵抗小）
　*　⇒　/
・Rf2ref/Rf1refの初期値を変更
　1　⇒　1e-1
Ver.0.3.1 ・spec追記
Ver.0.2.1 ・濃度過電圧dphieを液間電位差として追加考慮
・Deの実効値を修正（lee2/lee1を追加）
・収束判定条件を追加
Ver.0.1.1 ・不備の校正
・css追加
Ver.0.1.0 ・draft版

0D-model（反応と拡散を強連成）の特殊な実装例
Ver.0.3.1追記
仮定
Ver.0.2.1更新
Ver.0.11.3追記 		・反応分布が均一（容量型の大電流でmodelling誤差大）
・濃度過電圧を無視dphieを液間電位差として追加考慮（capacitance成分より長時定数で緩和）
前提
Ver.0.10.3追記
Ver.0.12.3更新
Ver.0.12.3追記
Ver.0.13.3更新 		・端子間の充電電流の符号を正
・初回充電||過充放電時に発生する副反応（不可逆な性能低下の要因）は未考慮
・通常反応のcapacitance成分（活物質界面の電気二重層>>圧力変動の圧電効果）を考慮
・ばらつきの均等化（活物質～電極～cell間のSOC等）は未考慮
備考
Ver.0.11.3追記
Ver.0.12.3更新 		・時定数：反応<<capacitance成分<拡散<<ばらつきの均等化
　通常使用範囲内：反応と拡散は弱連成（順次処理の時間発展が一般的）
・単極電位が金属電位より低い場合（psmpe<0）、金属析出の副反応が発生
・析出/溶解の副反応でdendrite成長（拡散律速凝集）の場合、定量化が困難
結果
Ver.0.4.2追記
Ver.0.5.2追記
Ver.0.7.2追記
Ver.0.8.2追記
Ver.0.8.3追記
Ver.0.9.3追記 		・capacitance成分の効果（C大で瞬時の副反応を抑制）：Graphing Calculator-test case-cell0d/-C*10参照
・反応と拡散を強連成の効果（natural）：Graphing Calculator-test case-cell0d-heavy参照
・外部強制の温度変化を弱連成（浮動小数点演算のrobust性を確認）：Graphing Calculator-test case-cell0d-Tout参照
・内部発熱と外部放熱の温度変化を弱連成（強連成の優先度低）：Graphing Calculator-test case-cell0d-Tin/-Trise参照
　電流条件：通常使用範囲外の超大電流（10秒以内に塩枯れ後、未考慮の副反応が優勢かつ反応分布が極端化）
　内部発熱：電圧変化量Vt-OCVtから非厳密の発熱量を仮定
　外部放熱：熱伝達係数alphaを調整（熱伝導を含めて換算）
・まとめ（反応と拡散を強連成＋温度変化を弱連成）：Graphing Calculator-test case-cell0d-Ta=-30~60参照
・濃度変化（修正後、高温で緩和速度大を確認）：Graphing Calculator-test case-cell0d-ce2/cse2参照
・ce濃度変化なしdphie≒0の全固体想定（低温特性の改善を確認）：Graphing Calculator-test case-cell0d-AS参照
　模擬条件：Deref*1e8（陰解法）

unknowns
j= 2: negative electrode s: separator 1: positive electrode
jLi jLi2R jLi1R
eta eta2 eta1
ce ce2 ces ce1
cse cse2 cse1
T
Ver.0.5.2追加 		T

constants
item unit define
TK K $TK=273.15;
Rgas J/(mol*K) $Rgas=8.31446261815324;
F C/mol $F=96485.33212331;

input parameters
item unit define
valence
Ver.0.1.1統合 		n=1;
anodic transfer coefficient aa2=0.5;
aa1=0.5;
transference number of ion t0=0.4;
component thickness cm L2=1e-3;
Ls=1e-3;
L1=1e-3;
thickness
Ver.0.5.2追加 		cm L=L2+Ls+L1;
radius of active material cm rs2=L2/2;
rs1=L1/2;
volume fraction of electrolyte epse2=0.4;
epses=0.7;
epse1=0.4;
volume fraction of active material epss2=(4*pi/3)/2**3;
epss1=(4*pi/3)/2**3;
diffusion length of electrolyte
Ver.0.2.1追加 		cm lee2=0.5*(epse2**1.5*L2+epses**1.5*Ls)/((epse2*epses)**1.5);
lee1=0.5*(epses**1.5*Ls+epse1**1.5*L1)/((epses*epse1)**1.5);
diffusion length of active material cm lse2=rs2/5;
lse1=rs1/5;
specific surface area cm^2/cm^3 as2=3*epss2/rs2;
as1=3*epss1/rs1;
volume averaged concentration of electrolyte mol/cm^3 ce0=1e-3;
volume averaged maximum concentration of active material mol/cm^3 csmax2=1e-2;
csmax1=1e-2;
volumetric heat capacity
Ver.0.5.2追加 		J/(cm^3*K) rhocp=1;
heat transfer coefficient
Ver.0.5.2追加 		W/(cm^2*K) alpha=1e-3;
reference temperature K Tref=TK+15;
reference values
Ver.0.1.1変更
Ver.0.3.2変更 		F/cm^2
F/cm^2
Ohm*cm^2
A/cm^2
A/cm^2
Ohm*cm^2
Ohm*cm^2
cm^2/s
cm^2/s
cm^2/s 		C2=1;
C1=1;
R0ref=1;
i02ref=1e-3;
i01ref=1e-3;
Rf2ref=1e-1;
Rf1ref=1e-1;
Deref=1e-6;
Ds2ref=1e-8;
Ds1ref=1e-8;
activation energies
Ver.0.1.1訂正 		J/mol ER0=20e3;
Ei02=50e3;
Ei01=50e3;
ERf2=50e3;
ERf1=50e3;
EDe=30e3;
EDs2=10e3;
EDs1=10e3;

input initial conditions
item unit define
time step
Ver.0.1.1追加
Ver.0.10.3追加 		s dt=1e-2;
get_dt(t)=t<=1e-2&&&1e-3|||1e-2;
time
Ver.0.1.1追加 		s t=0;
temperature K T=Tref;
ambient temperature
Ver.0.5.2追加 		K Ta=Tref;
current
Ver.0.6.2追加 		A/cm^2 It=1e-3;
It_rest=1e-8;
previous values
Ver.0.1.1変更
Ver.0.5.2追加 		mol/cm^3
mol/cm^3
mol/cm^3
mol/cm^3
V
V
K 		ce20=ce0;
ce10=ce0;
cs20=csmax2*1e-3;
cs10=csmax1*0.999;
psmpe20=U2(cs20/csmax2)=>;
psmpe10=U1(cs10/csmax1)=>;
T0=Ta;
unknowns
Ver.0.5.2追加 		A/cm^3
A/cm^3
V
V
mol/cm^3
mol/cm^3
mol/cm^3
mol/cm^3
K 		jLi2R=0;
jLi1R=0;
eta2=0;
eta1=0;
ce2=ce20;
ce1=ce10;
cse2=cs20;
cse1=cs10;
T=T0;

equations
item define
dphie
Ver.0.2.1仮定 		dphie=ln(ce1/ce2)*2*Rgas*T*(1-t0)/(n*F)
output voltage
Ver.0.2.1更新 		Vt=psmpe1-psmpe2+R0*It+dphie
convergence conditions of absolute residual norm
Ver.0.2.1追加 		isNaN(norm)||norm<1e-10
break conditions
Ver.0.1.1追加 		le(ce2,0)||le(ce1,0)||le(cse2,0)||ge(cse2,csmax2)||le(cse1,0)||ge(cse1,csmax1)
time evolution
Ver.0.1.1追加 		t=t+dt
temperature factor ratio_factor(T,Eact)=<exp((1/Tref-1/T)*Eact/Rgas);
Ohmic resistance R0=R0ref*(0.5+0.5/ratio_factor(T,ER0)=>)
exchange current density
Ver.0.1.1訂正 		i02=i02ref*ratio_factor(T,Ei02)=>*(ce2/ce0)**aa2*(1-cse2/csmax2)**aa2*(cse2/csmax2)**(1-aa2)
i01=i01ref*ratio_factor(T,Ei01)=>*(ce1/ce0)**aa1*(1-cse1/csmax1)**aa1*(cse1/csmax1)**(1-aa1)
SEI film resistance
Ver.0.3.2訂正 		Rf2=Rf2ref/ratio_factor(T,ERf2)=>
Rf1=Rf1ref/ratio_factor(T,ERf1)=>
diffusion coefficient De De=Deref*ratio_factor(T,EDe)=>
diffusion coefficient Ds Ds2=Ds2ref*ratio_factor(T,EDs2)=>
Ds1=Ds1ref*ratio_factor(T,EDs1)=>
Open Circuit Potential f(t,a,b,c)=<b+a*exp(c*t);
g(t,a1,b1,c1,a2,b2,c2)=<f(1-t,a1,b1,c1)=>+f(t,a2,b2,c2)=>;
U2(t)=<g(t,-0.1,0.1,-100,2.5,0.1*(1-t),-30)=>;
U1(t)=<g(t,-(4.5-1),4.5-t,-50,0.3,-0.3*(1-t),-20)=>;
potential psmpe2=U2(cse2/csmax2)=>+eta2+Rf2*jLi2R/as2
psmpe1=U1(cse1/csmax1)=>+eta1+Rf1*jLi1R/as1
double-layer capacitance
Ver.0.3.1削除 		jLi2C=C2*(psmpe2-psmpe20)/dt
jLi1C=C1*(psmpe1-psmpe10)/dt
volume averaged current density jLi2=-It/L2
jLi1=It/L1
charge conservation
Ver.0.3.2訂正
Ver.0.8.3追加 		jLi2=jLi2R+jLi2C
jLi1=jLi1R+jLi1C
-It/L2=jLi2R+C2*(U2(cse2/csmax2)=>+eta2+(Rf2ref/ratio_factor(T,ERf2)=>)*jLi2R/as2-psmpe20)/dt;
It/L1=jLi1R+C1*(U1(cse1/csmax1)=>+eta1+(Rf1ref/ratio_factor(T,ERf1)=>)*jLi1R/as1-psmpe10)/dt;
Butler-Volmer equation
Ver.0.1.1訂正
Ver.0.12.3追記 		jLi2R/as2=i02*(exp(eta2*aa2*n*F/(Rgas*T))-exp(-eta2*(1-aa2)*n*F/(Rgas*T)))
jLi1R/as1=i01*(exp(eta1*aa1*n*F/(Rgas*T))-exp(-eta1*(1-aa1)*n*F/(Rgas*T)))
jLi2R/as2=i02ref*ratio_factor(T,Ei02)=>*(ce2/ce0)**aa2*(1-cse2/csmax2)**aa2*(cse2/csmax2)**(1-aa2)*(exp(eta2*aa2*n*F/(Rgas*T))-exp(-eta2*(1-aa2)*n*F/(Rgas*T)));
jLi1R/as1=i01ref*ratio_factor(T,Ei01)=>*(ce1/ce0)**aa1*(1-cse1/csmax1)**aa1*(cse1/csmax1)**(1-aa1)*(exp(eta1*aa1*n*F/(Rgas*T))-exp(-eta1*(1-aa1)*n*F/(Rgas*T)));
average concentration cet=ce0*(epse2*L2+epses*Ls+epse1*L1)
ces ces=(cet-(ce2*epse2*L2+ce1*epse1*L1))/(epses*Ls)
ce
Ver.0.1.1訂正
Ver.0.2.1更新
Ver.0.8.3修正 		epse2*(ce2-ce20)/dt=(De/lee2)*(ces-ce2)/L2+jLi2*(1-t0)/(n*F)
epse1*(ce1-ce10)/dt=(De/lee1)*(ces-ce1)/L1+jLi1*(1-t0)/(n*F)
epse2*(ce2-ce20)/dt=(Deref*ratio_factor(T,EDe)=>/lee2)*((ce0*(epse2*L2+epses*Ls+epse1*L1)-(ce2*epse2*L2+ce1*epse1*L1))/(epses*Ls)-ce2)/L2-(It/L2)*(1-t0)/(n*F);
epse1*(ce1-ce10)/dt=(Deref*ratio_factor(T,EDe)=>/lee1)*((ce0*(epse2*L2+epses*Ls+epse1*L1)-(ce2*epse2*L2+ce1*epse1*L1))/(epses*Ls)-ce1)/L1+(It/L1)*(1-t0)/(n*F);
cs
Ver.0.1.1訂正 		cs2=cs20-dt*jLi2R/(epss2*n*F)
cs1=cs10-dt*jLi1R/(epss1*n*F)
cse
Ver.0.1.1訂正 		(Ds2/lse2)*(cse2-cs2)=-jLi2R/(as2*n*F)
(Ds1/lse1)*(cse1-cs1)=-jLi1R/(as1*n*F)
(Ds2ref*ratio_factor(T,EDs2)=>/lse2)*(cse2-(cs20-dt*jLi2R/(epss2*n*F)))=-jLi2R/(as2*n*F);
(Ds1ref*ratio_factor(T,EDs1)=>/lse1)*(cse1-(cs10-dt*jLi1R/(epss1*n*F)))=-jLi1R/(as1*n*F);
Open Circuit Voltage
Ver.0.5.2追加 		OCVt=U1(cs1/csmax1)=>-U2(cs2/csmax2)=>
T
Ver.0.5.2仮定 		L*rhocp*(T-T0)/dt=It*(Vt-OCVt)-alpha*(T-Ta)
T=(L*rhocp*T0+(It*(Vt-OCVt)+alpha*Ta)*dt)/(L*rhocp+alpha*dt);