Bugfix Concentration Indices, Updt WLTP example, Support Julia 1.10

BradyPlanden · Jan 6, 2024 · f80b764 · f80b764
1 parent a1ae688
commit f80b764
Show file tree

Hide file tree

Showing 4 changed files with 103 additions and 83 deletions.
diff --git 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
@@ -4,11 +4,16 @@ authors = ["Brady Planden <[email protected]>"]
 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
@@ -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ₛ)
 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, Ŝ)
-Results = WLTP(Cell, Ŝ, SOC, WLTP_P, A, B, C, D)
+Results, tₑ = WLTP(Cell, Ŝ, 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
@@ -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,
-               ϕ_ẽ1 = Array{Float64}(undef, tlength, size(ϕ_ẽInd, 1)) .= 0.0,
-               ϕ_ẽ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,
+        ϕ_ẽ1 = Array{Float64}(undef, tlength, size(ϕ_ẽInd, 1)) .= 0.0,
+        ϕ_ẽ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.ϕ_ẽ1[i + 1, :] = Results.y[i + 1, ϕ_ẽ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.α))) .*
-Original file line number
+Diff line change
@@ Expand Up / @@ -17,6 +17,7 @@ jobs: @@
               - '1.7'
               - '1.8'
               - '1.9'
+              - '1.10'
               - 'nightly'
         steps:
           - uses: actions/checkout@v2
@@ Expand Down @@