-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathconstant_matrices.m
More file actions
112 lines (95 loc) · 7.9 KB
/
Copy pathconstant_matrices.m
File metadata and controls
112 lines (95 loc) · 7.9 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
function [ matrices ] = constant_matrices( params )
%MATRICES Summary of this function goes here
% Detailed explanation goes here
%% Basic differentiation matrices for each electrode
matrices.n = constant_matrices_single_electrode(params.n);
matrices.p = constant_matrices_single_electrode(params.p);
%% Basic differentiation matrices for the separator
[matrices.s.adim.D1_noBC,matrices.s.adim.x] = cheb(params.s.dscrtzn.N_e);
matrices.s.adim.D2_noBC = matrices.s.adim.D1_noBC ^2;
%% Boundary operators
matrices.bdry_operators = bdry_operators(...
matrices.p.adim.D1_noBC(end,:),...
matrices.p.adim.D1_noBC(1,:),...
matrices.s.adim.D1_noBC(1,:),...
matrices.s.adim.D1_noBC(end,:),...
matrices.n.adim.D1_noBC(end,:),...
matrices.n.adim.D1_noBC(1,:),...
params...
);
end
function [matrices_single_electrode] = constant_matrices_single_electrode(params_electrode)
matrices_single_electrode.N_s = params_electrode.dscrtzn.N_s;
matrices_single_electrode.N_elyte = params_electrode.dscrtzn.N_e;
L = params_electrode.L;
Rs = params_electrode.R_s;
%% Scaled differentiation operators (no Boundary Conditions)
[ephase.D1_cheb,ephase.x_cheb] = cheb(matrices_single_electrode.N_elyte);
[sphase.D1_cheb,sphase.r_cheb] = cheb(matrices_single_electrode.N_s);
% For the dimensional system
matrices_single_electrode.ephase.D1_noBC = (2/L)*ephase.D1_cheb;
matrices_single_electrode.sphase.D1_noBC = (2/Rs)*sphase.D1_cheb;
matrices_single_electrode.ephase.x = (1 + ephase.x_cheb)*L/2;
matrices_single_electrode.sphase.r = (1 + sphase.r_cheb) * Rs/2;
% For the adimensional system
matrices_single_electrode.adim.ephase.D1_noBC = 2 * ephase.D1_cheb;
matrices_single_electrode.adim.sphase.D1_noBC = 2* sphase.D1_cheb;
matrices_single_electrode.adim.ephase.D2_noBC = matrices_single_electrode.adim.ephase.D1_noBC ^ 2;
matrices_single_electrode.adim.sphase.D2_noBC = matrices_single_electrode.adim.sphase.D1_noBC ^ 2;
matrices_single_electrode.adim.ephase.x = (1 + ephase.x_cheb)/2;
matrices_single_electrode.adim.sphase.r = (1 + sphase.r_cheb)/2;
%% Boundary conditions : electrolyte phase
% Getting the modifications that should be made to the _noBC matrices to get
% the _BC matrices :
[matrices_single_electrode.adim.ephase.elyte.Dmod,matrices_single_electrode.adim.ephase.elyte.Emod,matrices_single_electrode.adim.ephase.elyte.Umod,matrices_single_electrode.adim.ephase.elyte.Vmod] = ...
matrices_bdrys_wrt_internal(...
matrices_single_electrode.adim.ephase.D1_noBC(1,:),...
matrices_single_electrode.adim.ephase.D1_noBC(end,:),...
0,...
0,...
[0,0]...
);
% Remark : since there is no exogeneous boundary condition, matrices.adim.ephase.elyte.Umod and matrices.adim.ephase.elyte.Vmod are 0.
% Actually computing the _BC matrix, and the all the terms on the right side of dc_e/dt = ... and ddl/dt = ... :
% These terms are
% dc_e/dt = params.adim.theta_d * (elyte.D2_BC * ce + elyte.exogeneous_bdry ) ...
% + params.adim.theta_c * (dl.D2_BC * dl + dl.exogeneous_bdry) ...
% + params.adim.K * params.adim.theta_c *elyte.logterm_BC(ce)
% ddl/dt = params.adim.theta_c * (dl.D2_BC * dl + dl.exogeneous_bdry) ...
% - lambda(dl,ussurf) ...
% + params.adim.K * params.adim.theta_c *elyte.logterm_BC(ce)
% The matrices depending on dl depend on time (because they depend on i(t)),
% and then are defined in the file matrices_linearized_t_dep.m
matrices_single_electrode.adim.ephase.elyte.D2_BC = matrices_single_electrode.adim.ephase.D2_noBC(2:end-1,2:end-1)...
+ matrices_single_electrode.adim.ephase.D2_noBC(2:end-1,1) * matrices_single_electrode.adim.ephase.elyte.Dmod...
+ matrices_single_electrode.adim.ephase.D2_noBC(2:end-1,end) * matrices_single_electrode.adim.ephase.elyte.Emod;
matrices_single_electrode.adim.ephase.elyte.exogeneous_bdry_left = matrices_single_electrode.adim.ephase.D2_noBC(2:end-1,1) * matrices_single_electrode.adim.ephase.elyte.Umod; % This is 0 bcz we do not have an exogeneous term at the bdry for ce, but we still define it to have consistent notations
matrices_single_electrode.adim.ephase.elyte.exogeneous_bdry_right = matrices_single_electrode.adim.ephase.D2_noBC(2:end-1,end) * matrices_single_electrode.adim.ephase.elyte.Vmod; % Same
matrices_single_electrode.adim.ephase.elyte.exogeneous_bdry = matrices_single_electrode.adim.ephase.elyte.exogeneous_bdry_left + matrices_single_electrode.adim.ephase.elyte.exogeneous_bdry_right; % Same
matrices_single_electrode.adim.ephase.elyte.logterm_BC = ...
@(ELYTE2N) (matrices_single_electrode.adim.ephase.D2_noBC(2:end-1,2:end-1) * log(ELYTE2N)...
+ matrices_single_electrode.adim.ephase.D2_noBC(2:end-1,1) * log(matrices_single_electrode.adim.ephase.elyte.Umod + matrices_single_electrode.adim.ephase.elyte.Dmod * ELYTE2N)...
+ matrices_single_electrode.adim.ephase.D2_noBC(2:end-1,end) * log(matrices_single_electrode.adim.ephase.elyte.Vmod + matrices_single_electrode.adim.ephase.elyte.Emod * ELYTE2N)...
);
%% Boundary conditions : solid phase
%The operators defined below give the dynamics for u_s, with automatically
%satisfied boundary conditions. The equation is (k is a parameter that designs the x-absciss, k is in [2, N_elyte]:
% du_s/dt = D2_BC us + exogeneous_linear_dl * dl(k)...
% + exogeneous_nonlinear(dl(k),us)
matrices_single_electrode.adim.sphase.us.D2_BC = matrices_single_electrode.adim.sphase.D2_noBC(2:end-1,2:end-1);
matrices_single_electrode.adim.sphase.us.exogeneous_nonlinear = ...
@(USSURF) (matrices_single_electrode.adim.sphase.D2_noBC(2:end-1,1) * USSURF');
%matrices.adim.sphase.us.compute_ussurf = ...
% @(DLK,US2R) (DLK - numeric_inverse_x_min_sinh(params.misc.xmin,params.misc.xmax,coef,...
% DLK...
% + matrices.adim.sphase.D1_noBC(1,2:end-1) * US2R / (matrices.adim.sphase.D1_noBC(1) - 1)...
% )...
% );
%% Algebraic equation
matrices_single_electrode.adim.sphase.ussurf.linear_term_ussurf = - eye(matrices_single_electrode.N_elyte - 1);
matrices_single_electrode.adim.sphase.ussurf.linear_term_us = kron(eye(matrices_single_electrode.N_elyte - 1),matrices_single_electrode.adim.sphase.D1_noBC(1,2:end-1) / (1 - matrices_single_electrode.adim.sphase.D1_noBC(1)));
matrices_single_electrode.adim.ephase.dl.butler_volmer = ...
@(CE,DL2N, USSURF2N) (...
((CE .* (1 - USSURF2N) .* USSURF2N) .^ params_electrode.alpha) .* 2.*sinh(params_electrode.alpha * params_electrode.E * (DL2N - ocp_dualfoil(USSURF2N)/params_electrode.V_0) )...
);
end