PDF
Capacitive Deionization of Saline Water
Introduction
Desalination is the removal of salt (NaCl) from, for example, brine or seawater to obtain water suitable for human use. Efficient production of potable water from seawater is useful in parts of the world where freshwater access is limited. One method for desalination is capacitive deionization, where the sodium and chloride ions in seawater are driven into porous electrodes by an external applied electrical potential. This model shows how to implement the improved modified Donnan model for capacitive deionization in a “flow-between” capacitive deionization cell.
Model Definition
Geometry
The geometry of the system is shown in Figure 1.
Figure 1: Geometry of the capacitive deionization system. For clarity, the y-axis has been scaled by a factor of 30.
The inlet is located at the left-most vertical boundary and the outlet at the corresponding right-hand side vertical boundary. The main flow channel, from the inlet to the outlet, is bracketed by the two horizontal porous electrodes. The upper, positive, porous electrode is connected to the external circuit along its upper boundary, and the lower, negative, porous electrode is connected to the external circuit along its lower boundary.
Fluid Flow
The flow through the open main channel, and through the porous electrodes, is described using a Brinkman Equations interface. The flow field is solved for independently of the deionization process, and the flow field is then used in a Tertiary Current Distribution, Nernst–Planck interface.
Electrochemistry
The Tertiary Current Distribution, Nernst–Planck interface solves for the electrode and electrolyte electric potentials, as well as for the concentrations of the Na+ and Cl- ions. Electroneutrality is assumed in the electrolyte. A Porous Electrode node is used to describe the transport of species and charge in the porous electrode domains. A more detailed description of the governing equations of this interface can be found in the documentation for the Wire Electrode tutorial.
Improved Modified Donnan Model
In this model, the ions are driven, by the applied potential, and by a nonelectrostatic attractive force, into micropores inside the porous electrodes. This process reduces the salt concentration in the water exiting the cell (deionization). The deionization process is modeled using the “improved modified Donnan” model of Bazant and coworkers (Ref. 1). The improved modified Donnan model extends the Donnan equilibrium model by introducing a concentration- and micropore size-dependent attractive potential between the micropore walls and the ions inside the micropores of the porous electrode. The micropores are assumed to be homogeneously distributed inside the porous electrodes. It is assumed that all capacitive effects are localized to these micropores.
To derive the relationship between the ionic concentrations inside the micropores and the macropores of the porous electrodes, electrochemical equilibrium between the micro- and macropores is assumed:
(1)
In Equation 1, the subscript “i” refers to a species index, “mi” refers to “micropore”, and the subscript “l” refers to the liquid electrolyte. In Equation 1, is the electrochemical potential of species i (J/mol). Applying the definition of the electrochemical potential allows for rewriting Equation 1 in terms of concentrations (ci, mol/m3), electrical potentials (φ, V), and the previously mentioned nonelectrostatic attractive potential (μatt,i, J/mol). For the macropore electrolyte,
(2)
and for the micropore electrolyte,
(3)
In Equation 2 and Equation 3, μi is the chemical potential of species i, zi is the charge of species i, and F is Faraday’s constant (C/mol). The attractive potential is assumed to be the same for all ions, and to be given by
(4)
Here, the attractive volumetric energy density E (J/m3) is given by
(5)
where z is the ionic charge, kB is Boltzmann’s constant (J/K), T is the temperature (K), λB is the Bjerrum length in the micropore (m) and λp is the micropore size (m).
By assuming unit activity (ai, dimensionless), the chemical potential can be rewritten in terms of concentration as
(6)
In Equation 6, R is the molar gas constant (J/K/mol). Substituting these expressions into Equation 1, and rearranging, gives the Donnan condition
(7)
An additional expression is needed to determine the micropore electrolyte potential ϕmi. That expression is the macropore-micropore current balance
(8)
In Equation 8, C is the volumetric capacitance (F/m3)
(9)
where Cref (F/m3) and α (the pore charge to surface capacitance ratio, F·m3/mol2) are experimentally determined parameters, and εmi is the volume fraction of micropores (dimensionless).
Finally, the volume-averaged source/sink term Ri (mol/(m3·s)), defining the flux of ions from macropore to micropore, is written as
(10)
Current Efficiency
The charge efficiency, H (dimensionless), of the process is determined, for the charging part of the process (between 1980 s and 5220 s), using the following expression
(11)
where Γ is the cumulative amount of salt stored in the micropores (mol) and Σ is the cumulative electrical charge passed through the external circuit (C). Specifically,
(12)
where Q is the inlet volumetric flow rate (m3/s), and
(13)
In Equation 13, the surface integral is evaluated only for the top boundary of the upper porous electrode, and is is the electrode current density (A/m2).
Results and Discussion
The flow through the cell is shown in Figure 2. The flow velocity in the porous electrodes is much lower than that in the main channel.
Figure 2: The flow through the cell.
The simulation starts with both terminals being connected without any external load (0 V applied potential). Next, at 1800 s, a positive voltage of 0.4 V is applied on the top (positive) porous electrode. The lower electrode is assumed to be connected to ground. Finally, the external applied voltage is again turned off. Figure 3 shows the flow averaged outlet concentration during the capacitive deionization process. The outlet concentration is highly time dependent.
Figure 3: The mixing-cup average outlet salt concentration during the simulation.
Figure 3 shows that during the charging process, as the micropore concentrations increase toward the value given by the Donnan condition, the outlet concentration gradually approaches the inlet value. The variation in the micropore concentration of Cl- during the simulation is shown in Figure 4. During the charging part of the simulation, the amount of Cl- increases in the upper (positive) electrode, while the lower (negative) electrode is emptied. The concentrations gradually approach a saturation value. When the voltage is again set to 0 V, the micropore concentrations approach a value similar to the starting concentrations.
Figure 4: The micropore amount of Cl- ions during the process.
Figure 5 shows the desalination ratio throughout the macropores and main channel after 3600 s. At the outlet, the salt concentration has been reduced to about half of the inlet concentration.
Figure 5: The desalination ratio in the electrolyte, in macropores and the main channel, after 3600 s. The desalination ratio varies significantly between different parts of the porous electrodes.
The maximum and minimum values of the local attractive potential, μatt, and capacitance are shown in Figure 6 and Figure 7, respectively, both for the upper electrode. The two quantities vary significantly during the simulation.
Figure 6: The maximum and minimum values for the local attractive potential in the upper electrode.
Figure 7: The maximum and minimum values for the local capacitance in the upper electrode.
Finally, the charge efficiency during the charging part of the process is shown in Figure 8. The efficiency increases toward a value of 70%.
Figure 8: Charge efficiency during desalination.
Notes About the COMSOL Implementation
As the deionization process is assumed to be completely non-Faradaic, the Porous Electrode Reaction subnode under the Tertiary Current Distribution, Nernst–Planck, Porous Electrode 1 node is disabled.
The coupling between micro- and macropores is carried out through a Reactions node, where the flux according to Equation 10 is entered.
Equation 5 is solved for in the Domain ODEs and DAEs: Attractive Potential interface, and Equation 8 is solved for in the Domain ODEs and DAEs: Micropore Electrolyte Potential interface. The charge efficiency (Equation 11) is computed as a postprocessing operation, in the Evaluation Group: Charge Efficiency evaluation group.
References
1. P.M. Biesheuvel, S. Porada, M. Levi, and M.Z. Bazant, “Attractive forces in microporous carbon electrodes for capacitive deionization,” J. Solid State Electrochem., vol. 18, pp. 1365–1376, 2014.
2. A. Hemmatifar, M. Stadermann, and J.G. Santiago, “Two-Dimensional Porous Electrode Model for Capacitive Deionization,” J. Phys. Chem. C, vol. 119, pp. 24681–24694, 2015.
Application Library path: Electrochemistry_Module/Electrochemical_Engineering/capacitive_deionization
Modeling Instructions
From the File menu, choose New.
New
In the New window, click  Model Wizard.
Model Wizard
1
In the Model Wizard window, click  2D.
2
In the Select Physics tree, select Electrochemistry > Tertiary Current Distribution, Nernst–Planck > Tertiary, Electroneutrality (tcd).
3
Click Add.
4
In the Concentrations (mol/m³) table, enter the following settings:
cNap
cClm
5
In the Select Physics tree, select Fluid Flow > Porous Media and Subsurface Flow > Brinkman Equations (br).
6
Click Add.
7
In the Select Physics tree, select Mathematics > ODE and DAE Interfaces > Domain ODEs and DAEs (dode).
8
Click Add.
9
In the Field name (1) text field, type phil_micro.
10
In the Dependent variables (1) table, enter the following settings:
phil_micro
11
Click  Define Dependent Variable Unit.
12
In the Dependent variable quantity table, enter the following settings:
Dependent variable quantity
Unit
Custom unit
V
13
In the Source term quantity table, enter the following settings:
Source term quantity
Unit
Custom unit
A/m^3
14
In the Select Physics tree, select Mathematics > ODE and DAE Interfaces > Domain ODEs and DAEs (dode).
15
Click Add.
16
In the Source term quantity table, enter the following settings:
Source term quantity
Unit
Custom unit
1
17
In the Field name (1) text field, type mu.
18
In the Dependent variables (1) table, enter the following settings:
mu
19
Click  Study.
20
In the Select Study tree, select General Studies > Stationary.
21
Click  Done.
Global Definitions
Geometry
1
In the Model Builder window, under Global Definitions click Parameters 1.
2
In the Settings window for Parameters, type Geometry in the Label text field.
3
Locate the Parameters section. Click  Load from File.
4
Browse to the model’s Application Libraries folder and double-click the file capacitive_deionization_geometry.txt.
Electrolyte
1
In the Home toolbar, click  Parameters and choose Add > Parameters.
2
In the Settings window for Parameters, type Electrolyte in the Label text field.
3
Locate the Parameters section. Click  Load from File.
4
Browse to the model’s Application Libraries folder and double-click the file capacitive_deionization_electrolyte.txt.
Porous Electrode
1
In the Home toolbar, click  Parameters and choose Add > Parameters.
2
In the Settings window for Parameters, type Porous Electrode in the Label text field.
3
Locate the Parameters section. Click  Load from File.
4
Browse to the model’s Application Libraries folder and double-click the file capacitive_deionization_porous_electrode.txt.
Porous Electrode, Attractive Potential
1
In the Home toolbar, click  Parameters and choose Add > Parameters.
2
In the Settings window for Parameters, type Porous Electrode, Attractive Potential in the Label text field.
3
Locate the Parameters section. Click  Load from File.
4
Browse to the model’s Application Libraries folder and double-click the file capacitive_deionization_porous_electrode_attractive_potential.txt.
Flow
1
In the Home toolbar, click  Parameters and choose Add > Parameters.
2
In the Settings window for Parameters, type Flow in the Label text field.
3
Locate the Parameters section. Click  Load from File.
4
Browse to the model’s Application Libraries folder and double-click the file capacitive_deionization_flow.txt.
Process
1
In the Home toolbar, click  Parameters and choose Add > Parameters.
2
In the Settings window for Parameters, type Process in the Label text field.
3
Locate the Parameters section. Click  Load from File.
4
Browse to the model’s Application Libraries folder and double-click the file capacitive_deionization_process.txt.
Mesh
1
In the Home toolbar, click  Parameters and choose Add > Parameters.
2
In the Settings window for Parameters, type Mesh in the Label text field.
3
Locate the Parameters section. Click  Load from File.
4
Browse to the model’s Application Libraries folder and double-click the file capacitive_deionization_mesh.txt.
Next, define the variables, a function, nonlocal couplings and view settings needed for the model. Some of the variables and nonlocal couplings will show warnings, since we have not yet defined the geometry.
Definitions
Donnan Condition and Current Balance
1
In the Home toolbar, click  Variables and choose Local Variables.
2
In the Settings window for Variables, type Donnan Condition and Current Balance in the Label text field.
3
Locate the Variables section. Click  Load from File.
4
Browse to the model’s Application Libraries folder and double-click the file capacitive_deionization_donnan_condition_current_balance.txt.
Outlet Concentration
1
In the Home toolbar, click  Variables and choose Local Variables.
2
In the Settings window for Variables, type Outlet Concentration in the Label text field.
3
Locate the Variables section. In the table, enter the following settings:
Name
Expression
Unit
Description
c_out_Clm
intop_outlet(tcd.ntflux_cClm*d_z)/intop_outlet(u*d_z)
 
Outlet concentration, Cl-
c_out_Nap
intop_outlet(tcd.ntflux_cNap*d_z)/intop_outlet(u*d_z)
 
Outlet concentration, Na+
Porous Electrode, Attractive Potential
1
In the Home toolbar, click  Variables and choose Local Variables.
2
In the Settings window for Variables, type Porous Electrode, Attractive Potential in the Label text field.
3
Locate the Variables section. In the table, enter the following settings:
Name
Expression
Unit
Description
csum
cNap_micro+cClm_micro
mol/m³
Sum of micropore concentrations
mu_att
E/csum/root.k_B_const/T/root.N_A_const
 
Nonelectrostatic adsorption potential
Porous Electrode, Stern Layer Properties
1
In the Home toolbar, click  Variables and choose Local Variables.
2
In the Settings window for Variables, type Porous Electrode, Stern Layer Properties in the Label text field.
3
Locate the Variables section. In the table, enter the following settings:
Name
Expression
Unit
Description
C
C_ref+alpha*charge_sum^2
s4·A²/(kg·m5)
Stern layer capacitance
Storage
1
In the Home toolbar, click  Variables and choose Local Variables.
2
In the Settings window for Variables, type Storage in the Label text field.
3
Locate the Variables section. In the table, enter the following settings:
Name
Expression
Unit
Description
n_Clm_upper
intop_upper_electrode(cClm_micro*d_z)
 
Amount of Cl-, upper electrode
n_Clm_lower
intop_lower_electrode(cClm_micro*d_z)
 
Amount of Cl-, lower electrode
Define two integration operators, one for the upper and one for the lower electrode. They will be used to compute the total amount of stored Cl- in each electrode.
Integration 1 (intop1)
1
In the Definitions toolbar, click  Nonlocal Couplings and choose Integration.
2
In the Settings window for Integration, type intop_upper_electrode in the Operator name text field.
Integration 2 (intop_upper_electrode2)
1
Right-click Integration 1 (intop_upper_electrode) and choose Duplicate.
2
In the Settings window for Integration, type intop_lower_electrode in the Operator name text field.
Define a maximum and a minimum operator. They will be used to plot the maximum and minimum value of the local attractive potential and local capacitance in the upper electrode.
Maximum 1 (maxop1)
1
In the Definitions toolbar, click  Nonlocal Couplings and choose Maximum.
2
In the Settings window for Maximum, type maxop_upper in the Operator name text field.
Minimum 1 (minop1)
1
In the Definitions toolbar, click  Nonlocal Couplings and choose Minimum.
2
In the Settings window for Minimum, type minop_upper in the Operator name text field.
Define an integration operator for the outlet. It will be used to compute the mixing-cup average outlet concentration of salt.
Integration, outlet boundary
1
In the Definitions toolbar, click  Nonlocal Couplings and choose Integration.
2
In the Settings window for Integration, type Integration, outlet boundary in the Label text field.
3
In the Operator name text field, type intop_outlet.
Axis
1
In the Model Builder window, expand the Component 1 (comp1) > Definitions > View 1 node, then click Axis.
2
In the Settings window for Axis, locate the Axis section.
3
From the View scale list, choose Manual.
4
In the y scale text field, type 30.
Geometry 1
Rectangle 1 (r1)
1
In the Geometry toolbar, click  Rectangle.
2
In the Settings window for Rectangle, locate the Size and Shape section.
3
In the Width text field, type L_electrode.
4
In the Height text field, type h_porous.
Rectangle 2 (r2)
1
In the Geometry toolbar, click  Rectangle.
2
In the Settings window for Rectangle, locate the Size and Shape section.
3
In the Width text field, type L_electrode.
4
In the Height text field, type h_gap.
5
Locate the Position section. In the y text field, type h_porous.
Rectangle 3 (r3)
1
In the Geometry toolbar, click  Rectangle.
2
In the Settings window for Rectangle, locate the Size and Shape section.
3
In the Width text field, type L_electrode.
4
In the Height text field, type h_porous.
5
Locate the Position section. In the y text field, type h_porous+h_gap.
Rectangle 4 (r4)
1
In the Geometry toolbar, click  Rectangle.
2
In the Settings window for Rectangle, locate the Size and Shape section.
3
In the Width text field, type L_pipe.
4
In the Height text field, type h_gap.
5
Locate the Position section. In the x text field, type -L_pipe.
6
In the y text field, type h_porous.
Rectangle 5 (r5)
1
In the Geometry toolbar, click  Rectangle.
2
In the Settings window for Rectangle, locate the Size and Shape section.
3
In the Width text field, type L_pipe.
4
In the Height text field, type h_gap.
5
Locate the Position section. In the x text field, type L_electrode.
6
In the y text field, type h_porous.
Ignore Edges 1 (ige1)
1
In the Geometry toolbar, click  Virtual Operations and choose Ignore Edges.
2
In the Settings window for Ignore Edges, locate the Input section.
3
Click the  Paste Selection button for Edges to ignore.
4
In the Paste Selection dialog, type 6, 12 in the Selection text field.
5
Click OK.
6
In the Geometry toolbar, click  Build All.
7
Click the  Zoom Extents button in the Graphics toolbar.
Now that the geometry has been defined, set the correct selections on the nonlocal coupling operators.
Definitions
Integration 1 (intop_upper_electrode)
1
In the Model Builder window, under Component 1 (comp1) > Definitions click Integration 1 (intop_upper_electrode).
2
Select Domain 3 only.
Integration 2 (intop_lower_electrode)
1
In the Model Builder window, click Integration 2 (intop_lower_electrode).
2
Select Domain 2 only.
Maximum 1 (maxop_upper)
1
In the Model Builder window, click Maximum 1 (maxop_upper).
2
Select Domain 3 only.
Minimum 1 (minop_upper)
1
In the Model Builder window, click Minimum 1 (minop_upper).
2
Select Domain 3 only.
Integration, outlet boundary (intop_outlet)
1
In the Model Builder window, click Integration, outlet boundary (intop_outlet).
2
In the Settings window for Integration, locate the Source Selection section.
3
From the Geometric entity level list, choose Boundary.
4
Select Boundary 14 only.
All the warnings should now have disappeared.
Next, add material properties for the saline solution.
Add Material
1
In the Home toolbar, click  Add Material to open the Add Material window.
2
Go to the Add Material window.
3
In the tree, select Built-in > Water, liquid.
4
Click the Add to Component button in the window toolbar.
5
In the Home toolbar, click  Add Material to close the Add Material window.
Materials
Water, liquid (mat1)
A few properties are not defined in the water material properties, but we will define those later in the corresponding physics interfaces instead.
Now, define the settings in the physics interfaces
Tertiary Current Distribution, Nernst–Planck (tcd)
1
In the Model Builder window, under Component 1 (comp1) click Tertiary Current Distribution, Nernst–Planck (tcd).
2
In the Settings window for Tertiary Current Distribution, Nernst–Planck, locate the Out-of-Plane Thickness section.
3
In the dz text field, type d_z.
Species Charges 1
1
In the Model Builder window, under Component 1 (comp1) > Tertiary Current Distribution, Nernst–Planck (tcd) click Species Charges 1.
2
In the Settings window for Species Charges, locate the Charge section.
3
In the zcNap text field, type z_Nap.
4
In the zcClm text field, type z_Clm.
Electrolyte 1
1
In the Model Builder window, click Electrolyte 1.
2
In the Settings window for Electrolyte, locate the Convection section.
3
From the u list, choose Velocity field (br).
4
Locate the Diffusion section. In the DcNap text field, type DNa.
5
In the DcClm text field, type DCl.
Initial Values 1
1
In the Model Builder window, click Initial Values 1.
2
In the Settings window for Initial Values, locate the Initial Values section.
3
In the cClm text field, type c_NaCl.
Porous Electrode 1
1
In the Physics toolbar, click  Domains and choose Porous Electrode.
2
Select Domains 2 and 3 only.
3
In the Settings window for Porous Electrode, locate the Convection section.
4
From the u list, choose Velocity field (br).
5
Locate the Diffusion section. In the DcNap text field, type DNa.
6
In the DcClm text field, type DCl.
7
Locate the Electrode Current Conduction section. From the σs list, choose User defined.
Define an arbitrary, but high, electrical conductivity for the porous electrode. In the associated text field, type 100.
8
Locate the Porous Matrix Properties section. In the εs text field, type eps_s.
9
In the εl text field, type 1-eps_s.
Porous Electrode Reaction 1
Disable the Porous Electrode Reaction subnode. Faradaic reactions will not be taken into account in the present model.
1
In the Model Builder window, right-click Porous Electrode Reaction 1 and choose Disable.
Inflow 1
1
In the Physics toolbar, click  Boundaries and choose Inflow.
2
Select Boundary 1 only.
3
In the Settings window for Inflow, locate the Concentration section.
4
In the c0,cClm text field, type root.c_NaCl.
Outflow 1
1
In the Physics toolbar, click  Boundaries and choose Outflow.
2
Select Boundary 14 only.
Electric Ground 1
1
In the Physics toolbar, click  Boundaries and choose Electric Ground.
2
Select Boundary 5 only.
Load Cycle 1
1
In the Physics toolbar, click  Boundaries and choose Load Cycle.
2
Select Boundary 9 only.
3
In the Settings window for Load Cycle, locate the Load Type section.
4
From the list, choose Potentiostatic.
5
Locate the Continuation Conditions section. Select the Use elapsed time only checkbox.
Voltage 1
1
In the Physics toolbar, click  Attributes and choose Voltage.
2
In the Settings window for Voltage, locate the Continuation Conditions section.
3
Select the Elapsed time checkbox.
4
In the tmax text field, type 30[min].
Load Cycle 1
In the Model Builder window, click Load Cycle 1.
Voltage 2
1
In the Physics toolbar, click  Attributes and choose Voltage.
2
In the Settings window for Voltage, locate the Voltage section.
3
In the Eset text field, type V_max.
4
Locate the Continuation Conditions section. Select the Elapsed time checkbox.
Load Cycle 1
In the Model Builder window, click Load Cycle 1.
Voltage 3
In the Physics toolbar, click  Attributes and choose Voltage.
Reactions 1
1
In the Physics toolbar, click  Domains and choose Reactions.
To couple the micropore and macropore concentrations in the porous electrodes, add a Reactions node as follows:
2
Select Domains 2 and 3 only.
3
In the Settings window for Reactions, locate the Reaction Rates section.
4
In the RcNap text field, type R_Nap_macro.
5
In the RcClm text field, type R_Clm_macro.
6
Click to expand the Reacting Volume section. From the list, choose Total volume.
Brinkman Equations (br)
Porous Matrix 1
1
In the Model Builder window, under Component 1 (comp1) > Brinkman Equations (br) > Porous Medium 1 click Porous Matrix 1.
2
In the Settings window for Porous Matrix, locate the Matrix Properties section.
3
From the εp list, choose User defined. In the associated text field, type eps_s.
4
From the κ list, choose User defined. In the associated text field, type kappa_s.
Initial Values 1
1
In the Model Builder window, under Component 1 (comp1) > Brinkman Equations (br) click Initial Values 1.
2
In the Settings window for Initial Values, locate the Initial Values section.
3
Specify the u vector as
u_in
x
Fluid Properties 1
1
In the Physics toolbar, click  Domains and choose Fluid Properties.
2
Select Domain 1 only.
Inlet 1
1
In the Physics toolbar, click  Boundaries and choose Inlet.
2
Select Boundary 1 only.
3
In the Settings window for Inlet, locate the Boundary Condition section.
4
From the list, choose Fully developed flow.
5
Locate the Fully Developed Flow section. In the Uav text field, type u_in.
Outlet 1
1
In the Physics toolbar, click  Boundaries and choose Outlet.
2
Select Boundary 14 only.
Enable Normal flow to simulate that the pipe continues on the other side of the system boundary.
3
In the Settings window for Outlet, locate the Pressure Conditions section.
4
Select the Normal flow checkbox.
Define the settings of the Global ODEs and DAEs interface used to compute the micropore electrolyte potential as follows:
Domain ODEs and DAEs: Micropore Electrolyte Potential
1
In the Model Builder window, under Component 1 (comp1) click Domain ODEs and DAEs (dode).
2
In the Settings window for Domain ODEs and DAEs, type Domain ODEs and DAEs: Micropore Electrolyte Potential in the Label text field.
3
Select Domains 2 and 3 only.
4
Click to expand the Dependent Variables section.
Distributed ODE 1
1
In the Model Builder window, under Component 1 (comp1) > Domain ODEs and DAEs: Micropore Electrolyte Potential (dode) click Distributed ODE 1.
2
In the Settings window for Distributed ODE, locate the Source Term section.
3
In the f text field, type iv_dl_helmholtz-iv_dl_diffuse.
4
Locate the Damping or Mass Coefficient section. In the da text field, type 0.
Define the settings of the Domain ODEs and DAEs interface used to compute the attractive potential.
Domain ODEs and DAEs: Attractive Potential
1
In the Model Builder window, under Component 1 (comp1) click Domain ODEs and DAEs 2 (dode2).
2
In the Settings window for Domain ODEs and DAEs, type Domain ODEs and DAEs: Attractive Potential in the Label text field.
3
Select Domains 2 and 3 only.
1
In the Model Builder window, under Component 1 (comp1) > Domain ODEs and DAEs: Attractive Potential (dode2) click Distributed ODE 1.
2
In the Settings window for Distributed ODE, locate the Source Term section.
3
In the f text field, type mu_att-mu.
4
Locate the Damping or Mass Coefficient section. In the da text field, type 0.
Initial Values 1
1
In the Model Builder window, click Initial Values 1.
2
In the Settings window for Initial Values, locate the Initial Values section.
3
In the mu text field, type mu_0.
As the geometry has only 90 degree angles, set up a mapped mesh.
Mesh 1
1
In the Model Builder window, under Component 1 (comp1) click Mesh 1.
2
In the Settings window for Mesh, locate the Sequence Type section.
3
From the list, choose User-controlled mesh.
Size
1
In the Model Builder window, under Component 1 (comp1) > Mesh 1 click Size.
2
In the Settings window for Size, locate the Element Size section.
3
From the Calibrate for list, choose General physics.
4
From the Predefined list, choose Normal.
Boundary Layers 1, Corner Refinement 1, Free Triangular 1, Size 1
1
In the Model Builder window, under Component 1 (comp1) > Mesh 1, Ctrl-click to select Size 1, Corner Refinement 1, Free Triangular 1, and Boundary Layers 1.
2
Right-click and choose Delete.
Mapped 1
In the Mesh toolbar, click  Mapped.
Distribution 1
1
In the Mesh toolbar, click  Distribution.
2
Select Boundaries 7 and 12 only.
3
In the Settings window for Distribution, locate the Distribution section.
4
From the Distribution type list, choose Predefined.
5
In the Number of elements text field, type elem_porous_count.
6
In the Element ratio text field, type 2.
Distribution 2
1
In the Mesh toolbar, click  Distribution.
2
Select Boundaries 4 and 10 only.
3
In the Settings window for Distribution, locate the Distribution section.
4
From the Distribution type list, choose Predefined.
5
In the Number of elements text field, type elem_porous_count.
6
In the Element ratio text field, type 2.
7
Select the Reverse direction checkbox.
Distribution 3
1
In the Mesh toolbar, click  Distribution.
2
Select Boundaries 5 and 9 only.
3
In the Settings window for Distribution, locate the Distribution section.
4
From the Distribution type list, choose Predefined.
5
In the Number of elements text field, type elem_length_count.
Distribution 4
1
In the Mesh toolbar, click  Distribution.
2
Select Boundaries 1 and 14 only.
3
In the Settings window for Distribution, locate the Distribution section.
4
From the Distribution type list, choose Predefined.
5
In the Number of elements text field, type elem_pipe_count.
6
In the Element ratio text field, type 2.
7
Select the Symmetric distribution checkbox.
Study 1
Step 1: Stationary
1
In the Model Builder window, under Study 1 click Step 1: Stationary.
2
In the Settings window for Stationary, locate the Physics and Variables Selection section.
3
In the Solve for column of the table, under Component 1 (comp1), clear the checkboxes for Tertiary Current Distribution, Nernst–Planck (tcd), Domain ODEs and DAEs: Micropore Electrolyte Potential (dode), and Domain ODEs and DAEs: Attractive Potential (dode2).
Step 2: Time Dependent
1
In the Study toolbar, click  Time Dependent.
2
In the Settings window for Time Dependent, locate the Study Settings section.
3
In the Output times text field, type 0 range(30[min],1[min], 2[h]).
4
Locate the Physics and Variables Selection section. In the Solve for column of the table, under Component 1 (comp1), clear the checkbox for Brinkman Equations (br).
To reduce the computation time, change to the Constant (Newton) nonlinear method, on the Fully Coupled node.
Solution 1 (sol1)
1
In the Study toolbar, click  Show Default Solver.
2
In the Model Builder window, expand the Solution 1 (sol1) node.
3
In the Model Builder window, expand the Study 1 > Solver Configurations > Solution 1 (sol1) > Time-Dependent Solver 1 node, then click Fully Coupled 1.
4
In the Settings window for Fully Coupled, click to expand the Method and Termination section.
5
From the Nonlinear method list, choose Constant (Newton).
The model setup is now complete. Compute the study.
6
In the Study toolbar, click  Compute.
Next, add a single-select configuration, for 3600 s, and improve the default plots.
Results
Single-Select Solution 1
1
In the Results toolbar, click  Configurations and choose Single-Select Solution.
2
In the Settings window for Single-Select Solution, locate the Solution section.
3
From the Time (s) list, choose Interpolation.
4
In the Time text field, type 3600.
Electrolyte Potential (tcd)
1
In the Model Builder window, under Results click Electrolyte Potential (tcd).
2
In the Settings window for 2D Plot Group, locate the Data section.
3
From the Solution parameters list, choose From configuration.
Surface 1
1
In the Model Builder window, expand the Electrolyte Potential (tcd) node, then click Surface 1.
2
In the Settings window for Surface, locate the Coloring and Style section.
3
From the Color table list, choose Cynanthus.
Arrow Surface 1
1
In the Model Builder window, click Arrow Surface 1.
2
In the Settings window for Arrow Surface, locate the Coloring and Style section.
3
Select the Scale factor checkbox. In the associated text field, type 0.005.
4
In the Electrolyte Potential (tcd) toolbar, click  Plot.
Electrolyte Current Density (tcd)
In the Model Builder window, under Results right-click Electrolyte Current Density (tcd) and choose Delete.
Electrode Potential with Respect to Ground (tcd)
1
In the Model Builder window, expand the Results > Electrode Potential with Respect to Ground (tcd) node, then click Electrode Potential with Respect to Ground (tcd).
2
In the Settings window for 2D Plot Group, locate the Data section.
3
From the Solution parameters list, choose From configuration.
Surface 1
1
In the Model Builder window, click Surface 1.
2
In the Settings window for Surface, locate the Coloring and Style section.
3
From the Color table list, choose Cynanthus.
Arrow Surface 1
1
In the Model Builder window, click Arrow Surface 1.
2
In the Settings window for Arrow Surface, locate the Coloring and Style section.
3
Select the Scale factor checkbox. In the associated text field, type 0.005.
4
In the Electrode Potential with Respect to Ground (tcd) toolbar, click  Plot.
Electrode Current Density (tcd)
In the Model Builder window, under Results right-click Electrode Current Density (tcd) and choose Delete.
Concentration, Nap (tcd)
1
In the Settings window for 2D Plot Group, locate the Data section.
2
From the Solution parameters list, choose From configuration.
Surface 1
1
In the Model Builder window, expand the Concentration, Nap (tcd) node, then click Surface 1.
2
In the Settings window for Surface, locate the Coloring and Style section.
3
From the Color table list, choose Baptisia.
Streamline 1
1
In the Model Builder window, click Streamline 1.
2
In the Settings window for Streamline, locate the Streamline Positioning section.
3
From the Positioning list, choose On selected boundaries.
4
In the Number text field, type 16.
5
Locate the Selection section. Click to select the  Activate Selection toggle button.
6
Select Boundary 1 only.
7
Locate the Coloring and Style section. Find the Point style subsection.
8
Select the Scale factor checkbox. In the associated text field, type 0.5.
9
In the Concentration, Nap (tcd) toolbar, click  Plot.
Concentration, Clm (tcd)
1
In the Model Builder window, under Results click Concentration, Clm (tcd).
2
In the Settings window for 2D Plot Group, locate the Data section.
3
From the Solution parameters list, choose From configuration.
Surface 1
1
In the Model Builder window, expand the Concentration, Clm (tcd) node, then click Surface 1.
2
In the Settings window for Surface, locate the Coloring and Style section.
3
From the Color table list, choose Baptisia.
Streamline 1
1
In the Model Builder window, click Streamline 1.
2
In the Settings window for Streamline, locate the Streamline Positioning section.
3
From the Positioning list, choose On selected boundaries.
4
In the Number text field, type 16.
5
Locate the Selection section. Click to select the  Activate Selection toggle button.
6
Select Boundary 1 only.
7
Locate the Coloring and Style section. Find the Point style subsection.
8
Select the Scale factor checkbox. In the associated text field, type 0.5.
9
In the Concentration, Clm (tcd) toolbar, click  Plot.
Velocity (br)
1
In the Model Builder window, under Results click Velocity (br).
2
In the Settings window for 2D Plot Group, locate the Data section.
3
From the Solution parameters list, choose From configuration.
Surface
1
In the Model Builder window, expand the Velocity (br) node, then click Surface.
2
In the Settings window for Surface, locate the Coloring and Style section.
3
From the Color table list, choose Kyanite.
Velocity (br)
In the Velocity (br) toolbar, click  Arrow Line.
Arrow Line 1
1
In the Settings window for Arrow Line, click Replace Expression in the upper-right corner of the Expression section. From the menu, choose Component 1 (comp1) > Brinkman Equations > Velocity and pressure > u,v - Velocity field.
2
Locate the Arrow Positioning section. From the Placement list, choose Uniform anisotropic.
3
In the Number of arrows text field, type 20.
4
Locate the Coloring and Style section.
5
Select the Scale factor checkbox. In the associated text field, type 20.
6
From the Color list, choose Gray.
Selection 1
1
In the Velocity (br) toolbar, click  Selection.
2
In the Settings window for Selection, locate the Selection section.
3
Click  Paste Selection.
4
In the Paste Selection dialog, type 1 14 in the Selection text field.
5
Click OK.
6
In the Velocity (br) toolbar, click  Plot.
Pressure (br)
1
In the Model Builder window, under Results click Pressure (br).
2
In the Settings window for 2D Plot Group, locate the Data section.
3
From the Solution parameters list, choose From configuration.
Domain ODEs and DAEs: Micropore Electrolyte Potential
1
In the Model Builder window, click Domain ODEs and DAEs: Micropore Electrolyte Potential.
2
In the Settings window for 2D Plot Group, locate the Data section.
3
From the Solution parameters list, choose From configuration.
Surface 1
1
In the Model Builder window, expand the Domain ODEs and DAEs: Micropore Electrolyte Potential node, then click Surface 1.
2
In the Settings window for Surface, locate the Coloring and Style section.
3
From the Color table list, choose Cynanthus.
Domain ODEs and DAEs: Attractive Potential
1
In the Model Builder window, under Results click Domain ODEs and DAEs: Attractive Potential.
2
In the Settings window for 2D Plot Group, locate the Data section.
3
From the Solution parameters list, choose From configuration.
Surface 1
1
In the Model Builder window, expand the Domain ODEs and DAEs: Attractive Potential node, then click Surface 1.
2
In the Settings window for Surface, locate the Coloring and Style section.
3
From the Color table list, choose Arctium.
Salinity Ratio
1
In the Results toolbar, click  2D Plot Group.
2
In the Settings window for 2D Plot Group, type Salinity Ratio in the Label text field.
3
Locate the Data section. From the Solution parameters list, choose From configuration.
Surface 1
1
In the Salinity Ratio toolbar, click  Surface.
2
In the Settings window for Surface, locate the Expression section.
3
In the Expression text field, type cClm/c_NaCl.
4
Locate the Coloring and Style section. From the Color table list, choose MetasepiaBlue.
5
In the Salinity Ratio toolbar, click  Plot.
Mixing-Cup Average Outlet Salt Concentration
1
In the Home toolbar, click  Add Plot Group and choose 1D Plot Group.
2
In the Settings window for 1D Plot Group, type Mixing-Cup Average Outlet Salt Concentration in the Label text field.
3
Click to expand the Title section. From the Title type list, choose None.
4
Locate the Plot Settings section.
5
Select the y-axis label checkbox. In the associated text field, type Average outlet concentration (mol/m<sup>3</sup>), green.
6
Select the Two y-axes checkbox.
7
Select the Secondary y-axis label checkbox. In the associated text field, type Applied potential (V), blue.
8
Locate the Legend section. Clear the Show legends checkbox.
Applied Potential
1
In the Mixing-Cup Average Outlet Salt Concentration toolbar, click  Global.
2
In the Settings window for Global, type Applied Potential in the Label text field.
3
Locate the y-Axis section. Select the Plot on secondary y-axis checkbox.
4
Click Replace Expression in the upper-right corner of the y-Axis Data section. From the menu, choose Component 1 (comp1) > Tertiary Current Distribution, Nernst–Planck > Load Cycle 1 > tcd.lc1.E_app - Applied voltage - V.
Mixing-Cup Average Outlet Salt Concentration
In the Mixing-Cup Average Outlet Salt Concentration toolbar, click  Global.
Outlet Concentration
1
In the Settings window for Global, type Outlet Concentration in the Label text field.
2
Locate the y-Axis Data section. Click  Clear Table.
3
In the table, enter the following settings:
Expression
Unit
Description
c_out_Clm
mol/m^3
Outlet concentration, Cl-
4
In the Mixing-Cup Average Outlet Salt Concentration toolbar, click  Plot.
Mixing-Cup Average Outlet Salt Concentration and mu
1
In the Model Builder window, right-click Mixing-Cup Average Outlet Salt Concentration and choose Duplicate.
2
In the Model Builder window, click Mixing-Cup Average Outlet Salt Concentration 1.
3
In the Settings window for 1D Plot Group, type Mixing-Cup Average Outlet Salt Concentration and mu in the Label text field.
4
Locate the Plot Settings section. In the y-axis label text field, type Average outlet concentration (mol/m<sup>3</sup>).
5
In the Secondary y-axis label text field, type Attractive potential (1).
6
Locate the Legend section. Select the Show legends checkbox.
7
From the Position list, choose Middle left.
Applied Potential
In the Model Builder window, right-click Applied Potential and choose Delete.
Outlet Concentration
1
In the Settings window for Global, click to expand the Legends section.
2
From the Legends list, choose Manual.
3
In the table, enter the following settings:
Legends
c<sub>out, Cl-</sub>
Mixing-Cup Average Outlet Salt Concentration and mu
In the Mixing-Cup Average Outlet Salt Concentration and mu toolbar, click  Global.
Max and min \mu_att, upper electrode
1
In the Settings window for Global, type Max and min \mu_att, upper electrode in the Label text field.
2
Locate the y-Axis section. Select the Plot on secondary y-axis checkbox.
3
Locate the y-Axis Data section. In the table, enter the following settings:
Expression
Unit
Description
maxop_upper(mu_att)
1
max(\mu_att)
minop_upper(mu_att)
1
min(\mu_att)
4
Click to expand the Coloring and Style section. Find the Line style subsection. From the Line list, choose Dashed.
5
In the Mixing-Cup Average Outlet Salt Concentration and mu toolbar, click  Plot.
Mixing-Cup Average Outlet Salt Concentration and C
1
In the Model Builder window, right-click Mixing-Cup Average Outlet Salt Concentration and mu and choose Duplicate.
2
In the Settings window for 1D Plot Group, type Mixing-Cup Average Outlet Salt Concentration and C in the Label text field.
3
Locate the Plot Settings section. In the Secondary y-axis label text field, type Capacitance (MF/m<sup>3</sup>).
Max and min C, upper electrode
1
In the Model Builder window, expand the Mixing-Cup Average Outlet Salt Concentration and C node, then click Max and min \mu_att, upper electrode.
2
In the Settings window for Global, type Max and min C, upper electrode in the Label text field.
3
Locate the y-Axis Data section. In the table, enter the following settings:
Expression
Unit
Description
maxop_upper(C)
MF/m^3
max(C)
minop_upper(C)
MF/m^3
min(C)
4
In the Mixing-Cup Average Outlet Salt Concentration and C toolbar, click  Plot.
Define an evaluation group, used to compute the amount of salt stored and the charge passed through the external circuit. Together, these two quantities define the charge efficiency of the cell.
Evaluation Group: Charge Efficiency
1
In the Results toolbar, click  Evaluation Group.
2
In the Settings window for Evaluation Group, type Evaluation Group: Charge Efficiency in the Label text field.
Cumulative Amount of Salt Entering the Cell
1
In the Evaluation Group: Charge Efficiency toolbar, click  Global Evaluation.
2
In the Settings window for Global Evaluation, type Cumulative Amount of Salt Entering the Cell in the Label text field.
3
Locate the Expressions section. In the table, enter the following settings:
Expression
Unit
Description
Q*c_NaCl
mol/s
Inlet molar flow
4
Locate the Data Series Operation section. From the Transformation list, choose Integral.
5
Select the Cumulative checkbox.
Evaluation Group: Charge Efficiency
In the Evaluation Group: Charge Efficiency toolbar, click  Integration and choose Line Integration.
Cumulative Amount of Salt Exiting the Cell
1
In the Settings window for Line Integration, type Cumulative Amount of Salt Exiting the Cell in the Label text field.
2
Select Boundary 14 only.
3
Locate the Expressions section. Click  Clear Table.
4
In the table, enter the following settings:
Expression
Unit
Description
tcd.ntflux_cClm*d_z
mol/s
Outlet molar flow
5
Locate the Data Series Operation section. From the Transformation list, choose Integral.
6
Select the Cumulative checkbox.
Evaluation Group: Charge Efficiency
In the Evaluation Group: Charge Efficiency toolbar, click  Integration and choose Line Integration.
Cumulative Charge
1
In the Settings window for Line Integration, type Cumulative Charge in the Label text field.
2
Select Boundary 9 only.
3
Locate the Expressions section. Click  Clear Table.
4
In the table, enter the following settings:
Expression
Unit
Description
tcd.nIs*d_z
A
Boundary current
5
Locate the Data Series Operation section. From the Transformation list, choose Integral.
6
Select the Cumulative checkbox.
Evaluation Group: Charge Efficiency
1
In the Model Builder window, click Evaluation Group: Charge Efficiency.
2
In the Settings window for Evaluation Group, locate the Transformation section.
3
From the Transformation type list, choose General.
4
In the Expression text field, type -(gev1-int1)*F_const/int2*100.
5
In the Column header text field, type Charge Efficiency.
6
In the Evaluation Group: Charge Efficiency toolbar, click  Evaluate.
Evaluation Group: Charge Efficiency
1
Go to the Evaluation Group: Charge Efficiency window.
2
Click the Table Graph button in the window toolbar.
Results
Table Graph 1
1
In the Settings window for Table Graph, locate the Data section.
2
From the x-axis data list, choose Time (s).
3
From the Plot columns list, choose Manual.
4
In the Columns list box, select Charge Efficiency (1).
Charge Efficiency
1
In the Model Builder window, under Results click 1D Plot Group 17.
2
In the Settings window for 1D Plot Group, type Charge Efficiency in the Label text field.
3
Locate the Plot Settings section.
4
Select the y-axis label checkbox. In the associated text field, type Charge Efficiency (%).
5
Locate the Axis section. Select the Manual axis limits checkbox.
6
In the x minimum text field, type 1980.
7
In the x maximum text field, type 5220.
8
In the y minimum text field, type 0.
9
In the y maximum text field, type 100.
10
In the Charge Efficiency toolbar, click  Plot.
Amount of Cl- Stored in Micropores
1
In the Results toolbar, click  1D Plot Group.
2
In the Settings window for 1D Plot Group, type Amount of Cl- Stored in Micropores in the Label text field.
3
Locate the Title section. From the Title type list, choose None.
4
Locate the Plot Settings section.
5
Select the x-axis label checkbox. In the associated text field, type Time (s).
6
Select the y-axis label checkbox. In the associated text field, type Amount of Cl- stored in micropores (\mu mol).
7
Select the Two y-axes checkbox.
8
Select the Secondary y-axis label checkbox. In the associated text field, type Applied potential (V).
9
Locate the Legend section. From the Layout list, choose Outside graph axis area.
10
From the Position list, choose Bottom.
Applied Potential
1
In the Amount of Cl- Stored in Micropores toolbar, click  Global.
2
In the Settings window for Global, type Applied Potential in the Label text field.
3
Locate the y-Axis section. Select the Plot on secondary y-axis checkbox.
4
Click Replace Expression in the upper-right corner of the y-Axis Data section. From the menu, choose Component 1 (comp1) > Tertiary Current Distribution, Nernst–Planck > Load Cycle 1 > tcd.lc1.E_app - Applied voltage - V.
Amount of Cl- Stored in Micropores
In the Amount of Cl- Stored in Micropores toolbar, click  Global.
Amount Stored
1
In the Settings window for Global, type Amount Stored in the Label text field.
2
Locate the y-Axis Data section. Click  Clear Table.
3
In the table, enter the following settings:
Expression
Unit
Description
n_Clm_upper
umol
Upper electrode
n_Clm_lower
umol
Lower electrode
4
In the Amount of Cl- Stored in Micropores toolbar, click  Plot.