-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathPNP_final_project_notebook_and_report.jl
More file actions
2904 lines (2336 loc) · 118 KB
/
Copy pathPNP_final_project_notebook_and_report.jl
File metadata and controls
2904 lines (2336 loc) · 118 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
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
### A Pluto.jl notebook ###
# v0.19.11
using Markdown
using InteractiveUtils
# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
quote
local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
local el = $(esc(element))
global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
el
end
end
# ╔═╡ 2eae94d7-4830-4a94-95d4-abbabef7a044
# ╠═╡ show_logs = false
using Printf,PlutoUI,Plots,ExtendableGrids,SimplexGridFactory,VoronoiFVM,GridVisualize,PyPlot ,Triangulate ,LaTeXStrings
# ╔═╡ 470fc128-929b-41fc-b5ce-ff6449ec3111
md"""
# Numerically Solving the One-Dimenional Poisson-Nernst-Planck Equations Using the Finite Volume Method
### Course Project in Scientific Computing at TU Berlin Aditya Kumar and Kay Töpfer
"""
# ╔═╡ 15a44712-fef6-4160-8087-4927b8bbad07
md"""
###
"""
# ╔═╡ aef4afe5-f7b1-4946-9592-42f065b856ce
md"""
## Introduction
"""
# ╔═╡ a0fefd44-853b-4bb9-94fa-8f7bedeacfb1
md"""
Charge transport is ubiquitous in various applications, such as semiconductors, micro-and-nano-fluidics, and biological ion channels. Charge transportation can be modelled by the Poisson-Nernst-Planck (PNP) equations—a system of coupled nonlinear partial differential equations. The PNP equations describe the diffusion of charged particles (ions) in solution resulting from an electric field, so-called electrodiffusion. The PNP equations can be derived from the gradient flow of an electrostatic free energy, which is based on mean-field approximations of chemical ionic species interactions. Furthermore, the PNP model provides continuum descriptions of chemical species concentration (with the corresponding chemical potential of the species) and the electrostatic potential.
In the absence of chemical reactions, the two mechanisms that lead to flow of chemical species into or out of a control volume for our problem are diffusion and electromigration—this defines our flux. The electromigration term is defined by coupling the chemical potential generated by the movement of the charged species concentrations to the applied electrostatic field represented by the Poisson equation. This movement also affects the solute in which the charged species live in, i.e. the neutral solvent. Thus the solvent is transported (migrated) as well. The Nernst-Planck equation is an extension of Fick’s law describing the migration of charged chemical particles in a dilute solution. Fick’s law defines a flux density of species proportional to the gradient of the chemical species concentration and the diffusivity of the species in the neutral solvent. Fick’s law is a macroscopic way of representing the summed effect of the random motion of species owing to thermal fluctuations. As a result, we see that the Nernst-Planck equation forms a continuity equation for the time-dependent species concentration that gives us a conservation law for the underlying physical phenomena of electrodiffusion in dilute solutions. The Poisson and the Nernst-Planck equations together form a system of coupled partial differential equations. Whereas the theory of the Poisson equation is mainly rooted in electrostatics, the Nernst-Planck equation is an incarnation of a conservation law that is only valid for dilute solutions.
The combination of a mean-field approximation ion interactions and continuum theoretic description of chemical properties of the ions provides both quantitative and qualitative predictions of experimental measurements of chemical species transport problems in many areas of physics and electrical engineering involving semiconductor development, micro-and-nanofluidic device construction, and biologically relevant systems. Although a useful model, there are some limitations. In particular, the number of equations to be solved and the number of diffusion coefficient profiles to be determined for the calculation directly depend on the number of ion species in the system, since each ion species corresponds to one Nernst-Planck equation and one position-dependent diffusion coefficient profile. This leads to the numerical solution of PNP becoming computationally expensive. In this course project, the amount of charged species interacting is two, along with a neutral solvent. Finite-volumes of the charged species and solvent will be taken into account.
"""
# ╔═╡ 699f7fac-7a4b-4639-b469-d3e721a6baa3
md"""
###
"""
# ╔═╡ 078a97c5-d0f9-4db7-bbe5-127b54c12e2c
md"""
## Setting Up The Problem
As described in the introduction, the Poisson-Nernst-Planck (PNP) equations describe the motions of charged particles that are encountered in various biological, chemical, and physical processes. In a nondimensionalized and simple form, the PNP equations can be written as the following system of coupled nonlinear partial differential equations:
$\begin{align*}
−∇ · ∇φ = q = \sum_{i=1}^{N} z_ic_i \\
\vec J_i = −c_i(∇µ_i + z_i∇φ) \quad (i = 1 . . . N) \\
∂_tc_i + ∇ · \vec J_i = 0 \quad (i = 1 . . . N) \\
\end{align*}$
where the first equation describes the electrostatic potential via the Poisson equation, the second equation describes the flux term, and the third equation is the continuity equation. The second and third equation combined form the Nernst-Planck equation. Altogether we get the Poisson-Nernst-Planck equations. Here, $φ$ is the electrostatic potential, $c_i$ are the species concentrations, and $µ_i$ are the species chemical potentials—the energy necessary to keep a local species configuration. $\vec J_i$ are the species fluxes, $z_i$ the species charge numbers, $q$ is the charge density. The system describes the self-consistent electric field maintained by the charged species, and the species motion due to the electric field $∇φ$ and the gradient of the chemical potential $∇µ_i$. Furthermore, we require that the system be (thermodynamically) closed by establishing a relationship between chemical potential and species concentration. This physical constraint can be codified by the Boltzmann approximation:
$\begin{aligned}
µ_i = log \ c_i \quad (i = 1 . . . N) \\
\end{aligned}$
In electrolytes, a moving species $c_i$ needs to displace the solvent with concentration $c_0$ (and chemical potential µ_0), leading to a modification of the system:
$\begin{aligned}
\vec J_i = −D_ic_i(∇µ_i − ∇µ_0 + z_i∇φ) \quad (i = 1 . . . N)
\end{aligned}$
Once again, we can set $µ_i = log \ c_i$ with $(i = 0 . . . N)$ and add the incompressibility relationship $c_0 + c_1 + · · · + c_N = 1$ that ensures that no concentration can exceed 1. This is physically reasonable as we will have to assume that ions have finite sizes. This relationship allows to express $c_0$ via $c_1 . . . c_N$ and thus to keep just $N$ flux equations. Here the $D_i$ is the diffusion coefficient for each species. The determination of diffusion coefficient for each chemical species would require a combination of extra physical parameters to determine and adjust in order to define the diffusion coefficient. The diffusion coefficient is the amount of a particular substance that diffuses across a unit area in $1$ time step (units usually $1$ second) under the influence of a gradient of one unit. Since the problem is nondimensionalized, we can use the nondimensionalized form of this definition (akin to the notion of "natural units" in some mathematical physics problems) and say that all of the chemical species in our problem have this coefficient set to equal $1$.
For the project we discuss the problem in an $1D$ setting in a domain $Ω = [0, L]$ with $L = 20$. At $x = 0$, we pose the following boundary conditions:
$\begin{align*}
φ = Φ \\
\vec J_i· n = 0 \quad (i = 1 . . . N)
\end{align*}$
For $x = L$, we pose the following boundary conditions:
$\begin{align*}
φ = 0 \\
c_i = c^∗ \quad (i = 1 . . . N)
\end{align*}$
At $t = 0$ we set the initial conditions $c_i = c^∗$ for $(i = 1 . . . N)$. We look for the simulation in the interval $[0, T]$ with $T = 10^3$ time steps (clarification provided below in "Generate the grid" and "Set values for time evolution"). We assume $N = 2$ and $z_1 = 1$, $z_2 = −1$. Note that for the both the boundary conditions and the initial value, the charge density, $q$, is 0. This situation corresponds to the charging of an ideally polarized electrode, where the positive applied potential
attracts negative ions and repels positive ones, thus creating a double layer with nonzero charge density. With the problem stated, we will now proceed to discuss the finite volume space discretization approach we used for the numerical simulation.
"""
# ╔═╡ 6615b610-6157-4342-b9a0-519e44db5410
md"""
#### Generate the grid
Specifications:
- 1D
- Domain: $\Omega = [0,20]$
- stepsize = $0.2$
"""
# ╔═╡ acfaaff3-2a03-4e87-a093-c21403b23af3
# ╠═╡ show_logs = false
begin
# Set domain begin and end
domain_0 = 0
domain_L = 20
grid_step_size = 100
# Create grid using ExtendableGrids
grid=ExtendableGrids.simplexgrid(collect(range(domain_0,domain_L,length=grid_step_size+1)))
# Visulaize grid
GridVisualize.gridplot(grid,resolution=(600,200),Plotter=PyPlot)
end
# ╔═╡ dd7693a1-0e3d-4b82-8502-fb5b55169615
md"""
#### Set values for time evolution
"""
# ╔═╡ a73b7735-be1c-452c-924a-aa17a3167462
begin
ini_time = 0
fin_time = 1000
dt_step = 1
md"""
Specifications:
- The initial time: ini_time $= 0$
- The final time: fin_time $= 1000$
- The time step increment: dt_step $= 1$
"""
end
# ╔═╡ 2d219416-463d-4def-910b-22a06671f86a
# The Vornoi finite volume method: https://j-fu.github.io/VoronoiFVM.jl/stable/method/
md"""
## Finite Volume Space Discretization of the 1D PNP Equations
The finite volume method has two main ingredients: a geometry based approach to obtain a system describing communicating control volumes and a consistent description of the fluxes between two adjacent control volumes. Rather than reinventing the wheel, we will quote the [VoronoiVFM documentation](https://j-fu.github.io/VoronoiFVM.jl/stable/method/#The-discretization-approach) (which we relied upon heavily) to give a portion of the overview of the finite volume space discretization approach we used. We will further elaborate upon this brief overview. The whole finite volume space discretization of the problem being simulated is thoroughly written below the overview and elaboration.

"Given a continuity equation $∇⋅ \vec j=f$ in a domain $Ω$, integrate it over a control volume $\omega_k$ with associated node/collocation point $\vec x_k$, apply Gauss' Divergence Theorem, along with Newton-Leibniz Formula/The Fundamental Theorem of Calculus:
$\begin{aligned}
0&=\int_{\omega_k} (\nabla\cdot \vec j -f )\ d\omega
=\int_{\partial\omega_k} \vec j\cdot \vec n ds - \int_{\omega_k} f d\omega\\
&=\sum_{l\in N_k} \int_{\omega_k\cap \omega_l} \vec j\cdot \vec n ds + \int_{\partial\omega_k\cap \partial\Omega} \vec j\cdot \vec n ds - \int_{\omega_k} f d\omega \\
&\approx \sum_{l\in N_k} \frac{\sigma_{kl}}{h_{kl}}g(u_k, u_l) - |\omega_k| f_k + \text{boundary terms}
\end{aligned}$
Here, $N_k$ is the set of neighbor control volumes, $\sigma_{kl} = |\omega_k\cap \omega_l|$, $h_{kl}=|\vec x_k - \vec x_l|$, where $|⋅|$ denotes the measure (length resp. area) of a geometrical entity. In the approximation step, we replaced the normal flux integral over the interface between two control volumes by the measure of this interface multiplied by a function depending on the unknowns $u_k$, $u_l$ associated to the respective nodes divided by the distance between these nodes. The flux function can be derived from usual finite difference formulas discretizing a particular flux law."
For the boundary terms to be satisfied, we have to implement the various boundary condition types. "To implement a Robin boundary condition on $Γ=∂Ω$:
$- \vec j \cdot \vec n + a u = b,$
we note that by the very construction, the discretization nodes associated to control volumes adjacent to the domain boundary are located at the domain boundary, thus we can assume that the boundary condition is valid in the corresponding collocation node $u_k$. We assume that $\partial\omega_k\cap \partial_\Omega= \cup_{m\in\mathcal M_k} \gamma_{km}$ is the union of a finite number of line (plane) segments. For interior nodes, we set $\mathcal M_k = \emptyset$. Thus, for the boundary terms in the above equation, we have
$\begin{aligned}
\text{boundary terms}&=\sum_{m\in\mathcal M_k} \int_{\gamma_{km}} \vec j \cdot \vec n d s
&\approx \sum_{m\in\mathcal M_k} |\gamma_{km}| \vec j \cdot \vec n\\
&\approx\sum_{m\in\mathcal M_k} |\gamma_{km}| (au_k -b),
\end{aligned}$
We observe that for $\varepsilon\to 0$, the Robin boundary condition
$- \vec j \cdot \vec n + \frac{1}{\varepsilon}u = \frac{1}{\varepsilon}g$
tends to the Dirichlet bundary condition
$u=g.$
Therefore, a Dirichlet boundary condition can be approximated by choosing a small value of $\varepsilon$ and implying the aforementioned Robin boundary conditions. This approach called _penalty method_ is chosen for the implementation of Dirichlet boundary conditions in this package.
The entities describing the discrete system can be subdivided into two categories: Geometrical data and Physics-based data. The Geometrical data: $|\omega_k|$, $\gamma_k$, $\sigma_{kl}$, $h_{kl}$ together with the connectivity information simplex grid. These data are calculated from the discretization grid. For the Physics-based data, we have that the number of species and the functions $s$, $g$, $r$, $f$ etc. describing the particular problem. The solution of the nonlinear systems of equations is performed by Newton's method combined with various direct and iterative linear solvers. The Jacobi matrices used in Newton's method are assembled from the constitutive functions whith the help of forward mode automatic differentiation implemented in [ForwardDiff.jl](https://github.com/JuliaDiff/ForwardDiff.jl)" (which uses automatic differentiation via dual numbers and the chain rule).
Now to further elaborate: we want to find some unknowns $u(x,t)$ our solutions that are local species that are stored in $s$, with flux $\vec j$ discretized as $g$, and some reaction term $r$. When we take this type of equations, we can look at them in some representative elementary volume (REV), and apply Gauss' (Ostrogradsky's) Divergence theorem and Newton-Leibniz Rule/Fundamental Theorem of Calculus, in order to make some kind of species balance in this REV. This then gives the change of the amount of species between two moments of time steps that is proportional to what is flowing in or out and to what is created or destroyed during the reaction. This type of view on the problem is very closely linked to the basic physical principles that are used in order to derive these types of systems of partial differential equations (continuity equations/conservational law formulations). Thus, we can discretize such a system by subdividing our domain into artifical REVs (the $\omega_k$ above) with collocation points $\vec x_k$, where we can assign a value of our solution-unknown variable $u$ to certain collocation points associated to each of these REVs. Then we can watch the evolution of these $u$ evaluated at the collocation points due to the integral-balance. Essentially, we get a system of nonlinear equations on the neighborhood graph of the REV subdivision. Thus, the change of the amount of species during a time interval is proportional to the sum of the fluxes between neighboring REVs and to the amount of species created or destroyed due to the reaction. The time evolution is discretized using implicit Euler for stability purposes. An important point to keep in mind is that the integral formulation of the finite volume method requires that the fluxes are defined as the normal fluxes inside or outside these REVs, along with having the subdivision being formulated in such a way that the interface between two REVs is orthogonal to the two collocation points associated to these two neighboring REVs (also called control volumes). The REV subdivision is generated by first creating boundary conforming Delaunay triangulation of the domain $\Omega$ and then creating a dual Voronoi tessellation where the REVs are now the Voronoi cells.
The finite volume space discretization method has some theoretical advantages: the way this method is derived and formulated essentially exhibits local mass conservation (an important physical property of these systems like the $1D$ PNP); then for convection-diffusion problems, there are robust concepts of upwinding (for discretization of the flux) that allow to avoid, in a guaranteed way (via discrete maximum principle and the M-property of the discretization matrix), unphysical oscillations and we can guarantee positive concentrations in cases where this is important—for the problem of PNP discretization; it can be demonstrated in many cases (in particular for the PNP problem being examined in this project) that there is consistency to thermodynamic principles (a discrete second law of thermodynamics); and finally there are convergence theories for nonlinear systems based on compactness methods.
The flux used in this project is the [SEDAN flux](https://www.wias-berlin.de/preprint/2811/wias_preprints_2811.pdf), which is named after the SEDAN III semiconductor device code. This flux is based on using a discretization approach that is called the Scharfetter-Gummel scheme to approximate the flux at the control volume interface (this is a two-point flux approximation finite volume scheme). This is an exponential fitting upwind flux scheme, which means we have guaranteed sign pattern and M-property. There exist two other variants of fluxes based upon extensions of the Scharfetter-Gummel scheme: the [activity-based flux](https://hal.archives-ouvertes.fr/hal-02194604v3/document) and the [The Bessemoulin-Chatard flux](https://hal.archives-ouvertes.fr/hal-02194604v3/document). Another possible flux discretization approach is the [centered flux](https://www.wias-berlin.de/preprint/2811/wias_preprints_2811.pdf)/[central finite differences](https://www.wias-berlin.de/preprint/2263/wias_preprints_2263_20160708.pdf). This particular approach of using central finite differences for the flux is not used since the scheme may become unstable resulting under certain conditions to large oscillations. Furthermore, the maximum principle can be violated and boundary conditions may be unphysical.
Aside from other ways of defining the discretized flux for the project, it is also possible to use two other discretization solution methods: finite differences and finite elements. The comparison of these three methods is adequately summarized in [Finite Volume Methods by Robert Eymard, Thierry Gallouët, Raphaèle Herbin](https://hal.archives-ouvertes.fr/hal-02100732v2/document), "From the industrial point of view, the finite volume method is known as a robust and cheap method
for the discretization of conservation laws (by robust, we mean a scheme which behaves well even for particularly difficult equations, such as nonlinear systems of hyperbolic equations and which can easily be extended to more realistic and physical contexts than the classical academic problems). The finite volume method is cheap thanks to short and reliable computational coding for complex problems. It may be more
adequate than the finite difference method (which in particular requires a simple geometry). However, in some cases, it is difficult to design schemes which give enough precision. Indeed, the finite element
method can be much more precise than the finite volume method when using higher order polynomials, but it requires an adequate functional framework which is not always available in industrial problems. Other more precise methods are, for instance, particle methods or spectral methods but these methods can be more expensive and less robust than the finite volume method."
As stated at the beginning of this section on describing the finite volume space discretization approach in general, below is the full finite volume space discretization approach to our problem. We then used VoronoiVFM to implement the discretization and run the numerical simulations.
$\begin{align*}
\\
\\
\\
\partial _x \partial _x \phi - z_1 c_1 - z_2 c_2 &=& 0 \\
1-v_0c_0-v_1c_1-v_2c_2 &=& 0 \\
\partial _t c_i +\partial _x [-D_ic_i(\partial _x \log \frac{c_i}{\overline{c}} - \partial _x k_i \frac{c_0}{\overline{c}} + z_i \partial _x \phi)] &=& 0 \\
\\
\end{align*}$
Flux terms:
$\begin{align*}
J_{\phi} &=& \partial _x \phi \\
J_{c_0} &=& 0 \\
J_{c_1} &=& -D_1c_1(\partial _x \log \frac{c_1}{\overline{c}} - \partial _x k_1 \frac{c_0}{\overline{c}} + z_1 \partial _x \phi) \\
J_{c_2} &=& -D_2c_2(\partial _x \log \frac{c_2}{\overline{c}} - \partial _x k_2 \frac{c_0}{\overline{c}} + z_2 \partial _x \phi) \\
\end{align*}$
Reaction terms:
$\begin{align*}
R_{\phi} &=& - z_1 c_1 - z_2 c_2 \\
R_{c_0} &=& 1-v_0c_0-v_1c_1-v_2c_2 \\
R_{c_1} &=& 0 \\
R_{c_2} &=& 0
\end{align*}$
---
From $J_{c_{1/2}}$ to $g_{c_{1/2}}$:
$J_{c_{1/2}} = -D_ic_i(\partial _x \log \frac{c_i}{\overline{c}} - \partial _x k_i \frac{c_0}{\overline{c}} + z_i \partial _x \phi)$
Desired form:
$J_{des} = -D_i(\partial _x c_i - c_i \frac{v}{D_i})$
Since $\partial _x \log c_i = \frac{1}{c_i} \cdot \partial _x c_i$:
$J_{c_{1/2}} = -D_i [\partial _x c_i - c_i (k_i \partial _x \log(c_0) - (k_i-1)\log(\overline{c}) - z_i \partial _x \phi)]$
$v = D_i [k_i \partial _x \log(c_0) - (k_i-1)\log(\overline{c}) - z_i \partial _x \phi]$
$v_{kl} = \frac{D_i}{|\sigma|} \cdot[k_i \int _{\sigma} \partial _x \log(c_0) d \gamma - (k_i-1) \int _{\sigma} \partial _x \log(\overline{c}) d \gamma - z_i \int _{\sigma} \partial _x \phi d \gamma]$
$v_{kl} = \frac{D_i}{h_{kl}} [k_i \log(\frac{c_{0_R}}{c_{0_L}}) + (k_i-1) \log(\frac{\overline{c}_L}{\overline{c}_R})+ z_i (\phi _L - \phi _R)]$
$g_{c_{1/2}} = D_i(B(\frac{v_{kl}h_{kl}}{D_i})c_L - B(\frac{-v_{kl}h_{kl}}{D_i})c_R)$
---
Applied BC
$\phi (0) = \frac{1}{2}$
$\phi (L) = 0$
$c_i(L) = \frac{1}{3}$
$J_i \cdot n = 0$
---
Applied IC
$\phi = \frac{1}{2}$
$c_i = \frac{1}{3}$
---
Defined Variables
$D_i = 1$
$z_0 = 0$
$z_1 = -1$
$z_2 = 1$
$v_0 = 1$
$v_1 = 1$
$v_2 = 1$
---
Other Variables (for clearification)
$\begin{align*}
\overline{c} &=& \sum c_i \\
k_i &=& \frac{v_i}{v_0} = const. \\
\end{align*}$
$\newline$
"""
# ╔═╡ 0e996cc2-ce10-44ce-8698-bc77ee0819e4
md"""
## VoronoiFVM Implementation of FVM Discretization of 1D PNP Equations
"""
# ╔═╡ 2c0b8207-5f12-424d-81b8-d4adb2c202c7
begin
# Create the mutable struct "Data" which is able to hold all the needed variables
mutable struct Data
z_1::Float64
z_2::Float64
D::Float64
v_0::Float64
v_1::Float64
v_2::Float64
k_1::Float64
k_2::Float64
ic0::Int32
ic1::Int32
ic2::Int32
iphi::Int32
Data()=new()
end
# Set which species are under time derivate
function storage!(f,u,node,data)
# Import values from data
ic0 = data.ic0
ic1 = data.ic1
ic2 = data.ic2
iphi = data.iphi
# c_1 and c_2 are time derivable ones
f[iphi] = 0
f[ic0] = 0
f[ic1] = u[ic1]
f[ic2] = u[ic2]
end
# Set reaction terms
function reaction!(f,u,node,data)
# Import values from data
ic0 = data.ic0
ic1 = data.ic1
ic2 = data.ic2
iphi = data.iphi
z_1 = data.z_1
z_2 = data.z_2
v_0 = data.v_0
v_1 = data.v_1
v_2 = data.v_2
# Reaction terms
f[iphi] = -(z_1*u[ic1]+z_2*u[ic2])
f[ic0] = 1 - v_0*u[ic0] - v_1*u[ic1] - v_2*u[ic2]
f[ic1] = 0
f[ic2] = 0
end
# Set flux terms
function sedanflux!(f,u0,edge,data)
u=unknowns(edge,u0)
# Import values from data
ic0 = data.ic0
ic1 = data.ic1
ic2 = data.ic2
iphi = data.iphi
D = data.D
z_1 = data.z_1
z_2 = data.z_2
k_1 = data.k_1
k_2 = data.k_2
# Calculate the logarithmns needed for the flux terms
log_min = -5
log_res_c0 = log(abs(u[ic0,2]/u[ic0,1]))
log_res_c012 = log(abs((u[ic0,1]+u[ic1,1]+u[ic2,1])/(u[ic0,2]+u[ic1,2]+u[ic2,2])))
# If the logarithms smaller than log_min, set the logarithms to log_min to avoid errors caused by high negativ logarithmns
if log_res_c0 < log_min
log_res_c0 = log_min
end
if log_res_c012 < log_min
log_res_c012 = log_min
end
# Flux terms
f[iphi] = u[iphi,1]-u[iphi,2]
f[ic0] = 0
bp1, bm1 = fbernoulli_pm(k_1*log_res_c0+(k_1-1)*log_res_c012+z_1*(u[iphi,1]-u[iphi,2])) # bernoulli function
bp2, bm2 = fbernoulli_pm(k_2*log_res_c0+(k_2-1)*log_res_c012+z_2*(u[iphi,1]-u[iphi,2])) # bernoulli function
f[ic1] = D*(bm1*u[ic1, 1]-bp1*u[ic1, 2])
f[ic2] = D*(bm2*u[ic2, 1]-bp2*u[ic2, 2])
end
# Main function to solve the problem
function main_problem(;n=20,Plotter=nothing,dlcap=false,verbose=false,unknown_storage=:sparse,DiffEq=nothing, t_min, t_max, t_step, t_step_max, c1_dir, c2_dir, v_0, v_1, v_2)
# Specify h and grid
h=20/convert(Float64,n)
grid=VoronoiFVM.Grid(collect(0:h:20))
# Save needed variables in data
data = Data()
data.z_1 = -1
data.z_2 = 1
data.D = 1
data.v_0 = v_0
data.v_1 = v_1
data.v_2 = v_2
## Species number
data.iphi = 1
data.ic1 = 2
data.ic2 = 3
data.ic0 = 4
# Import values from data
ic1 = data.ic1
ic2 = data.ic2
iphi= data.iphi
ic0 = data.ic0
v_0 = data.v_0
v_1 = data.v_1
v_2 = data.v_2
# Calculate k_1 and k_2 and save it in data
data.k_1 = v_1 / v_0
data.k_2 = v_2 / v_0
# Set up the physics
physics=VoronoiFVM.Physics(data=data,
num_species=4,
flux=sedanflux!,
reaction=reaction!,
storage=storage!
)
# Create finite volume system
sys=VoronoiFVM.System(grid,physics,unknown_storage=unknown_storage)
# Enable species
enable_species!(sys,1,[1])
enable_species!(sys,2,[1])
enable_species!(sys,3,[1])
enable_species!(sys,4,[1])
# Set boundary conditions
## Dirichlet
boundary_dirichlet!(sys,iphi,1,0.5)
boundary_dirichlet!(sys,iphi,2,0.0)
boundary_dirichlet!(sys,ic0,2,1-c1_dir-c2_dir)
boundary_dirichlet!(sys,ic1,2,c1_dir)
boundary_dirichlet!(sys,ic2,2,c2_dir)
# Neumann
boundary_neumann!(sys,ic0,1,0)
boundary_neumann!(sys,ic1,1,0)
boundary_neumann!(sys,ic2,1,0)
# Set initial conditions
inival=unknowns(sys)
@views inival[iphi,:].=0.5
@views inival[ic0,:].=1-c1_dir-c2_dir
@views inival[ic1,:].=c1_dir
@views inival[ic2,:].=c2_dir
# Solve the system
## Create solver control info for constant time step size
control=VoronoiFVM.NewtonControl()
control.verbose=verbose
control.Δt_min=0
control.Δt=t_step
control.Δt_grow=1.2
control.Δt_max=t_step_max
control.Δu_opt=100
control.damp_initial=0.5
if isnothing(DiffEq)
tsol=solve(inival,sys,[t_min,t_max],control=control)
else
tsol=solve(DiffEq,inival,sys,[0.0,10],
initializealg=DiffEq.NoInit(),
dt=t_step)
end
# Return solution
return tsol
end
end
# ╔═╡ 0d1ac969-ce6a-45c6-9c7e-3675cc858055
md"""
### Calculate Solution
"""
# ╔═╡ cb94ba74-a654-4ccd-bfd0-7309e3899c21
# ╠═╡ show_logs = false
solution = main_problem(;n=grid_step_size,Plotter=PyPlot,dlcap=true,verbose=true,unknown_storage=:sparse,DiffEq=nothing, t_min=ini_time, t_max=fin_time, t_step=dt_step, t_step_max=0.1, c1_dir=1/3, c2_dir=1/3, v_0=1, v_1=1, v_2=1);
# ╔═╡ fdf88786-5691-422c-9484-fa905c47bcde
md"""
#
"""
# ╔═╡ d4ea8aa9-cec8-4cfa-8ed3-22a32f0c0fdc
md"""
## Presentation of Results
"""
# ╔═╡ 95ecd8a5-8d5f-48ef-a4c4-b4eb1ac131c5
# Define some functions to make plotting easier
begin
# Create visualizer and set x limit
visualizer=GridVisualizer(Plotter=PyPlot, resolution=(700,700),layout=(2,1), xlimits=(domain_0, domain_L))
# Function for plotting our scalar plots
function scalarPlots(solution, num_species, time_var)
t = round(solution.t[time_var]*100)/100
scalarplot!(visualizer[1,1],grid,solution[:,time_var][1,:],title="Potential at time t=$t",flimits=(minimum(solution),maximum(solution)*1.05), xlabel = latexstring("Domain \$\\Omega\$"), ylabel = latexstring("Potential \$\\Phi\$"))
scalarplot!(visualizer[2,1],grid,solution[:,time_var][2,:],title="Concentrations at time t=$t",flimits=(0,1),color=:blue,label="\$c_1\$", xlabel = latexstring("Domain \$\\Omega\$"), ylabel = latexstring("Concentration \$c\$"), legend=:best)
scalarplot!(visualizer[2,1],grid,solution[:,time_var][3,:],flimits=(0,1),color=:red,clear=false,label="\$c_2\$", xlabel = latexstring("Domain \$\\Omega\$"), ylabel = latexstring("Concentration \$c\$"))
if num_species == 4
scalarplot!(visualizer[2,1],grid,solution[:,time_var][4,:],flimits=(0,1),color=:green,clear=false,label="\$c_0\$", xlabel = latexstring("Domain \$\\Omega\$"), ylabel = latexstring("Concentration \$c\$"))
end
reveal(visualizer)
end
# Function to get the arrays needed for our contour plots
function contourArrays(input_array, num_species, grid_len)
time_len = length(input_array)
output_dict = Dict()
for i in 1:num_species
temp_array = input_array[:,1][i,:] # temporary array for cache
for j in 2:time_len
append!(temp_array, input_array[:,j][i,:])
end
output_dict[i] = temp_array
end
return output_dict, time_len
end
# Function for plotting the countor plots
function contourPlots(contour_dict, num_species, time_len, c_max, automatic_lims)
x = grid[Coordinates][:]
y = log10.(1:time_len)
plot_dict = Dict()
if automatic_lims == true
clims_phi = :auto
clims_c = :auto
else
clims_phi = (0,0.5)
clims_c = (0,c_max)
end
plot_dict[1] = Plots.contour(x, y, contour_dict[1], fill = true, title=latexstring("Potential \$\\phi\$"), xlabel = latexstring("Domain \$\\Omega\$"), ylabel = latexstring("timestep in \$10^x\$"), linewidth = 0, c = :redsblues, clims=clims_phi)
plot_dict[2] = Plots.contour(x, y, contour_dict[2], fill = true, title=latexstring("Concentration \$c_1\$"), xlabel = latexstring("Domain \$\\Omega\$"), ylabel = latexstring("timestep in \$10^x\$"), linewidth = 0, c = :roma, clims=clims_c)
plot_dict[3] = Plots.contour(x, y, contour_dict[3], fill = true, title=latexstring("Concentration \$c_2\$"), xlabel = latexstring("Domain \$\\Omega\$"), ylabel = latexstring("timestep in \$10^x\$"), linewidth = 0, c = :roma, clims=clims_c)
if num_species == 4
plot_dict[4] = Plots.contour(x, y, contour_dict[4], fill = true, title=latexstring("Concentration \$c_0\$"), xlabel = latexstring("Domain \$\\Omega\$"), ylabel = latexstring("timestep in \$10^x\$"), linewidth = 0, c = :roma, clims=clims_c)
end
return plot_dict
end
# Function for plotting our animations
function animations(solution, num_species, animation_step)
x = grid[Coordinates][:]
min_sol = round(minimum(solution))
max_sol = round(maximum(solution))
time_len = length(solution)
anim_solution = @animate for i = 1:animation_step:time_len
Plots.contour(x, min_sol:max_sol, append!(solution[:,i][1,:], solution[:,i][1,:]), fill = true, ticks = true, linewidth = 0, c = :redsblues, clims=(0,0.5), title="Potential and Concentrations over time", colorbar_title="Potential \$\\phi\$")
Plots.plot!(x, solution[:,i][2,:], ylims = (min_sol, max_sol), xlims = (domain_0, domain_L), linewidth=3, linecolor="aqua", xlabel = latexstring("Domain \$\\Omega\$"), ylabel = latexstring("Concentration \$c\$"), label="\$c_1\$")
Plots.plot!(x, solution[:,i][3,:], linewidth=3, linecolor="fuchsia", label="\$c_2\$")
if num_species == 4
Plots.plot!(x, solution[:,i][4,:], linewidth=3, linecolor="green", label="\$c_0\$")
end
end
return anim_solution
end
end
# ╔═╡ bfebf9f5-57a1-4d2d-b010-31c33f28c5d5
md"""
### Scalar Plots of Species Concentrations and Potential
#### Slider for time t
"""
# ╔═╡ a2e9f3cb-93bf-4e8d-82a4-bed79423a6ea
@bind time Slider(1:1:length(solution))
# ╔═╡ 23b59e41-c5a8-4de6-882c-019a87d74fef
println("t = "*string(round(solution.t[time]*100)/100))
# ╔═╡ 4ff1af93-59a7-4d0b-b19e-e85d216e4398
scalarPlots(solution, 4, time)
# ╔═╡ cb43d58e-454c-4e65-afe2-f14e876f127a
md"""
"""
# ╔═╡ 6e9cc937-bfde-4aa7-9615-fbb1a45bd0d1
md"""
### Contour Plots of the Time Evolution of Concentrations and Potential
"""
# ╔═╡ a06a6592-4206-41eb-b18f-8a5fa3dd87e8
begin
contour_dict, time_len = contourArrays(solution, 4, grid_step_size + 1)
contour_plots = contourPlots(contour_dict, 4, time_len, 0.5, false)
Plots.plot(contour_plots[2], contour_plots[3])
end
# ╔═╡ b282a67c-799e-406a-97c9-9f8952133c02
Plots.plot(contour_plots[4], contour_plots[1])
# ╔═╡ 914bbe12-14b2-461e-b678-be84e014c40e
md"""
###
"""
# ╔═╡ 69f6fd9b-8a37-46e4-a1cb-94f3d53cc841
md"""
### Animation of the Above Scalar and Contour Plots Combined
"""
# ╔═╡ c4c34aaa-c15a-478e-b86a-0479c3b0dcbf
# ╠═╡ show_logs = false
begin
animated_solution = animations(solution, 4, 10)
gif(animated_solution, "animated_solution.gif", fps = 30)
end
# ╔═╡ 62fdd60b-b26d-4802-91df-aea028296e85
md"""
####
"""
# ╔═╡ 82157790-e629-44b8-864d-09fbb61aa484
md"""
The plots and animations above show that the evolution of potential $\phi$, concentration of negative ($c_1$), positive ($c_2$) and solvent ($c_0$). Both the two charged species and the solvent have equal size/volume $v_0=v_1=v_2=1)$. The results show the following: there is an attraction of the negative charged species ($c_1$) to $x=0$ because $\phi$ at $x=0$ is larger than $\phi$ at $x=L$; and there is a repelling of the positive charged species ($c_2$) from $x=0$ because $\phi$ at $x=0$ is larger than $\phi$ at $x=L$. This leads to the results exhibiting double layer charging. Furthermore, the results show that at $x=L, c_i=c^*=1/3$ for $\ i=0,1,2$, we have that all of the species have the same concentration on the right side of the domain, i.e. $+\phi(x=0)=\phi _0 = 1/2$ and $\phi(x=L) = \phi _L = 0$. This shows that the dirichlet boundary conditions are fullfilled.
It almost makes sense physically—the positive species concentration increases by a marginal amount, and the incompressibility condition though satisfied ends up decreasing the solvent concentration, which affects the mass conservation—solvent concentration should not have decreased. The positive species concentration increased due to the collision flux at the boundary layer where the applied positive potential is located. As soon as the positive chemical species arrives at this boundary, the self-consistent field condition of the Poisson equation of the applied electric potential creates a sort-of-feedback with the positive charges of the chemical species. As stated above, at $x=0$ we have a double layer charging occuring at this boundary where there is an applied electrostatic potential.
This is almost physically reasonable, although not quite correct dynamics of the physical system, but there seems to be an issue with the conservation of concentration of our two species, $c_1$ and $c_2$, and our solvent $c_0$. The concentrations are all positive and conserved. From the graph in our code, it looks like there is a small creation of species concentration for $c_1$ as that particular species diffuses towards our "source" of the initial electric potential on the left-hand boundary at $x = 0$. Because of the self-consistent field approach to our system, the sudden creation of more species concentration around $x = 0$ is due to the reliance of the distribution of charges creating the electric field and the electric field relying upon the distribution of charges—seeing the attraction/repelling of species according to their charge around the boundary layer. This is an almost-physically-correct solution given that there is a slight issue with the rising of concentration values for our species $c_1$ and $c_2$, which is not physically correct but, it can be interpreted as physically correct by observing the self-consistent field approach and the way we set up the "reaction" terms for our discretization within the Voronoi API.
"""
# ╔═╡ 96b27dfd-b4a6-4531-bf55-a4ed717467c2
md"""
###
"""
# ╔═╡ 69df9d73-3733-4b1f-92da-a7cf2622a71e
# ╔═╡ f4675a1a-a6b9-4306-91ba-a0d6710c5afd
md"""
### Plots for Different Ion Volume
$v_1 = 2$
$v_2 = 1$
$v_0 = 1$
"""
# ╔═╡ 1056ab3e-edae-486c-bdcc-84c8c9765abb
# ╠═╡ show_logs = false
solution_diff_vol = main_problem(;n=grid_step_size,Plotter=PyPlot,dlcap=true,verbose=true,unknown_storage=:sparse,DiffEq=nothing, t_min=ini_time, t_max=fin_time, t_step=dt_step, t_step_max=0.1, c1_dir=1/3, c2_dir=1/3, v_0=1, v_1=1.5, v_2=1);
# ╔═╡ 971fc5e4-ce3c-4c39-a7cf-f4e35e2f8c63
begin
contour_dict_diff_vol, time_len_diff_vol = contourArrays(solution_diff_vol, 4, grid_step_size + 1)
contour_plots_diff_vol = contourPlots(contour_dict_diff_vol, 4, time_len_diff_vol, 0.5, false)
Plots.plot(contour_plots_diff_vol[2], contour_plots_diff_vol[3])
end
# ╔═╡ 582814b0-0410-4356-82b3-845fef6e7978
Plots.plot(contour_plots_diff_vol[4], contour_plots_diff_vol[1])
# ╔═╡ 5b283a49-510b-4d9e-b29d-7dd41d0bf3eb
md"""
The plots above involving the finite volume show the evolution of potential $\phi$, the concentrations of the negative ($c_1$), the positive ($c_2$), and the solvent ($c_0$) charged species with volumne of ions for $c_0, c_1,$ and $c_2$ given by $v_0=v_2 = 1$ and $v_1=2$. The volume of the negative charged species is double that of the volume of the positive charged species and the solvent. There is an attraction of the negative charged species ($c_1$) to the $x=0$ boundary due to the positive $\phi$ at $x=0$ being larger than $\phi$ at $x=L$. There is a repelling of the positive charged species ($c_2$) from $x=0$ because $\phi$ at $x=0$ is larger than $\phi$ at $x=L$. This leads to the observation of the double layer charging phenomena. Furthermore, the concentrations of both the negative charged species and the solvent are lower at $x=0$. The concentration of negative charged species is lower in the rest of the domain (but at $x=L$, $c_1=1/3$ for $v_1=1$ and $v_1=2$). The results also show that the concentration of positive charged species is a bit higher over the domain (but at $x=L$, $c_2=1/3$ for $v_1=1$ and $v_1=2$). This is because of similar reasons stated in the previous section on the interpretation and presentation of the plots and animations of the concentrations are all initially equal and the sizes/volumes of species and solvent are all equal.
"""
# ╔═╡ 3ff2cc4c-81a6-4484-a3a0-a0c1916cddf1
md"""
###
"""
# ╔═╡ 6dc1afff-12f4-4f22-a60c-5ba843facadd
# ╔═╡ ef9ce3be-f0bc-4cd9-931e-cfdd0cc559c2
md"""
### Results for Varying Values of $c^*$ for $c_1$ and $c_2$
"""
# ╔═╡ de92b301-eaec-4b10-8a3b-454c6393ef9e
md"""
##### Set values for stepsize, the smallest and the largest possible $c^*$
"""
# ╔═╡ bae36b4f-90d7-41a9-ba30-5607a085fa62
begin
c_start = 0.05
c_step = 0.05
c_end = 0.95
md"""
Specifications:
$c_{start} = 0.05$
$c_{step} = 0.05$
$c_{end} = 0.95$
"""
end
# ╔═╡ a37133fc-0f30-409a-a4aa-57e3de882963
# ╠═╡ show_logs = false
begin
# Going through every combination of c_1 and c_2
begin
global solution_array_c1_c2 = []
for c2_dir in c_start:c_step:c_end # go through every possibility for c_2
global solution_array_c1 = []
for c1_dir in c_start:c_step:c_end # gp through every possibility for c_1
try # try to get a solution for each combination
global solution_combination = main_problem(;n=grid_step_size,Plotter=PyPlot,dlcap=true,verbose=true,unknown_storage=:sparse,DiffEq=nothing, t_min=ini_time, t_max=1000, t_step=1, t_step_max=10, c1_dir, c2_dir, v_0=1, v_1=1, v_2=1)
catch combination_doesnt_work
global solution_combination = [0]
end
push!(solution_array_c1, solution_combination)
end
push!(solution_array_c1_c2, solution_array_c1) # save every combination in solution_array_c1_c2
end
end
md"""
In order to show the results for diffrent combinations of $c_1$ and $c_2$, the results for every combination of $c_1$ and $c_2$ were calculated.
"""
end
# ╔═╡ cee4bd55-1861-4efc-bf04-79521a6803de
md"""
#### Slider for $c_1$ and $c_2$
"""
# ╔═╡ 9848decf-5fca-4fc4-bef4-77463778e6cf
@bind c1_pos Slider(1:1:length(solution_array_c1))
# ╔═╡ 338930f9-0d99-48ec-803a-e4eeb2dac207
c_1 = round((c_start+(c1_pos-1)*c_step)*100)/100
# ╔═╡ d6841195-f4df-46c1-bfa3-2be5c2d54e77
@bind c2_pos Slider(1:1:length(solution_array_c1))
# ╔═╡ b6dc67de-9ba6-4301-93ef-6c88823f5269
c_2 = round((c_start+(c2_pos-1)*c_step)*100)/100
# ╔═╡ f85ed8b7-fede-4458-9e1a-0e2510a932b0
md"""
### Contour plots of the Time Evolution of Varying $c_1$ and $c_2$
"""
# ╔═╡ 25993f86-3457-49b8-bf21-9573a86471d7
begin
try
contour_dict_c1_c2, time_len_c1_c2 = contourArrays(solution_array_c1_c2[c2_pos][c1_pos], 4, grid_step_size + 1)
global contour_plots_c1_c2 = contourPlots(contour_dict_c1_c2, 4, time_len_c1_c2, max(c_1,c_2,1-c_1-c_2), true)
Plots.plot(contour_plots_c1_c2[2], contour_plots_c1_c2[3])
catch
println("Combination does not work")
end
end
# ╔═╡ 228a62af-47c3-485e-adde-d59e980e426f
begin
try
Plots.plot(contour_plots_c1_c2[4], contour_plots_c1_c2[1])
catch
println("Combination does not work")
end
end
# ╔═╡ e9198ceb-479b-4b95-93a8-6f3c1beb3c43
md"""
The above plots show the following:
When reducing the concentration of the positive charged species $c_2$ (in comparison to the case in the previous section above), it is observed that at end of time evolution $\phi$ at $x$ is larger than $\phi$ at $x=L$, which is $0$ for most part of the interior of the domain. The concentration of the positive charged species is low on both sides of the domain than in the case above in the previous section. The concentration of the negative charged species is higher at $c_1$ at $x=0$ than at $c_1$ at $x=L$. The concentration of the solvent is, in general (across the domain), higher throughout the evolution.
When increasing the concentration of the positive charged species $c_2$ (in comparison to the case in the previous section above), it is observed that at the end of the time evolution $\phi$ at $x$ is larger than $\phi$ at $x=L$, which is $0$ for all $x \in \Omega$. Results also show that $c_2$ at $x=0$ is smaller than $c_2$ at $x=L$ and $c_1$ at $x=0$ is larger than $c_1$ at $x=L$. The concentration of the solvent is, in general (across the domain), lower throughout the evolution.
When reducing the concentration of the negative charged species $c_1$ (in comparison to the case above), it is observed that at the end of the time evolution $\phi$ at $x$ is larger than $\phi$ at $x=L$, which is equal to $0$ for all $x \in \Omega$. Furthermore, the following is observed: $c_2$ at $x=0$ is smaller than $c_2$ at $x=L$; $c_1$ at $x=0$ is larger than $c_1$ at $x=L$; and the concentration of the solvent is, in general (across the domain), higher throughout the evolution.
When increasing the concentration of the negative charged species $c_1$ (in comparison to the case above), it can be seen that at the end of the time evolution
$\phi$ at $x$ is larger than $\phi$ at $x=L$, which is equal to $0$ for most of the interior of the domain. Furthermore, the following observations can be made: $c_2$ at $x=0$ is smaller than $c_2$ at $x=L$; the concentration of the negative charged species is higher on both the domain and the boundaries relative to the case above; $c_1$ at $x=0$ is larger than $c_1$ at $x=L$; and finally the concentration of the solvent is, in general (across the domain), lower throughout the evolution
Finally, it is observed that some combinations of concentrations are not feasible. This could be due to the combination of the incompressibility condition being imposed and the issues discussed in the first results section documenting the issues with concentration conservation and the applied electrostatic potential building a double layer that may be affecting the species concentration due to the self-consistent field.
"""
# ╔═╡ a3192f15-7ac5-485e-9a32-556226f1faa0
# ╔═╡ 6d70eea9-7ff9-4803-8c04-bf8b7b8301cc
md"""
# Simplification of Original Problem by Removal of Solvent $c_0$
"""
# ╔═╡ 070a943f-3e50-4c9b-bd55-d91ae124f85a
md"""
## Finite Volume Space Discretization of the Simplied 1D PNP Equations (removal of $c_0$)
We remove the solvent from the original 1D PNP equations in order to compare this simplification to the original 1D PNP equations above. By removing the solvent term, we have a simplified two chemical species system with the corresponding incompressibility condition and volume condition. The reaction and flux terms are defined below.
$\begin{align*}
\partial _x \partial _x \phi - z_1 c_1 - z_2 c_2 &=& 0 \\
1-v_1c_1-v_2c_2 &=& 0 \\
\partial _t c_i +\partial _x [-D_ic_i(\partial _x \log \frac{c_i}{\overline{c}} + z_i \partial _x \phi)] &=& 0 \\
\\
\end{align*}$
Flux terms:
$\begin{align*}
J_{\phi} &=& \partial _x \phi \\
J_{c_1} &=& -D_1c_1(\partial _x \log \frac{c_1}{\overline{c}} + z_1 \partial _x \phi) \\
J_{c_2} &=& 0 \\
\end{align*}$
Reaction terms:
$\begin{align*}
R_{\phi} &=& - z_1 c_1 - z_2 c_2 \\
R_{c_1} &=& 0 \\
R_{c_2} &=& 1-v_1c_1-v_2c_2 \\
\end{align*}$
---
From $J_{c_{1/2}}$ to $g_{c_{1/2}}$:
$J_{c_{1}} = -D_1c_1(\partial _x \log \frac{c_1}{\overline{c}} + z_1 \partial _x \phi)$
Desired form:
$J_{des} = -D_1(\partial _x c_1 - c_1 \frac{v}{D_1})$
Since $\partial _x \log c_1 = \frac{1}{c_1} \cdot \partial _x c_1$:
$J_{c_{1}} = -D_1 [\partial _x c_1 - c_1 (\log(\overline{c}) - z_1 \partial _x \phi)]$
$v = D_1 [\log(\overline{c}) - z_1 \partial _x \phi]$
$v_{kl} = \frac{D_1}{|\sigma|} \cdot[\int _{\sigma} \partial _x \log(\overline{c}) d \gamma - z_1 \int _{\sigma} \partial _x \phi d \gamma]$
$v_{kl} = \frac{D_1}{h_{kl}} [\log(\frac{\overline{c}_R}{\overline{c}_L})+ z_1 (\phi _L - \phi _R)]$
$g_{c_{1}} = D_i(B(\frac{v_{kl}h_{kl}}{D_1})c_L - B(\frac{-v_{kl}h_{kl}}{D_1})c_R)$
---
Applied BC
$\phi (0) = \frac{1}{2}$
$\phi (L) = 0$
$c_i(L) = \frac{1}{2}$
$J_i \cdot n = 0$
---
Applied IC
$\phi = \frac{1}{2}$
$c_i = \frac{1}{2}$
---
Defined Variables
$D_i = 1$
$z_1 = -1$
$z_2 = 1$
$v_1 = 1$
$v_2 = 1$
---
Other Variables (for clearification)
$\begin{align*}
\overline{c} &=& \sum c_i = c_1 + c_2 \\
\end{align*}$
$\newline$
"""
# ╔═╡ 675092ef-b725-4c00-ba0d-26b82ac4735c
md"""
## VoronoiFVM Implementation of FVM Discretization of the Simplified 1D PNP Equations
"""
# ╔═╡ 3aa234a5-2058-4715-af9b-f2a1d4bca2b4
# Works like above - but diffrent reaction & flux terms
begin
function storage_reduced!(f,u,node,data)
ic1 = data.ic1
ic2 = data.ic2
iphi = data.iphi
f[iphi] = 0
f[ic1] = u[ic1]
f[ic2] = 0
end
function reaction_reduced!(f,u,node,data)
ic1 = data.ic1
ic2 = data.ic2
iphi = data.iphi
z_1 = data.z_1
z_2 = data.z_2
v_1 = data.v_1
v_2 = data.v_2
f[iphi] = -(z_1*u[ic1]+z_2*u[ic2])
f[ic1] = 0
f[ic2] = 1 - v_1 * u[ic1] - v_2 * u[ic2]
end
function sedanflux_reduced!(f,u0,edge,data)
u=unknowns(edge,u0)
ic1 = data.ic1
ic2 = data.ic2
iphi = data.iphi
D = data.D
z_1 = data.z_1
z_2 = data.z_2
f[iphi] = u[iphi,1]-u[iphi,2]
bp1, bm1 = fbernoulli_pm(log(abs((u[ic1,2]+u[ic2,2])/(u[ic1,1]+u[ic2,1])))+z_1*(u[iphi,1]-u[iphi,2]))
bp2, bm2 = fbernoulli_pm(log(abs((u[ic1,2]+u[ic2,2])/(u[ic1,1]+u[ic2,1])))+z_2*(u[iphi,1]-u[iphi,2]))
f[ic1] = D*(bm1*u[ic1, 1]-bp1*u[ic1, 2])
f[ic2] = 0
end
function main_problem_reduced(;n=20,Plotter=nothing,dlcap=false,verbose=false,unknown_storage=:sparse,DiffEq=nothing, t_min, t_max, t_step)
h=20/convert(Float64,n)
grid=VoronoiFVM.Grid(collect(0:h:20))
data = Data()
data.z_1 = -1
data.z_2 = 1
data.D = 1
data.v_1 = 1
data.v_2 = 1
data.iphi = 1
data.ic1 = 2
data.ic2 = 3
ic1 = data.ic1
ic2 = data.ic2
iphi= data.iphi
v_1 = data.v_1
v_2 = data.v_2
physics=VoronoiFVM.Physics(data=data,
num_species=3,
flux=sedanflux_reduced!,
reaction=reaction_reduced!,
storage=storage_reduced!
)
sys=VoronoiFVM.System(grid,physics,unknown_storage=unknown_storage)
enable_species!(sys,1,[1])
enable_species!(sys,2,[1])
enable_species!(sys,3,[1])
boundary_dirichlet!(sys,iphi,1,0.5)
boundary_dirichlet!(sys,iphi,2,0.0)
boundary_dirichlet!(sys,ic1,2,1/2)
boundary_dirichlet!(sys,ic2,2,1/2)
boundary_neumann!(sys,ic1,1,0)
boundary_neumann!(sys,ic2,1,0)
inival=unknowns(sys)
@views inival[iphi,:].=0.5
@views inival[ic1,:].=1/2
@views inival[ic2,:].=1/2
tstep=t_step
control=VoronoiFVM.NewtonControl()
control.verbose=verbose
control.Δt_min=0
control.Δt=tstep
control.Δt_grow=1.2