% Program_Level = Canner*_Single_Cell_Discharge_&_Charge_with_Constant_Current ;
程式等級﹕製罐級*，只有定電流功能，用來估算單顆電池之充電和放電行為；
% Default_Battery_Type = Panasonic_NCR18650_3100mAh_NCA_LixC6 ;
預設電池型號﹕松下Panasonic_NCR18650_3100mAh，鎳鈷鋁氧化物/鋰碳化合物；
% Author = H5N1AIDS , National_Taiwan_University ^_^ ; 程式撰寫人﹕H5N1AIDS，國立台灣大學 ^_^；
% E-mail = [email protected] ; 電子郵件﹕[email protected]；
% Last_Update_Date = 2016 / 07 / 10 ; 最後更新日期﹕2016年7月10日；
% Main_References = Journal_of_Power_Sources_252_(2014)_214-228 && Journal_of_Power_Sources_247_(2014)_365-376； 主要參考文獻；
% @ Copyright 請勿任意翻印
% $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ clc 清空執行視窗
clear 清除暫存工作空間之變數
% ________________________________________________
C_rate = 1 ; 設定充放電速率；
Lower_Cut_off_Voltage = 3 ; 最低截止電壓；
Upper_Cut_off_Voltage = 4.3 ; 最高截止電壓；
% ________________________________________________
T_constant = 1 ; 定溫功能，1表示啟動，0代表關閉；
T_ambient = 298.15 ; 環境溫度；
Convection_coefficient = 10 ; 對流係數；
doi:10.6342/NTU201602734
129
% ________________________________________________
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%基本參數設定開始 Standard_Ah = 3.1 ; 標準安時數；
C_solid_n_max = 28000 ; 負電極最大之鋰原子濃度；
C_solid_p_max = 51000 ; 正電極最大之鋰原子濃度；
Initial_Available_ratio_n = 0.63 ; 負電極最大之鋰原子濃度，其放電初始之可使用比例；
Initial_Available_ratio_p = 0.395 ; 正電極最大之鋰原子濃度，其放電初始之可使用比例；
C_electrolyte_n_avg = 1200 ; 負極區電解液濃度；
C_electrolyte_s_avg = 1200 ; 隔離膜區電解液濃度；
C_electrolyte_p_avg = 1200 ; 正極區電解液濃度；
Activation_energy_D_Solid_n = 5000 ; 負電極內鋰原子擴散之活化能；
Activation_energy_D_Solid_p = 25000 ; 正電極內鋰原子擴散之活化能；
Activation_energy_kn_reaction = 20000 ; 負電極之電化學反應活化能；
Activation_energy_kp_reaction = 35000 ; 正電極之電化學反應活化能；
Electrode_area_n = 0.132 ; 負電極之投影面積；
Electrode_area_p = 0.132 ; 正電極之投影面積；
Rn = 12.5*10^-6 ; 負電極顆粒半徑；
Rp = 4*10^-9 ; 正電極顆粒半徑；
Ln_space = 93*10^-6 ; 負極空間整體長度；
Ls_space = 17*10^-6 ; 隔離膜空間整體長度；
Lp_space = 69*10^-6 ; 正極空間整體長度；
Ln_electrolyte = 93*10^-6 ; 負極區之電解液空間長度；
Ls_electrolyte = 17*10^-6 ; 隔離膜區之電解液空間長度；
Lp_electrolyte = 69*10^-6 ; 正極區之電解液空間長度；
Ln_solid = 93*10^-6 ; 負電極長度；
Lp_solid = 69*10^-6 ; 正電極長度；
doi:10.6342/NTU201602734
130
dln_activity_coefficient_n_over_dln_C_electrolyte_n_avg = 1 ; = 1；
dln_activity_coefficient_s_over_dln_C_electrolyte_s_avg = 1 ; = 1；
dln_activity_coefficient_p_over_dln_C_electrolyte_p_avg = 1 ; = 1；
tn_plus = 0.38 ; 負極區電解液之轉移係數；
ts_plus = 0.38 ; 隔離膜區電解液之轉移係數；
tp_plus = 0.38 ; 正極區電解液之轉移係數；
D_solid_p_0 = NaN ; 正電極內的鋰原子擴散係數= Not a Number。由於帶入隨時間演變之經驗公式，故不在此處給值設定；
D_solid_n_0 = 3.9.*10.^-14 ; 負電極內鋰原子的擴散係數；
D_electrolyte_p_0 = 2.6.*10.^-10 ; 正極區電解液之擴散係數；
D_electrolyte_s_0 = 2.6.*10.^-10 ; 隔離膜區電解液之擴散係數；
D_electrolyte_n_0 = 2.6.*10.^-10 ; 負極區電解液之擴散係數；
D_tuning_parameter = 4/3 ; 電極內鋰原子擴散通道係數；
kp_0 = 1*10^-10 ; 正極區之電化學反應速率常數；
kn_0 = 2.334*10^-11 ; 負極區之電化學反應速率常數；
Solid_epsilon_n = 0.5802 ; 負電極之體積分率；
Solid_epsilon_p = 0.4483 ; 正電極之體積分率；
Solid_epsilon_n_0 = 0.5802 ; 負電極之體積分率初始值；
Solid_epsilon_p_0 = 0.4483 ; 正電極之體積分率初始值；
Electorlyte_epsilon_n = 0.364 ; 負極區電解液之體積分率；
Electorlyte_epsilon_s = 0.5 ; 隔離膜區電解液之體積分率；
Electorlyte_epsilon_p = 0.3382 ; 正極區電解液之體積分率；
Filler_Volume_fraction_n = 1 - Solid_epsilon_n - Electorlyte_epsilon_n ; 負極區粘著物之體積分率；
Filler_Volume_fraction_p = 1 - Solid_epsilon_p - Electorlyte_epsilon_p ; 正極區粘著物之體積分率；
Bruggman_coefficient_n = 1.5 ; 負極區的Bruggman指數；
Bruggman_coefficient_s = 1.5 ; 隔離膜區的Bruggman指數；
doi:10.6342/NTU201602734
131
Bruggman_coefficient_p = 1.5 ; 正極區的Bruggman指數；
sigma_n = 100 ; 負電極單位長度電導；
sigma_p = 100 ; 正電極單位長度電導；
Length_collector_n_Cu = 25.*10.^-6 ; 負極集電器長度；
Length_collector_p_Al = 25.*10.^-6 ; 正極集電器長度；
Conductivity_collector_n_Cu = 5.998.*10.^7 ; 負極集電器之單位長度電導率；
Conductivity_collector_p_Al = 3.774.*10.^7 ; 正極集電器之單位長度電導率；
alpha = 0.5 ; 傳遞係數，正、負極同值；
IR_contact_n = 0 ; 負極之單位面積接觸電阻值；
IR_contact_p = 0 ; 正極之單位面積接觸電阻值；
R_SEI = 0.04 ; SEI膜之單位面積電阻值；
Cp = 750 ; 電池之定壓熱容值；
Battery_density = 2017 ; 電池密度；
Battey_volume = pi.*9.^2.*65e-9 ; 電池體積；
Battery_surface_area = 2.*pi.*9e-3.*65e-3 + 2.*pi.*9e-3.*9e-3 ; 電池表面積；
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
t = 0 ; 初始時序；
dt = 1 ; 時間步階長度；
cyc = 1 ; 初始循環序號；
Power_output = 0 ; 初始之電池輸出功率值；
Phi_solid_n_terminal = 0 ; 定義負電極為零電位面；
current = Standard_Ah.*C_rate ; 定義電流，放電時為正值，充電時為負值；
F = 96487 ; 法拉第常數；
R = 8.314 ; 理想氣體常數；
Delta_n = Ln_electrolyte./3 ; 負極區電解液之微量擴散長度；
Delta_p = Lp_electrolyte./3 ; 正極區電解液之微量擴散長度；
doi:10.6342/NTU201602734
132
Diffusion_length_s = Ls_electrolyte ; 隔離膜區電解液之微量擴散長度；
N_n = (Solid_epsilon_n.*Electrode_area_n.*Ln_solid)./(4./3.*pi.*Rn.^3) ; 負電極顆粒總數；
N_p = (Solid_epsilon_p.*Electrode_area_p.*Lp_solid)./(4./3.*pi.*Rp.^3) ; 正電極顆粒總數；
T_reference = 298.15 ; 電化學參數之參考溫度標準；
T_core = T_ambient ; 初始電池溫度等於環境溫度；
R_collector_n_Cu = Length_collector_n_Cu/Conductivity_collector_n_Cu ; 負極集電器之單位面積電阻值；
R_collector_p_Al = Length_collector_p_Al/Conductivity_collector_p_Al ; 正極集電器之單位面積電阻值；
Modify_sigma_n = sigma_n.*Solid_epsilon_n ; 負電極單位長度電導值之多孔性修正；
Modify_sigma_p = sigma_p.*Solid_epsilon_p ; 正電極單位長度電導值之多孔性修正；
as_n = (3.*Solid_epsilon_n)./Rn ; 負電極每單位體積之有效表面積；
as_p = (3.*Solid_epsilon_p)./Rp ; 正電極每單位體積之有效表面積；
aLp = as_p.*Lp_solid ; 正電極每單位面積之有效表面積；
aLn = as_n.*Ln_solid ; 負電極每單位面積之有效表面積；
tp_plus_bar = tp_plus - 1 ; 正極區電解液之轉移係數減一；
tn_plus_bar = tn_plus - 1 ; 負極區電解液之轉移係數減一；
ts_plus_bar = ts_plus - 1 ; 隔離膜區電解液之轉移係數減一；
C_solid_n_surface = C_solid_n_max*Initial_Available_ratio_n ; 初始負電極表面之鋰原子濃度；
C_solid_p_surface = C_solid_p_max*Initial_Available_ratio_p ; 初始正電極表面之鋰原子濃度；
U_SOC_p = C_solid_p_surface/C_solid_p_max ; 初始正電極鋰原子濃度之儲量比值；
U_SOC_n = C_solid_n_surface/C_solid_n_max ; 初始負電極鋰原子濃度之儲量比值；
dUp_over_dT = -2743.152.*U_SOC_p.^10 + 16813.562.*U_SOC_p.^9 - 45715.78.*U_SOC_p.^8 + 72552.932.*U_SOC_p.^7 -
74367.873.*U_SOC_p.^6 + 51403.626.*U_SOC_p.^5 - 24247.543.*U_SOC_p.^4 + 7702.849.*U_SOC_p.^3 - 1576.405.*U_SOC_p.^2 + 187.617.*U_SOC_p - 9.859 ; 正電極開路電壓隨溫度變化；
doi:10.6342/NTU201602734
133
dUn_over_dT = (0.00527 + 3.29927.*U_SOC_n - 91.79326.*U_SOC_n.^2 + 1004.91101.*U_SOC_n.^3 - 5812.27813.*U_SOC_n.^4 + 19329.7549.*U_SOC_n.^5 - 37147.8947.*U_SOC_n.^6 + 38379.18127.*U_SOC_n.^7 - 16515.05308.*U_SOC_n.^8)./(1 -
48.09287.*U_SOC_n + 1017.2348.*U_SOC_n.^2 - 10481.80419.*U_SOC_n.^3 + 59431.30001.*U_SOC_n.^4 -
195881.6488.*U_SOC_n.^5 + 374577.3152.*U_SOC_n.^6 - 385821.1607.*U_SOC_n.^7 + 165705.8597.*U_SOC_n.^8) ; 負電極開路電壓隨溫度變化；
if 0.36 <= U_SOC_p && U_SOC_p <= 0.41 正電極開路電壓隨U_SOC_p變化；
Up = 8.535 - 17.059.*U_SOC_p + 21.038.*U_SOC_p.^2 - 9.153.*U_SOC_p.^3 + 9.875.*(U_SOC_p - 0.7).^3 - 2.176.*
(U_SOC_p - 0.55).^3 - 1331.866.*(U_SOC_p - 0.41).^3 ; elseif 0.41 <= U_SOC_p && U_SOC_p <= 0.55
Up = 8.535 - 17.059.*U_SOC_p + 21.038.*U_SOC_p.^2 - 9.153.*U_SOC_p.^3 + 9.875.*(U_SOC_p - 0.7).^3 - 2.176.*
(U_SOC_p - 0.55).^3 ;
elseif 0.55 <= U_SOC_p && U_SOC_p <= 0.7
Up = 8.535 - 17.059.*U_SOC_p + 21.038.*U_SOC_p.^2 - 9.153.*U_SOC_p.^3 + 9.875.*(U_SOC_p - 0.7).^3 ; elseif 0.7 <= U_SOC_p && U_SOC_p <= 0.935
Up = 8.535 - 17.059.*U_SOC_p + 21.038.*U_SOC_p.^2 - 9.153.*U_SOC_p.^3 + 9.875.*(U_SOC_p - 0.7).^3;
elseif 0.935 <= U_SOC_p && U_SOC_p <= 0.959
Up = 8.535 - 17.059.*U_SOC_p + 21.038.*U_SOC_p.^2 - 9.153.*U_SOC_p.^3 + 9.875.*(U_SOC_p - 0.7).^3 - 5370.872.*(U_SOC_p - 0.935).^3 ;
elseif 0.959 <= U_SOC_p && U_SOC_p <= 0.98
Up = 8.535 - 17.059.*U_SOC_p + 21.038.*U_SOC_p.^2 - 9.153.*U_SOC_p.^3 + 9.875.*(U_SOC_p - 0.7).^3 - 5370.872.*(U_SOC_p - 0.935).^3 - 47690.304.*(U_SOC_p - 0.959).^3 ;
end
Up = Up + 0.001.*(T_core - T_reference).*dUp_over_dT ; 正電極開路電壓隨U_SOC_p與溫度變化；
Un = 0.6379 + 0.5416.*exp(-305.5309.*U_SOC_n) + 0.044.*tanh((0.1958 - U_SOC_n)./0.1088) - 0.1978.*tanh((U_SOC_n - 1.0571)./0.0854) - 0.6875.*tanh((U_SOC_n + 0.0117)./0.0529) - 0.0175.*tanh((U_SOC_n - 0.5692)./0.0875) + 0.001.*(T_core - T_reference).*dUn_over_dT ; 負電極開路電壓隨U_SOC_n與溫度變化；
doi:10.6342/NTU201602734
134
Battery_Working_Potential = Up - Un ; 電池之開路電壓；
Battery_density_dot_Cp = Battery_density.*Cp ; 電池密度與定壓熱容之乘積；
Eigenvalues_square =
[4.4934094579091;7.7252518369377;10.9041216594289;14.0661939128315;17.2207552719308;20.3713029592876;23.51945249 86890;26.6660542588127;29.8115987908930;36.1006222443756;39.2444323611642;42.3879135681319;45.5311340139913;48.6
741442319544;51.8169824872797].^2 ; 球座標擴散之級數解其前15項特徵值的平方；
Ei_length = length(Eigenvalues_square) ; 特徵值矩陣其行向量之長度；
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x1 = zeros(1,3700/C_rate) ; 資料保存空間儲位；
x2 = zeros(1,3700/C_rate) ; 資料保存空間儲位；
% x3 = zeros(1,3700/C_rate) ; 資料保存空間儲位；
% x4 = zeros(1,3700/C_rate) ; 資料保存空間儲位；
% x5 = zeros(1,3700/C_rate) ; 資料保存空間儲位；
% x6 = zeros(1,3700/C_rate) ; 資料保存空間儲位；
%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%基本參數設定結束 tic 運算計時開始
while Battery_Working_Potential >= Lower_Cut_off_Voltage
只要電池工作電壓仍然大於或等於最低截止電壓，運算就繼續執行 t = t + dt ; 時間演進；
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
iapp_p = current./Electrode_area_p ; 正極工作電流密度；
iapp_n = current./Electrode_area_n ; 負極工作電流密度；
Jp = iapp_p./aLp ; 正電極之Butler-Volmer電流；
Jn = -iapp_n./aLn ; 負電極之Butler-Volmer電流；
IR_SEI = -Jn.*R_SEI ; SEI膜造成之電位損失；
RT = R.*T_core ; 為減少重複運算而設定之參數；
doi:10.6342/NTU201602734
135
R_diff_T = (1./T_reference - 1./T_core)./R ; 為減少重複運算而設定之參數；
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %% % D_solid_p =3.*10.^-15.*((1 + tanh(-20.*(U_SOC_p - 0.73))) + 0.002).*exp(Activation_energy_D_Solid_p.*R_diff_T) ; 正電極內鋰原子之擴散係數；
D_solid_n = D_solid_n_0.*exp(Activation_energy_D_Solid_n.*R_diff_T) ; 負電極內鋰原子之擴散係數；
D_electrolyte_p = (10.^(-8.443-54./(T_core-229-0.005.*C_electrolyte_p_avg)-2.2.*10.^
-4.*C_electrolyte_p_avg)).*(Electorlyte_epsilon_p.^Bruggman_coefficient_p);
正極區內鋰離子之擴散係數；
D_electrolyte_s = (10.^(-8.443-54./(T_core-229-0.005.*C_electrolyte_s_avg)-2.2.*10.^
-4.*C_electrolyte_s_avg)).*(Electorlyte_epsilon_s.^Bruggman_coefficient_s);
隔離膜區內鋰離子之擴散係數；
D_electrolyte_n = (10.^(-8.443-54./(T_core-229-0.005.*C_electrolyte_n_avg)-2.2.*10.^
-4.*C_electrolyte_n_avg)).*(Electorlyte_epsilon_n.^Bruggman_coefficient_n) ; 負極區內鋰離子之擴散係數；
Diffusion_length_p = Delta_p.*(1 - exp(-D_tuning_parameter.*sqrt(D_electrolyte_p.*t)./Delta_p)) ; 正極區內鋰離子之擴散長度；
Diffusion_length_n = Delta_n.*(1 - exp(-D_tuning_parameter.*sqrt(D_electrolyte_p.*t)./Delta_n)) ; 負極區內鋰離子之擴散長度；
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % Chicken_p = 0 ; 正電極內鋰原子濃度擴散級數解之儲位設定與歸零；
Chicken_n = 0 ; 負電極內鋰原子濃度擴散級數解之儲位設定與歸零；
其中因為第一次的設定代碼為G，因此事後整理時，將G改為Chicken。
for i = 1 : Ei_length 級數解演算疊加開始；
Dinosaur_p = 2.*(exp(-Eigenvalues_square(i).*D_solid_p.*t./Rp.^2)) ; 又因為雞走路跟暴龍很像，所以又令一代碼為Dinosaur。
Chicken_p = Chicken_p + Dinosaur_p ;
doi:10.6342/NTU201602734
136
Dinosaur_n = 2.*(exp(-Eigenvalues_square(i).*D_solid_n.*t./Rn.^2)) ; Chicken_n = Chicken_n + Dinosaur_n ;
end 級數解演算疊加結束；
U_SOC_p = U_SOC_p + Jp.*dt./(Rp.*C_solid_p_max.*F).*(3 + Chicken_p) ; 正電極內鋰原子濃度儲藏比值之變化；
U_SOC_n = U_SOC_n + Jn.*dt./(Rn.*C_solid_n_max.*F).*(3 + Chicken_n) ; 負電極內鋰原子濃度儲藏比值之變化；
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % kp = kp_0.*exp(Activation_energy_kp_reaction.*R_diff_T) ;
正極電化學反應速率常數之溫度修正；
kn = kn_0.*exp(Activation_energy_kn_reaction.*R_diff_T) ; 負極電化學反應速率常數之溫度修正；
kappa_p = (-8.2488 + 0.053248.*T_core - 0.000029871.*T_core.^2 + 0.26235.*C_electrolyte_p_avg.*10^-3 - 0.0093063.*C_electrolyte_p_avg.*10^-3.*T_core + 0.000008069.*C_electrolyte_p_avg.*10^-3.*T_core.^2 + 0.22002.*(C_electrolyte_p_avg.*10^-3).^2 - 0.0001765.*
(C_electrolyte_p_avg.*10^-3).^2.*T_core ).^2.*C_electrolyte_p_avg.*10^-3.*0.1 ; 正極區內電解液之導電度；
kappa_s = (-8.2488 + 0.053248.*T_core - 0.000029871.*T_core.^2 + 0.26235.*C_electrolyte_s_avg.*10^-3 - 0.0093063.*C_electrolyte_s_avg.*10^-3.*T_core + 0.000008069.*C_electrolyte_s_avg.*10^-3.*T_core.^2 + 0.22002.*(C_electrolyte_s_avg.*10^-3).^2 - 0.0001765.*
(C_electrolyte_s_avg.*10^-3).^2.*T_core ).^2.*C_electrolyte_s_avg.*10^-3.*0.1 ; 隔離膜區內電解液之導電度；
kappa_n = (-8.2488 + 0.053248.*T_core - 0.000029871.*T_core.^2 + 0.26235.*C_electrolyte_n_avg.*10^-3 - 0.0093063.*C_electrolyte_n_avg.*10^-3.*T_core + 0.000008069.*C_electrolyte_n_avg.*10^-3.*T_core.^2 + 0.22002.*(C_electrolyte_n_avg.*10^-3).^2 - 0.0001765.*
(C_electrolyte_n_avg.*10^-3).^2.*T_core ).^2.*C_electrolyte_n_avg.*10^-3.*0.1 ; 負極區內電解液之導電度；
doi:10.6342/NTU201602734
137
dUp_over_dT = -2743.152.*U_SOC_p.^10 + 16813.562.*U_SOC_p.^9 - 45715.78.*U_SOC_p.^8 + 72552.932.*U_SOC_p.^7 - 74367.873.*U_SOC_p.^6 + 51403.626.*U_SOC_p.^5 - 24247.543.*U_SOC_p.^4 + 7702.849.*U_SOC_p.^3
- 1576.405.*U_SOC_p.^2 + 187.617.*U_SOC_p - 9.859 ; 正電極開路電壓隨溫度變化；
dUn_over_dT = (0.00527 + 3.29927.*U_SOC_n - 91.79326.*U_SOC_n.^2 + 1004.91101.*U_SOC_n.^3
- 5812.27813.*U_SOC_n.^4 + 19329.7549.*U_SOC_n.^5 - 37147.8947.*U_SOC_n.^6 + 38379.18127.*U_SOC_n.^7 - 16515.05308.*U_SOC_n.^8)./(1 - 48.09287.*U_SOC_n + 1017.2348.*U_SOC_n.^2 - 10481.80419.*U_SOC_n.^3 + 59431.30001.*U_SOC_n.^4 - 195881.6488.*U_SOC_n.^5 + 374577.3152.*U_SOC_n.^6 - 385821.1607.*U_SOC_n.^7 + 165705.8597.*U_SOC_n.^8) ; 負電極開路電壓隨溫度變化；
if 0.36 <= U_SOC_p && U_SOC_p <= 0.41 正電極開路電壓隨U_SOC_p變化；
Up = 8.535 - 17.059.*U_SOC_p + 21.038.*U_SOC_p.^2 - 9.153.*U_SOC_p.^3 + 9.875.*(U_SOC_p - 0.7).^3 - 2.176.*(U_SOC_p - 0.55).^3 - 1331.866.*(U_SOC_p - 0.41).^3 ;
elseif 0.41 <= U_SOC_p && U_SOC_p <= 0.55
Up = 8.535 - 17.059.*U_SOC_p + 21.038.*U_SOC_p.^2 - 9.153.*U_SOC_p.^3 + 9.875.*(U_SOC_p - 0.7).^3 - 2.176.*(U_SOC_p - 0.55).^3 ;
elseif 0.55 <= U_SOC_p && U_SOC_p <= 0.7
Up = 8.535 - 17.059.*U_SOC_p + 21.038.*U_SOC_p.^2 - 9.153.*U_SOC_p.^3 + 9.875.*(U_SOC_p - 0.7).^3 ; elseif 0.7 <= U_SOC_p && U_SOC_p <= 0.935
Up = 8.535 - 17.059.*U_SOC_p + 21.038.*U_SOC_p.^2 - 9.153.*U_SOC_p.^3 + 9.875.*(U_SOC_p - 0.7).^3;
elseif 0.935 <= U_SOC_p && U_SOC_p <= 0.959
Up = 8.535 - 17.059.*U_SOC_p + 21.038.*U_SOC_p.^2 - 9.153.*U_SOC_p.^3 + 9.875.*(U_SOC_p - 0.7).^3 - 5370.872.*(U_SOC_p - 0.935).^3 ;
elseif 0.959 <= U_SOC_p && U_SOC_p <= 0.98
Up = 8.535 - 17.059.*U_SOC_p + 21.038.*U_SOC_p.^2 - 9.153.*U_SOC_p.^3 + 9.875.*(U_SOC_p - 0.7).^3 - 5370.872.*(U_SOC_p - 0.935).^3 - 47690.304.*(U_SOC_p - 0.959).^3 ;
end
Up = Up + 0.001.*(T_core - T_reference).*dUp_over_dT ; 正電極開路電壓隨U_SOC_p與溫度變化；
doi:10.6342/NTU201602734
138
Un = 0.6379 + 0.5416.*exp(-305.5309.*U_SOC_n) + 0.044.*tanh((0.1958 - U_SOC_n)./0.1088) - 0.1978.*
tanh((U_SOC_n - 1.0571)./0.0854) - 0.6875.*tanh((U_SOC_n + 0.0117)./0.0529) - 0.0175.*
tanh((U_SOC_n - 0.5692)./0.0875) + 0.001.*(T_core - T_reference).*dUn_over_dT ; 負電極開路電壓隨U_SOC_p與溫度變化；
% % % % % % % % % % % % % % % % % % % % % % %電池電位演算開始，詳見論文第三章。
Alpha_n = D_electrolyte_n./Diffusion_length_n ; Beta_sep = D_electrolyte_s./Diffusion_length_s ; Gamma_p = D_electrolyte_p./Diffusion_length_p ; Omega = Gamma_p.*Alpha_n./Beta_sep ; GAO = Gamma_p + Alpha_n + Omega ; GA_over_GAO = Gamma_p.*Alpha_n./GAO ; An = -GA_over_GAO./Ln_electrolyte ; Bn = -An ;
Ap = GA_over_GAO./Lp_electrolyte ; Bp = -Ap ;
Ce_interface1 = ((Alpha_n + Omega).*C_electrolyte_n_avg+ Gamma_p.*C_electrolyte_p_avg)./ GAO ; Ce_interface2 = (Alpha_n.*C_electrolyte_n_avg + (Gamma_p + Omega).*C_electrolyte_p_avg)./ GAO ; Phi_solid_n_avg = -iapp_n.*Ln_solid./(3.*Modify_sigma_n) ;
Phi_electrolyte_n_avg = Phi_solid_n_avg - Un - IR_SEI - RT./(alpha.*F).*asinh(-Jn.*(2.*F.*kn).^-1.*
((C_solid_n_max -C_solid_n_surface).*C_solid_n_surface.*C_electrolyte_n_avg).^-alpha) ;
Kappa_eff_n = kappa_n.*RT.*tn_plus_bar.*(1+ dln_activity_coefficient_n_over_dln_C_electrolyte_n_avg)./F.^2 ; S_interface1 = iapp_n + 3.*Kappa_eff_n.*(1 - C_electrolyte_n_avg./Ce_interface1)./Ln_electrolyte ;
Phi_interface1 = S_interface1.*Ln_electrolyte./3 + Phi_electrolyte_n_avg ;
Kappa_eff_s = kappa_s.*RT.*ts_plus_bar.*(1+ dln_activity_coefficient_s_over_dln_C_electrolyte_s_avg)./F.^2 ; Phi_interface2 = Phi_interface1 - Kappa_eff_s./kappa_s.*log(Ce_interface2./Ce_interface1) ;
Kappa_eff_p = kappa_p.*RT.*tp_plus_bar.*(1+ dln_activity_coefficient_p_over_dln_C_electrolyte_p_avg)./F.^2 ;
doi:10.6342/NTU201602734
139
S_interface2 = iapp_p + 3.*Kappa_eff_p.*(C_electrolyte_p_avg./Ce_interface2 - 1)./Lp_electrolyte + Kappa_eff_p ; Phi_electrolyte_p_avg = S_interface2.*Lp_electrolyte./3 + Phi_interface2 ;
Phi_solid_p_avg = Phi_electrolyte_p_avg + Up + RT./(alpha.*F).*asinh(-Jp.*(2.*F.*kp).^-1.*
((C_solid_p_max - C_solid_p_surface).*C_solid_p_surface.*C_electrolyte_p_avg).^-alpha) ; Phi_solid_p_terminal = (-Lp_solid/3).*iapp_p./Modify_sigma_p + Phi_solid_p_avg ;
IR_collector = (R_collector_n_Cu + IR_contact_n).*abs(iapp_n) + (R_collector_p_Al + IR_contact_p).*abs(iapp_p) ; C_electrolyte_n_avg = (An.*C_electrolyte_n_avg + Bn.*C_electrolyte_p_avg
+ as_n.*tn_plus_bar.*Jn./F).*dt./Electorlyte_epsilon_n + C_electrolyte_n_avg ; C_electrolyte_p_avg = (Ap.*C_electrolyte_n_avg + Bp.*C_electrolyte_p_avg + as_p.*tp_plus_bar.*Jp./F).*dt./Electorlyte_epsilon_p + C_electrolyte_p_avg ; C_electrolyte_s_avg = (Ce_interface1 + Ce_interface2)./2 ;
Battery_Working_Potential = Phi_solid_p_terminal - Phi_solid_n_terminal - IR_collector ; Battery_n_SOC = U_SOC_n./Initial_Available_ratio_n ;
if T_constant == 0 如果關閉定溫功能，則進行發熱與溫度計算；
Heat = current.*(Un + Up - T_core.*(dUn_over_dT + dUp_over_dT).*0.001 + Phi_solid_n_terminal - Phi_electrolyte_n_avg + Phi_electrolyte_p_avg - Phi_solid_p_terminal + current.*
((IR_contact_n + R_collector_n_Cu).*Electrode_area_n + (IR_contact_p.*R_collector_p_Al).*Electrode_area_p));
T_core = (Heat - (Convection_coefficient.*Battery_surface_area ).*
(T_core - T_ambient)).*dt./(Battery_density_dot_Cp.*Battey_volume) + T_core ; end
Power_output = Power_output + current*Battery_Working_Potential ; 電池之輸出功率統計；
% % % % % % % % % % % % % % % % % % % % % % % % % %電池電位演算結束。
x1(1,t) = current ; 保存電流資料於x1中；
x2(1,t) = Battery_Working_Potential ; 保存電池工作電壓資料於x2中；
% % % % % % % % % % % % % % % % % % % % % % % % % %
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
doi:10.6342/NTU201602734
140
if imag(Battery_Working_Potential) ~= 0 || isnan(Battery_Working_Potential)
程式保護設定，當電池電壓出現虛部或非數字時
disp('Wrong Battery_Working_Potential') 顯示｢錯誤工作電壓｣
break 立刻結束運算 end 保護設定指令結束 end 時間演進邏輯指令結束 toc 運算計時結束
doi:10.6342/NTU201602734
141
因為有您在而分擔不少國慶老師給予的關注 壓力 ，十分謝謝您!此外，謝謝施華儒學長，雖然我都不 太把你當學長對待，三不五時跟您嬉鬧，還會口出酸言嗆聲兩句，但您總是笑而化之包容我。在此要 感謝您為 228 和平實驗室注入滿滿的活力，從中秋節烤肉發跡，接著是火鍋、烤雞塊、脆皮燒肉、親 子丼、屏東豬腳及雪花糕、醉雞、定期聚餐…等等，說不完的驚艷與有趣，還要謝謝您在我想偷懶不運 動時，拉著我一起健身，謝謝!另外，謝謝同窗好戰友翁偉澤，從每周的例行作戰會議，到最後的口試
因為有您在而分擔不少國慶老師給予的關注 壓力 ，十分謝謝您!此外，謝謝施華儒學長，雖然我都不 太把你當學長對待，三不五時跟您嬉鬧，還會口出酸言嗆聲兩句，但您總是笑而化之包容我。在此要 感謝您為 228 和平實驗室注入滿滿的活力，從中秋節烤肉發跡，接著是火鍋、烤雞塊、脆皮燒肉、親 子丼、屏東豬腳及雪花糕、醉雞、定期聚餐…等等，說不完的驚艷與有趣，還要謝謝您在我想偷懶不運 動時，拉著我一起健身，謝謝!另外，謝謝同窗好戰友翁偉澤，從每周的例行作戰會議，到最後的口試