diff --git a/Project.toml b/Project.toml index 5aec7dd..e93134c 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ElectrochemicalKinetics" uuid = "a2c6e634-85ca-418a-9c67-9b5417ce2d04" authors = ["Rachel Kurchin ", "Holden Parks ", "Dhairya Gandhi "] -version = "0.1.4" +version = "0.1.5" [deps] DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" diff --git a/README.md b/README.md index 696ba7c..ebdf329 100644 --- a/README.md +++ b/README.md @@ -71,6 +71,7 @@ models = [fit_model(exp_data, model_type; dos_file=dosfile) for model_type in mo plot_exp_and_models(exp_data, models) ``` The resulting plot looks like: + models_expts But we might also want to see the higher-voltage behavior and compare the models absent the experimental data. There's another function for that (by default it plots to +/- 1V)... @@ -86,7 +87,7 @@ For more on this, see upcoming paper: [arxiv link placeholder] ### Butler-Volmer Probably the most basic (and largely empirical) kinetic model. The rate constants are given by: -$$k_{\text{ox}}^{\text{BV}}(\eta)=A\exp\left(\frac{\alpha\eta}{k_{\text B}T}\right)\\k_{\text{red}}^{\text{BV}}(\eta)=A\exp\left(\frac{(1-\alpha)\eta}{k_{\text B}T}\right)$$ +$$k_{\text{ox}}^{\text{BV}}(\eta)=A\exp\left(\frac{\alpha\eta}{k_{\text B}T}\right)\\ k_{\text{red}}^{\text{BV}}(\eta)=A\exp\left(\frac{(1-\alpha)\eta}{k_{\text B}T}\right)$$ It is implemented in the package as `ButlerVolmer`, possessing two fields: the prefactor `A` and the transfer coefficient `α`. diff --git a/src/phase_diagrams.jl b/src/phase_diagrams.jl index e39cd5d..6178ca4 100644 --- a/src/phase_diagrams.jl +++ b/src/phase_diagrams.jl @@ -43,6 +43,7 @@ function g_kinetic(I, km::KineticModel; intercalate=true, warn=true, T=room_T, k return g end +# TODO: figure out T passthrough issue # zeros of this function correspond to pairs of x's satisfying the common tangent condition for a given µ function # case where we just feed in two points (x should be of length two) function common_tangent(x::Vector, I, km::KineticModel; intercalate=true, kwargs...) @@ -72,18 +73,22 @@ NOTE 1: appropriate values of `I_step` depend strongly on the prefactor of your NOTE 2: at lower temperatures (<=320K or so), ButlerVolmer models with the default thermodynamic parameters have a two-phase region at every current, so setting a finite value of I_max is necessary for this function to finish running. """ -function phase_diagram(km::KineticModel; I_start=0.0, I_step=1.0, I_max=Inf, verbose=false, intercalate=true, start_guess=[0.05, 0.95], kwargs...) +function phase_diagram(km::KineticModel; I_start=0.0, I_step=1.0, I_max=Inf, verbose=false, intercalate=true, start_guess=[0.05, 0.95], tol=5e-3, kwargs...) I = I_start pbs_here = find_phase_boundaries(I, km; intercalate=intercalate, guess=start_guess, kwargs...) pbs = pbs_here' I_vals = [I_start] - while abs(pbs_here[2] - pbs_here[1]) > 1e-3 && I < I_max + pb_diff = [0.0, 0.0] + while abs(pbs_here[2] - pbs_here[1]) > tol && I < I_max I = I + I_step if verbose println("Solving at I=", I, "...") end try - pbs_here = find_phase_boundaries(I, km; intercalate=intercalate, guess=pbs_here, kwargs...) + pbs_old = pbs_here + pbs_here = find_phase_boundaries(I, km; intercalate=intercalate, guess=pbs_old .+ pb_diff, kwargs...) + pb_diff = pbs_here .- pbs_old + # TODO: check that they haven't crossed over pbs = vcat(pbs, pbs_here') push!(I_vals, I) if verbose diff --git a/test/phase_diagram_tests.jl b/test/phase_diagram_tests.jl index ef9c0f4..5a56988 100644 --- a/test/phase_diagram_tests.jl +++ b/test/phase_diagram_tests.jl @@ -117,24 +117,60 @@ end end @testset "Phase Diagram" begin - km = bv # TODO: expand this - - # simplest case, just one pair of x values (this function is still pretty slow though) - v1 = find_phase_boundaries(100, km) - - @test all(isapprox.(common_tangent(v1, 100, km), Ref(0.0), atol=1e-6)) - v2 = find_phase_boundaries(100, km, T=350) - @test all(isapprox.(common_tangent(v2, 100, km, T=350), Ref(0.0), atol=1e-5)) - # they should get "narrower" with temperature - @test v2[1] > v1[1] - @test v2[2] < v1[2] - v3 = find_phase_boundaries(400, km) - # ...and also with current - @test v3[1] > v1[1] - @test v3[2] < v1[2] - - # test actual numerical values too - @test all(isapprox.(v1, [0.045547, 0.841501], atol=1e-5)) - @test all(isapprox.(v2, [0.090478, 0.77212], atol=1e-4)) - @test all(isapprox.(v3, [0.105747, 0.6809], atol=1e-4)) + v1_vals = Dict(:ButlerVolmer=>[0.045547, 0.841501], + :Marcus=>[0.04901, 0.81884], + :AsymptoticMarcusHushChidsey=>[0.04958, 0.81897]) + v2_vals = Dict(:ButlerVolmer=>[0.090478, 0.77212], + :Marcus=>[0.078934, 0.804275], + :AsymptoticMarcusHushChidsey=>[0.07609, 0.81652]) + v3_vals = Dict(:ButlerVolmer=>[0.105747, 0.6809], + :Marcus=>[0.154934, 0.54218], + :AsymptoticMarcusHushChidsey=>[0.147034,0.566617]) + I_max_T300 = Dict(:ButlerVolmer=>700.0, + :Marcus=>600.0, + :AsymptoticMarcusHushChidsey=>650.0) + I_max_T400 = Dict(:ButlerVolmer=>250.0, + :Marcus=>400.0, + :AsymptoticMarcusHushChidsey=>500.0) + for km in kms + model_type = nameof(typeof(km)) + @testset "$model_type" begin + # simplest case, just one pair of x values (this function is still pretty slow though) + v1 = find_phase_boundaries(100, km) + + @test all(isapprox.(common_tangent(v1, 100, km), Ref(0.0), atol=1e-6)) + v2 = find_phase_boundaries(100, km, T=350) + @test all(isapprox.(common_tangent(v2, 100, km, T=350), Ref(0.0), atol=1e-5)) + # they should get "narrower" with temperature + @test v2[1] > v1[1] + @test v2[2] < v1[2] + v3 = find_phase_boundaries(400, km, guess=[0.1,0.8]) + # ...and also with current + @test v3[1] > v1[1] + @test v3[2] < v1[2] + + # test actual numerical values too + @test all(isapprox.(v1, v1_vals[model_type], atol=1e-5)) + @test all(isapprox.(v2, v2_vals[model_type], atol=1e-4)) + @test all(isapprox.(v3, v3_vals[model_type], atol=1e-4)) + + # and actually build the full map for a couple temps + pbs, I = phase_diagram(km, I_max=700, I_step=50) + @test maximum(I) == I_max_T300[model_type] + end + end + + @testset "MHC" begin + # do some lighterweight stuff here since this takes awhile + mhc = MarcusHushChidsey(amhc.A, amhc.λ) + # also give it a good initial guess so it only takes one or two optimizer steps + v1 = find_phase_boundaries(100, mhc, guess=v1_vals[:AsymptoticMarcusHushChidsey]) + @test all(isapprox.(v1, [0.04967204036, 0.81822937543], atol=1e-5)) + + v2 = find_phase_boundaries(100, mhc, T=350, guess=v2_vals[:AsymptoticMarcusHushChidsey]) + @test all(isapprox.(v2, [0.075615, 0.8171636], atol=1e-5)) + + v3 = find_phase_boundaries(400, mhc, guess=v3_vals[:AsymptoticMarcusHushChidsey]) + @test all(isapprox.(v3, [0.1467487, 0.5704125], atol=1e-5)) + end end