% Program_Level = Ultimate*******_Driving_Cycle_Test_with_Constant_Current_||_Constant_Voltage_Charger ;
程式等級﹕極致級*******，能估算大規模串、並聯電池組之老化程度，並且擁有〈定電流 − 定電壓〉之充電功能，可用
來估算電動車電池組於道路駕駛時之狀態及未來老化走勢；
% Default_Battery_Type = Panasonic_NCR18650_3100mAh_NCA_LixC6 ;
預設電池型號﹕松下Panasonic_NCR18650_3100mAh，鎳鈷鋁氧化物/鋰碳化合物；
% Default_Vehicle = Telsa_Model_S_90kWh
預設電動車之車型﹕Telsa_Model_S_90kWh；
% Default_Driving_Cycle = Common_ARTEMIS_Driving_Cycles
預設駕駛路況﹕Common_ARTEMIS_Driving_Cycles；
% 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/5 ; 設定充放電速率；
doi:10.6342/NTU201602734
143
Cycles = 1 ; 設定循環次數；
Lower_Cut_off_Voltage = 320 ; 最低截止電壓；
Upper_Cut_off_Voltage = 410 ; 最高截止電壓；
Guess_Tolerance = 0.1 ; 設置定電壓充電時，其電壓誤差之容忍度；
Current_Adjust = 0.0005; 設置定電壓充電時，其電流之最小調控精度；
Shen_Ji_Miao_Suan = 1.001 ; 設置定電壓充電時，猜測電流之起始值； ｢神機妙算｣的漢語拼音
Start_CV_Process = 400 ; 定電壓充電啟動時機；
Constant_Voltage = 400 ; 定電壓充電值；
Cut_off_Current = -0.155*75 ; 定電壓充電結束時機；
Capacity_Test = 0 ; 電容量測定功能，1表示啟動，0代表關閉；
First_capacity_test_cycles = 1 ; 第一次電容量測定時機；
capacity_test_per_cycles = 50 ; 電容量測定間隔；
% ________________________________________________
T_ambient = 298.15 ; 環境溫度；
Convection_coefficient = 10 ; 對流係數；
% ________________________________________________
Balance_Weight = 0.08 ; 並聯電池組之平衡電流權重；
Series_Number = 100 ; 串聯數目；
Parallel_Number = 75 ; 並聯數目；
% ________________________________________________
Manufacture_Deviation = 0.05 ; 電池之製造偏差，0.05表示±5%；
Initial_average_SOC = 0.9 ; 電池單元起始蓄電狀態，0.9表示SOC = 90%；
SOC_Range = 0.1 ; 電池組內部之電池單元起始蓄電狀態差異，0.1表示對Initial_average_SOC ± 10%；
% ________________________________________________
load cycles_artemis.mat 讀取ARTEMIS駕駛路況之資料；
Driving_Cycle_Type = [ ArtRoad.data ; ArtMw150.data ; ArtUrban.data ] ; 組合ARTEMIS駕駛路況；
doi:10.6342/NTU201602734
144
Dr_length = length(Driving_Cycle_Type) ; 組合駕駛路況之矩陣向量長度；
Gravity = 9.8 ; 重力加速度；
Air_density = 1.2 ; 空氣密度；
Road_inclination = 0 ; 路面坡度；
Inertia_factor = 1.1 ; 慣性因子；
Carring_mass = 65 ; 載運質量；
Vehicle_mass = 2144 ; 電動車質量；
Cross_sectional_area = 2.34 ; 電動車的風阻面積；
Drag_coefficient = 0.24 ; 空氣阻力係數；
Rolling_coefficient = 0.012 ; 滾動阻力係數；
Efficient_work = 0.8 ; 電池組管理系統之能源轉換效率；
Efficient_braking = 0.5 ; 發電機之回充效率；
Power_consumption_air_condition = 2000 ; 車內空調之耗能；
Power_consumption_cooling_system = 1000 ; 電池組冷卻系統之耗能；
Power_consumption_management_computer = 500 ; 行車控制系統之耗能；
Power_consumption_sum_others = Power_consumption_air_condition + Power_consumption_cooling_system + Power_consumption_management_computer ; 上述電器設備耗能總和；
Mass_total = (Vehicle_mass + Carring_mass)* Inertia_factor ; 整體總質量；
Weight_total = Mass_total*Gravity ; 整體總重量；
Vehicle_velocity = Driving_Cycle_Type./3.6 ; 行車速率單位換算；
Vehicle_velocity(Dr_length+1) = Vehicle_velocity(Dr_length) ; 假設駕駛路況結束後的下一秒，該車之速率與最後一秒相同；
% ________________________________________________
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%
Standard_Ah = 3.1 ; 標準安時數；
Battery_Type_C_solid_n_max = 28000 ; 負電極最大之鋰原子濃度；
doi:10.6342/NTU201602734
145
Battery_Type_C_solid_p_max = 51000 ; 正電極最大之鋰原子濃度；
C_solid_n_max_0 = 28000 ; 負電極最大之鋰原子濃度初始值；
C_solid_p_max_0 = 51000 ; 正電極最大之鋰原子濃度初始值；
Initial_Available_ratio_n = 0.63 ; 負電極最大之鋰原子濃度，其放電初始之可使用比例；
Initial_Available_ratio_p = 0.395 ; 正電極最大之鋰原子濃度，其放電初始之可使用比例；
Battery_Type_C_electrolyte_n_avg = 1200 ; 負極區電解液濃度；
Battery_Type_C_electrolyte_s_avg = 1200 ; 隔離膜區電解液濃度；
Battery_Type_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.*ones(Parallel_Number,Series_Number) ; 每個電池單元其負電極之顆粒半徑；
Rp = 4*10^-9.*ones(Parallel_Number,Series_Number) ; 每個電池單元其正電極之顆粒半徑；
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 ; 正電極長度；
dln_activity_coefficient_n_over_dln_C_electrolyte_n_avg = 1 ; = 1；
doi:10.6342/NTU201602734
146
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.*ones(Parallel_Number,Series_Number) ; 每個電池單元其負電極之體積分率；
Solid_epsilon_p = 0.4483.*ones(Parallel_Number,Series_Number) ; 每個電池單元其正電極之體積分率；
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指數；
doi:10.6342/NTU201602734
147
Bruggman_coefficient_s = 1.5 ; 隔離膜區的Bruggman指數；
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 ; 正極之單位面積接觸電阻值；
Y = 0.68 ; 電池老化經驗式中的γ；
U_SEI = 0.4 ; SEI膜之開路電壓；
M_SEI = 0.162 ; SEI膜之莫耳質量；
Density_SEI = 1690 ; SEI膜之密度；
R_SEI = 0.09.*ones(Parallel_Number,Series_Number) ; 每個電池單元其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 ; 電池表面積；
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
I = 0 ; 起始電池組的工作總電流；
t = 0 ; 初始時序；
t_Drive = 0 ; 駕駛路況用之時序；
CV = 0 ; 初始定電壓功能為關閉；
cyc = 1 初始循環序號；
doi:10.6342/NTU201602734
148
check = 0 ; 初始電容量測定功能為關閉；
Go_home = 0 初始位置在家中，因此初始設定值為0
Go_company = 1 從家裡出發去公司上班，因此初始設定值為1
Power_output = 0 ; 初始之電池組輸出功率值；
Start_charged = 0 ; 初始充電功能設為關閉；
dt = 1 ; 時間步階長度；
F = 96487 ; 法拉第常數；
R = 8.314 ; 理想氣體常數；
Delta_n = Ln_electrolyte./3 ; 負極區電解液之微量擴散長度；
Delta_p = Lp_electrolyte./3 ; 正極區電解液之微量擴散長度；
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.*ones(Parallel_Number,Series_Number) ; 電化學參數之參考溫度標準；
T_ambient = T_ambient.*ones(Parallel_Number,Series_Number) ; 初始電池溫度等於環境溫度；
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 ; 負電極單位長度電導值之多孔性修正；
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 ; 負電極每單位面積之有效表面積；
doi:10.6342/NTU201602734
149
tp_plus_bar = tp_plus - 1 ; 正極區電解液之轉移係數減一；
tn_plus_bar = tn_plus - 1 ; 負極區電解液之轉移係數減一；
ts_plus_bar = ts_plus - 1 ; 隔離膜區電解液之轉移係數減一；
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) ; 特徵值矩陣其行向量之長度；
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Minions_Kevin = ones(Parallel_Number,Series_Number) ;
預先建立內容全是1之矩陣，待後續運算使用，避免於時間演進時重複開儲位之動作，以加速運算。其中Minions_Kevin意即 小小兵中的凱文，是身長高的小小兵類型，像數字1一樣；
Minions_Bob = zeros(Parallel_Number,Series_Number) ;
預先建立內容全是0之矩陣，待後續運算使用，避免於時間演進時重複開儲位之動作，以加速運算。其中Minions_Bob意即 小小兵中的鮑伯，是身長矮的小小兵類型，像數字0一樣；
current = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
C_solid_n_max = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
C_solid_p_max = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
C_solid_n_surface = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
C_solid_p_surface = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
C_electrolyte_n_avg = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
C_electrolyte_s_avg = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
C_electrolyte_p_avg = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
Ah_Discharge_plus_Charge = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
Series_Voltage = zeros(Parallel_Number,1) ; 預先配置記憶體；
doi:10.6342/NTU201602734
150
Balance_Current = zeros(Parallel_Number,1) ; 預先配置記憶體；
Equal_Resistance = zeros(Parallel_Number,1) ; 預先配置記憶體；
power = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
iapp_n = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
iapp_p = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
Jn = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
Jp = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
IR_SEI = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
RT = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
R_diff_T = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
D_solid_n = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
D_solid_p = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
D_electrolyte_n = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
D_electrolyte_p = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
D_electrolyte_s = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
Diffusion_length_n = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
Diffusion_length_p = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
U_SOC_n = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
U_SOC_p = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
dUn_over_dT = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
dUp_over_dT = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
Un = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
Up = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
Alpha_n = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
Beta_sep = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
Gamma_p = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
doi:10.6342/NTU201602734
151
Omega = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
GAO = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
GA_over_GAO = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
An = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
Bn = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
Ap = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
Bp = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
kn = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
kp = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
kappa_n = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
kappa_s = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
kappa_p = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
Kappa_eff_n = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
Kappa_eff_s = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
Kappa_eff_p = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
Ce_interface1 = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
Ce_interface2 = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
Phi_solid_n_avg = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
Phi_solid_n_terminal = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
Phi_electrolyte_n_avg = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
S_interface1 = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
S_interface2 = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
Phi_interface1 = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
Phi_interface2 = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
Phi_solid_p_avg = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
Phi_solid_p_terminal = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
doi:10.6342/NTU201602734
152
Phi_electrolyte_p_avg = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
IR_collector = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
Battery_Working_Potential = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
Heat = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；
Heat = zeros(Parallel_Number,Series_Number) ; 預先配置記憶體；