From f80b764c0fb1b1fbdc4e40f821f03721bd7dc58f Mon Sep 17 00:00:00 2001 From: Brady Planden Date: Sat, 6 Jan 2024 12:21:08 +0000 Subject: [PATCH] Bugfix Concentration Indices, Updt WLTP example, Support Julia 1.10 --- .github/workflows/CI.yml | 1 + Project.toml | 7 ++- examples/WLTP.jl | 94 ++++++++++++++++++++++----------------- src/Functions/Simulate.jl | 84 +++++++++++++++++----------------- 4 files changed, 103 insertions(+), 83 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 4c7a5a5..12efd4f 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -17,6 +17,7 @@ jobs: - '1.7' - '1.8' - '1.9' + - '1.10' - 'nightly' steps: - uses: actions/checkout@v2 diff --git a/Project.toml b/Project.toml index 6158c4b..f81c153 100644 --- a/Project.toml +++ b/Project.toml @@ -4,11 +4,16 @@ authors = ["Brady Planden "] version = "0.3.3" [deps] +Colors = "5ae59095-9a9b-59fe-a467-6f913c188581" FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" +LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +MAT = "23992714-dd62-5051-b70f-ba57cb901cac" +PGFPlotsX = "8314cec4-20b6-5062-9cdb-752b83310925" Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a" +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" TSVD = "9449cd9e-2762-5aa3-a617-5413e99d722e" @@ -23,4 +28,4 @@ Parameters = "0.12" Roots = "1.3, 2" TSVD = "0.4" UnitSystems = "0.3" -julia = "1.7, 1.8" +julia = "1.7, 1.8, 1.9, 1.10" diff --git a/examples/WLTP.jl b/examples/WLTP.jl index 94aa870..f5a9bea 100644 --- a/examples/WLTP.jl +++ b/examples/WLTP.jl @@ -3,7 +3,10 @@ plotly() default(show = true) #---------- Cell Definition -----------------# +Sₑ = 4 # Spatial points in electrolyte +Sₛ = 4 # Spatial point in solid Cell = Construct("LG M50") +Spatial!(Cell, Sₑ, Sₛ) Ŝ = collect(1.0:-1:0.0) SOC = 0.717 WLTP_File = matopen("examples/WLTP/WLTP_M50_M3.mat") @@ -12,50 +15,61 @@ WLTP_P = read(WLTP_File, "P_Models") #---------- Generate & Simulate Model -----------------# Cell.RA.H1 = Cell.RA.H2 = [1:2000; 3000:3500; 4000:4250] A, B, C, D = Realise(Cell, Ŝ) -Results = WLTP(Cell, Ŝ, SOC, WLTP_P, A, B, C, D) +Results, tₑ = WLTP(Cell, Ŝ, SOC, WLTP_P, A, B, C, D) #----------- Plotting ---------------------------# plot(Results.t, Results.Cell_V; - legend = :topright, - color = :blue, - bottom_margin = 5Plots.mm, - left_margin = 5Plots.mm, - right_margin = 15Plots.mm, - ylabel = "Terminal Voltage (V)", - xlabel = "Time (s)", - title = "WLTP Voltage", - label = "Voltage", - size = (1280, 720)) + legend = :topright, + color = :blue, + bottom_margin = 5Plots.mm, + left_margin = 5Plots.mm, + right_margin = 15Plots.mm, + ylabel = "Terminal Voltage (V)", + xlabel = "Time (s)", + title = "WLTP Voltage", + label = "Voltage", + size = (1280, 720)) plot(Results.t, Results.Ce; - legend = :topright, - bottom_margin = 5Plots.mm, - left_margin = 5Plots.mm, - right_margin = 15Plots.mm, - ylabel = "Electrolyte Concen. (mol/m³)", - xlabel = "Time (s)", - title = "Electrolyte Concentration", - label = ["Neg. Separator Interface" "Neg. Current Collector" "Pos. Current Collector" "Pos. Separator Interface"], - size = (1280, 720)) + legend = :topright, + bottom_margin = 5Plots.mm, + left_margin = 5Plots.mm, + right_margin = 15Plots.mm, + ylabel = "Electrolyte Concen. (mol/m³)", + xlabel = "Time (s)", + title = "Electrolyte Concentration", + label = ["Neg. Separator Interface" "Neg. Current Collector" "Pos. Current Collector" "Pos. Separator Interface"], + size = (1280, 720)) -plot(Results.t, Results.Cse_Pos; - legend = :topright, - bottom_margin = 5Plots.mm, - left_margin = 5Plots.mm, - right_margin = 15Plots.mm, - ylabel = "Concentration (mol/m³)", - xlabel = "Time (s)", - title = "Positive Electrode Concentration", - label = ["Current Collector" "Separator Interface"], - size = (1280, 720)) +plot(Results.t, Results.Cseₚ; + legend = :topright, + bottom_margin = 5Plots.mm, + left_margin = 5Plots.mm, + right_margin = 15Plots.mm, + ylabel = "Concentration (mol/m³)", + xlabel = "Time (s)", + title = "Positive Electrode Concentration", + label = ["Current Collector" "Separator Interface"], + size = (1280, 720)) -plot(Results.t, Results.Cse_Neg; - legend = :topright, - bottom_margin = 5Plots.mm, - left_margin = 5Plots.mm, - right_margin = 15Plots.mm, - ylabel = "Concentration (mol/m³)", - xlabel = "Time [s]", - title = "Negative Electrode Concentration", - label = ["Current Collector" "Separator Interface"], - size = (1280, 720)) +plot(Results.t, Results.Cseₙ; + legend = :topright, + bottom_margin = 5Plots.mm, + left_margin = 5Plots.mm, + right_margin = 15Plots.mm, + ylabel = "Concentration (mol/m³)", + xlabel = "Time [s]", + title = "Negative Electrode Concentration", + label = ["Current Collector" "Separator Interface"], + size = (1280, 720)) + +plot(Results.t, Results.Cell_SOC; + legend = :topright, + bottom_margin = 5Plots.mm, + left_margin = 5Plots.mm, + right_margin = 15Plots.mm, + ylabel = "–", + xlabel = "Time [s]", + title = "Cell SOC", + label = "Cell SOC", + size = (1280, 720)) diff --git a/src/Functions/Simulate.jl b/src/Functions/Simulate.jl index bdc32cb..a0d5d06 100644 --- a/src/Functions/Simulate.jl +++ b/src/Functions/Simulate.jl @@ -55,40 +55,40 @@ function Simulate(Cell, Input, Def, Tk, SList, SOC, A₀, B₀, C₀, D₀, t) CeSepInd = findall(CeSep .- CeNegOffset .== 1) CePosInd = findall(CePos .- CeSepOffset .== 1) - csegain_neg = C[CseNegInd[1][1], 1] #First Column in C Array (zeros column) - csegain_pos = C[CsePosInd[1][1], 1] #First Column in C Array (zeros column) + csegain_neg = C[CseNegInd[1][1], end] #End Column in C Array (zeros column) + csegain_pos = C[CsePosInd[1][1], end] #End Column in C Array (zeros column) # Memory Allocation Results = (θₙ = Array{Float64}(undef, tlength, 1) .= 0.0, - θₚ = Array{Float64}(undef, tlength, 1) .= 0.0, - jeqₙ = Array{Float64}(undef, tlength, 1) .= 0.0, - jeqₚ = Array{Float64}(undef, tlength, 1) .= 0.0, - x = Array{Float64}(undef, tlength + 1, size(A, 1)) .= 0.0, - y = Array{Float64}(undef, tlength, size(C, 1)) .= 0.0, - Cseₙ = Array{Float64}(undef, tlength, size(CseNegInd, 1)) .= 0.0, - Cseₚ = Array{Float64}(undef, tlength, size(CsePosInd, 1)) .= 0.0, - Ce = Array{Float64}(undef, tlength, size(CeInd, 1)) .= Cell.Const.ce0, - η₀ = Array{Float64}(undef, tlength, 1) .= 0.0, - ηₙ = Array{Float64}(undef, tlength, size(FluxNegInd, 1)) .= 0.0, - ηL = Array{Float64}(undef, tlength, 1) .= 0.0, - ηₚ = Array{Float64}(undef, tlength, size(FluxPosInd, 1)) .= 0.0, - ϕ_ẽ1 = Array{Float64}(undef, tlength, size(ϕ_ẽInd, 1)) .= 0.0, - ϕ_ẽ2 = Array{Float64}(undef, tlength, size(CeInd, 1)) .= 0.0, - ϕse₀ⁿ = Array{Float64}(undef, tlength, 1) .= 0.0, #Replace with length of ϕ_seNegInd @ zero - jₙ = Array{Float64}(undef, tlength, size(FluxNegInd, 1)) .= 0.0, - j₀ = Array{Float64}(undef, tlength, 1) .= 0.0, - jₚ = Array{Float64}(undef, tlength, size(FluxPosInd, 1)) .= 0.0, - jL = Array{Float64}(undef, tlength, 1) .= 0.0, - Rₜⁿ = Array{Float64}(undef, tlength, 1) .= 0.0, - Rₜᵖ = Array{Float64}(undef, tlength, 1) .= 0.0, - Uocpⁿ = Array{Float64}(undef, tlength, 1) .= 0.0, - Uocpᵖ = Array{Float64}(undef, tlength, 1) .= 0.0, - Cell_V = Array{Float64}(undef, tlength, 1) .= 0.0, - ϕ_e = Array{Float64}(undef, tlength, size(CeInd, 1)) .= 0.0, - Cell_SOC = Array{Float64}(undef, tlength, 1) .= 0, - Iapp = Array{Float64}(undef, tlength + 1, 1) .= 0, - t = t, - tₑ = Int64(1)) + θₚ = Array{Float64}(undef, tlength, 1) .= 0.0, + jeqₙ = Array{Float64}(undef, tlength, 1) .= 0.0, + jeqₚ = Array{Float64}(undef, tlength, 1) .= 0.0, + x = Array{Float64}(undef, tlength + 1, size(A, 1)) .= 0.0, + y = Array{Float64}(undef, tlength, size(C, 1)) .= 0.0, + Cseₙ = Array{Float64}(undef, tlength, size(CseNegInd, 1)) .= 0.0, + Cseₚ = Array{Float64}(undef, tlength, size(CsePosInd, 1)) .= 0.0, + Ce = Array{Float64}(undef, tlength, size(CeInd, 1)) .= Cell.Const.ce0, + η₀ = Array{Float64}(undef, tlength, 1) .= 0.0, + ηₙ = Array{Float64}(undef, tlength, size(FluxNegInd, 1)) .= 0.0, + ηL = Array{Float64}(undef, tlength, 1) .= 0.0, + ηₚ = Array{Float64}(undef, tlength, size(FluxPosInd, 1)) .= 0.0, + ϕ_ẽ1 = Array{Float64}(undef, tlength, size(ϕ_ẽInd, 1)) .= 0.0, + ϕ_ẽ2 = Array{Float64}(undef, tlength, size(CeInd, 1)) .= 0.0, + ϕse₀ⁿ = Array{Float64}(undef, tlength, 1) .= 0.0, #Replace with length of ϕ_seNegInd @ zero + jₙ = Array{Float64}(undef, tlength, size(FluxNegInd, 1)) .= 0.0, + j₀ = Array{Float64}(undef, tlength, 1) .= 0.0, + jₚ = Array{Float64}(undef, tlength, size(FluxPosInd, 1)) .= 0.0, + jL = Array{Float64}(undef, tlength, 1) .= 0.0, + Rₜⁿ = Array{Float64}(undef, tlength, 1) .= 0.0, + Rₜᵖ = Array{Float64}(undef, tlength, 1) .= 0.0, + Uocpⁿ = Array{Float64}(undef, tlength, 1) .= 0.0, + Uocpᵖ = Array{Float64}(undef, tlength, 1) .= 0.0, + Cell_V = Array{Float64}(undef, tlength, 1) .= 0.0, + ϕ_e = Array{Float64}(undef, tlength, size(CeInd, 1)) .= 0.0, + Cell_SOC = Array{Float64}(undef, tlength, 1) .= 0, + Iapp = Array{Float64}(undef, tlength + 1, 1) .= 0, + t = t, + tₑ = Int64(1)) # Defining SOC SOCₙ = SOC * (Cell.Neg.θ_100 - Cell.Neg.θ_0) + Cell.Neg.θ_0 @@ -99,10 +99,10 @@ function Simulate(Cell, Input, Def, Tk, SList, SOC, A₀, B₀, C₀, D₀, t) # Loop through time - compute dependent variables (voltage, flux, etc.) # for i in 0:(tlength - 1) - cs_neg_avg = Results.x[i + 1, 1] * csegain_neg + SOCₙ * Cell.Neg.cs_max < 0.0 ? - 0.0 : Results.x[i + 1, 1] * csegain_neg + SOCₙ * Cell.Neg.cs_max #Zero if < 0 - cs_pos_avg = Results.x[i + 1, 1] * csegain_pos + SOCₚ * Cell.Pos.cs_max < 0.0 ? - 0.0 : Results.x[i + 1, 1] * csegain_pos + SOCₚ * Cell.Pos.cs_max #Zero if < 0 + cs_neg_avg = Results.x[i + 1, end] * csegain_neg + SOCₙ * Cell.Neg.cs_max < 0.0 ? + 0.0 : Results.x[i + 1, end] * csegain_neg + SOCₙ * Cell.Neg.cs_max #Zero if < 0 + cs_pos_avg = Results.x[i + 1, end] * csegain_pos + SOCₚ * Cell.Pos.cs_max < 0.0 ? + 0.0 : Results.x[i + 1, end] * csegain_pos + SOCₚ * Cell.Pos.cs_max #Zero if < 0 if cs_neg_avg > Cell.Neg.cs_max cs_neg_avg = Cell.Neg.cs_max @@ -170,7 +170,7 @@ function Simulate(Cell, Input, Def, Tk, SList, SOC, A₀, B₀, C₀, D₀, t) # Relinearise dependent on ν, σ, κ D = D_Linear(Cell, ν_neg, ν_pos, σ_eff_Neg, κ_eff_Neg, σ_eff_Pos, κ_eff_Pos, - κ_eff_Sep) + κ_eff_Sep) # Interpolate C & D Matrices C = interp(C₀, SList, Results.Cell_SOC[i + 1]) @@ -196,11 +196,11 @@ function Simulate(Cell, Input, Def, Tk, SList, SOC, A₀, B₀, C₀, D₀, t) # Potentials Results.Uocpⁿ[i + 1] = Cell.Const.Uocp("Neg", - Results.Cseₙ[i + 1, 1] / - Cell.Neg.cs_max) + Results.Cseₙ[i + 1, 1] / + Cell.Neg.cs_max) Results.Uocpᵖ[i + 1] = Cell.Const.Uocp("Pos", - Results.Cseₚ[i + 1, 1] / - Cell.Pos.cs_max) + Results.Cseₚ[i + 1, 1] / + Cell.Pos.cs_max) Results.ϕse₀ⁿ[i + 1] = Results.y[i + 1, ϕ_seNegInd[1]] + Results.Uocpⁿ[i + 1] #Location 0 Results.ϕ_ẽ1[i + 1, :] = Results.y[i + 1, ϕ_ẽInd] @@ -219,7 +219,7 @@ function Simulate(Cell, Input, Def, Tk, SList, SOC, A₀, B₀, C₀, D₀, t) # Neg j₀ⁿ = findmax([ones(size(Results.Cseₙ, 2)) * eps() (kₙ .* (Results.Cseₙ[i + 1, - :] .^ + :] .^ Cell.Neg.α) .* (Results.Ce[i + 1, 1] .^ (1 - Cell.Neg.α))) .* @@ -234,7 +234,7 @@ function Simulate(Cell, Input, Def, Tk, SList, SOC, A₀, B₀, C₀, D₀, t) # Pos j₀ᵖ = findmax([ones(size(Results.Cseₚ, 2)) * eps() (kₚ .* (Results.Cseₚ[i + 1, - :] .^ + :] .^ Cell.Pos.α) .* (Results.Ce[i + 1, 1] .^ (1 - Cell.Pos.α))) .*