Case Study — Electrochemical Library

Getting Started

This case study explores more advanced topics of building custom Simscape™ libraries. It uses an example library for modeling electrochemical systems. The library introduces a new electrochemical domain and defines all of the fundamental components required to build electrochemical models, including an electrochemical reference, through and across sensors, sources, and a cross-domain component. The example illustrates some of the salient features of Physical Networks modeling, such as selection of Through and Across variables and how power is converted between domains. We suggest that you work through the previous section, Case Study — Basic Custom Block Library, before looking at this more advanced example.

The example library comes built and on your path so that it is readily executable. However, it is recommended that you copy the source files to a new folder, for which you have write permission, and add that folder to your MATLAB^® path. This will allow you to make changes and rebuild the library for yourself. The source files for the example library are in the following namespace folder:

matlabroot/toolbox/physmod/simscape/supporting_files/example_libraries/+ElectroChem

where matlabroot is the MATLAB root folder on your machine, as returned by entering

matlabroot

in the MATLAB Command Window.

After copying the files, change the folder name +ElectroChem to another name, for example +MyElectroChem, so that your copy of the library builds with a unique name.

Building the Custom Library

To build the library, type

sscbuild MyElectroChem

in the MATLAB Command Window. If building from within the +MyElectroChem namespace folder, you can omit the argument and type just

sscbuild

When the build completes, open the generated library by typing

MyElectroChem_lib

For more information on the library build process, see Building Custom Block Libraries.

Defining a New Domain

Simscape software comes with several Foundation domains, such as mechanical translational, mechanical rotational, electrical, and so on. Where possible, use these predefined domains. For example, when creating new electrical components, use the Foundation electrical domain foundation.electrical.electrical. This ensures that your components can be connected to the standard Simscape blocks.

As an example of an application requiring the addition of a new domain, consider a battery where the underlying equations involve both electrical and chemical processes [1].

Electrochemical Battery Driving a Resistive Load R

Two half-cells are separated by a membrane that prevents the ions flowing between cells, and hence electrons flow from the solid lead anode to the platinum cathode. The two half-cell reactions are:

$P b \leftrightarrow P b^{2 +} + 2 e^{-}$

$F e^{2 +} \leftrightarrow F e^{3 +} + e^{-}$

The current results in the lead being oxidized and the iron being reduced, with the overall reaction given by:

$P b + 2 F e^{3 +} \leftrightarrow P b^{2 +} + 2 F e^{2 +}$

The chemical reaction can be modeled using the network concepts of Through and Across variables (for details, see Basic Principles of Modeling Physical Networks). The Through variable represents flow, and the Across variable represents effort. When selecting the Through and Across variables, you should use SI units and the product of the two variables is usually chosen to have units of power.

In the electrochemical reactions, an obvious choice for the Through variable is the molar flow rate $\dot{n}$ of ions, measured in SI units of mol/s. The corresponding Across variable is called chemical potential, and must have units of J/mol to ensure that the product of Through and Across variables has units of power, J/s. The chemical potential or Gibb’s free energy per mol is given by:

$μ = μ_{0} + R T \ln a$

where μ₀ is the standard state chemical potential, R is the perfect gas constant, T is the temperature, and a is the activity. In general, the activity can be a function of a number of different parameters, including concentration, temperature, and pressure. Here it is assumed that the activity is proportional to the molar concentration defined as number of moles of solute divided by the mass of solvent.

To see the electrochemical domain definition, open the Simscape file +MyElectroChem/ElectroChem.ssc.

domain ElectroChem
% Electrochemical Domain
% Define through and across variables for the electrochemical domain

% Copyright 2008-2014 The MathWorks, Inc.

    variables
        % Chemical potential
        mu = { 1.0 'J/mol' };
    end

    variables(Balancing = true)
        % Molar flow
        ndot = { 1.0 'mol/s' };
    end

end

The molar fundamental dimension and unit is predefined in the Simscape unit registry. If it had not been, then you could have added it with:

pm_adddimension('mole','mol')

Note

The case study does not involve modifying the domain definition. However, if at a later point you decide to modify this electrochemical domain definition and see the effect of these modifications on the custom components, make sure to change all the node declarations in the component files from ElectroChem.ElectroChem to MyElectroChem.ElectroChem. Otherwise, your component files still point to the domain definition in the original namespace.

Structuring the Library

It is good practice to structure a library by adding hierarchy. To do this, you can subdivide the namespace folder into subfolders, each subfolder name starting with the + character. If you look at the +MyElectroChem folder, you will see that it has subfolders +Elements, +Sensors, and +Sources. Open the library by typing MyElectroChem_lib, and you will see the three corresponding sublibraries.

Defining a Reference Component

A physical network must have a reference block, against which Across variables are measured. So, for example, the Foundation library contains the Electrical Reference block for the electrical domain, Mechanical Rotational Reference block for the rotational mechanical domain, and so on. The electrochemical zero chemical potential is defined by the component file +MyElectroChem/+Elements/Reference.ssc.

component Reference
% Chemical Reference
% Port A is a zero chemical potential reference port.

% Copyright 2008-2016 The MathWorks, Inc.

    nodes
        A = ElectroChem.ElectroChem; % A:top
    end

    connections
        connect(A, *);
    end

end

The component has one electrochemical port, named A, located at the top of the block icon.

The component uses a connection to an implicit reference node:

connect(A, *);

For more information on component connections and the implicit reference node syntax, see Connections to Implicit Reference Node.

Defining an Ideal Source Component

An ideal Across source provides a constant value for the Across variable regardless of the value of the Through variable. In the electrical domain, this corresponds to the DC Voltage Source block in the Foundation library. In the example library, the component file +MyElectroChem/+Sources/ChemPotentialSource.ssc implements the equivalent source for the chemical domain.

component ChemPotentialSource
% Constant Potential Source
% Provides a constant chemical potential between ports A and B.

% Copyright 2008-2013 The MathWorks, Inc.

    nodes
        A = ElectroChem.ElectroChem; % A:top
        B = ElectroChem.ElectroChem; % B:bottom
    end

    parameters
        mu0 = {0, 'J/mol'}; % Chemical potential
    end

    variables(Access=private)
        ndot = { 0, 'mol/s' }; % Molar flow rate
    end

    branches
        ndot: A.ndot -> B.ndot; % Through variable ndot from node A to node B
    end

    equations
        let
            mu = A.mu - B.mu; % Across variable from A to B
        in
            mu == mu0;
        end
    end

end

The dual of an ideal Across source is an ideal Through source, which maintains the Through variable to some set value regardless of the value of the Across variable. In the electrical domain, this corresponds to the DC Current Source block in the Foundation library. In the example library, this source is not implemented.

Defining Measurement Components

Every domain requires both a Through and an Across measurement block. In the example library, the component file +MyElectroChem/+Sensors/SensorThrough.ssc implements a molar flow rate sensor.

component SensorThrough
% Molar Flow Sensor
% Returns the value of the molar flow between the A and the B port
% to the physical signal port PS.

% Copyright 2008-2013 The MathWorks, Inc.

    nodes
        A = ElectroChem.ElectroChem; % A:top
        B = ElectroChem.ElectroChem; % B:bottom
    end

    outputs
        out  = { 0, 'mol/s' }; % PS:top
    end

    variables(Access=private)
        ndot = { 0, 'mol/s' }; % Molar flow rate
    end

    branches
        ndot: A.ndot -> B.ndot; % Through variable ndot from node A to node B
    end

    equations
        let
            mu = A.mu - B.mu; % Across variable from A to B
        in
            mu == 0;     % No potential drop
            out == ndot; % Equate value of molar flow to PS output
        end
    end

end

The flow rate is presented as a Physical Signal, which can then in turn be passed to Simulink via a PS-Simulink Converter block. The branches section and the let statement in the equation section define the relationship between Through and Across variables for the sensor. In this case, an ideal flow sensor has zero potential drop, that is mu == 0, where mu is the chemical potential. The second equation assigns the value of the Through variable to the Physical Signal output.

The component file +MyElectroChem/+Sensors/SensorAcross.ssc implements a chemical potential sensor.

component SensorAcross
% Chemical Potential Sensor
% Returns the value of the chemical potential across the A and B ports
% to the physical signal port PS.

% Copyright 2008-2013 The MathWorks, Inc.

    nodes
        A = ElectroChem.ElectroChem; % A:top
        B = ElectroChem.ElectroChem; % B:bottom
    end

    outputs
        out  = { 0, 'J/mol' }; % PS:top
    end

    variables(Access=private)
        ndot = { 0, 'mol/s' }; % Molar flow rate
    end

    branches
        ndot: A.ndot -> B.ndot; % Through variable ndot from node A to node B
    end

    equations
        let
            mu = A.mu - B.mu; % Across variable from A to B
        in
            ndot == 0; % Draws no molar flow
            out == mu; % Equate value of chemical potential difference to PS output
        end
    end

end

The chemical potential is presented as a Physical Signal, which can then in turn be passed to Simulink via a PS-Simulink Converter block. The branches section and the let statement in the equation section define the relationship between Through and Across variables for the sensor. In this case, an ideal chemical potential sensor draws no flow, that is ndot == 0, where ndot is the flow rate. The second equation assigns the value of the Across variable to the Physical Signal output.

Defining Basic Components

Having created the measurement and reference blocks, the next step is to create blocks that define behavioral relationships between the Through and Across variables. In the electrical domain, for example, such components are resistor, capacitor, and inductor.

As an example of a basic electrochemical component, consider the chemical reduction or oxidation of an ion, which can be thought of as the electrochemical equivalent of a nonlinear capacitor. The defining equations in terms of Through and Across variables ν and μ are:

$\dot{n} = ν$

$a = \frac{n}{C_{0} M}$

$μ = μ_{0} + R T \ln a$

where n is the number of moles of the ion, C₀ is the standard concentration of 1 mol/kg, and M is the mass of the solute.

To see the implementation of these equations, open the file +MyElectroChem/+Elements/ChemEnergyStore.ssc.

component ChemEnergyStore
% Chemical Energy Store :1 :fixed
% Represents a solution of dissolved ions. The port A presents the
% chemical potential defined by mu0 + log(n/(C0*M))*R*T where mu0 is the
% standard state oxidizing potential, n is the number of moles of the ion,
% C0 is the standard concentration of 1 mol/kg, M is the mass of solvent,
% R is the universal gas constant, and T is the temperature.

% Copyright 2008-2015 The MathWorks, Inc.

    nodes
        A = ElectroChem.ElectroChem; % A:top
    end

    parameters
        mu0 = {-7.42e+04, 'J/mol'}; % Standard state oxidizing potential
        m_solvent = {1, 'kg'};      % Mass of solvent
        T = {300, 'K'};             % Temperature
    end

    parameters (Access=private)
        R = {8.314472, '(J/K)/mol'}; % Universal gas constant
        C0 = {1, 'mol/kg'};          % Standard concentration
        n1 = {1e-10, 'mol'};         % Minimum number of moles
    end

    variables
        ndot = { 0, 'mol/s' }; % Molar flow rate
        n  = {value = { 0.01, 'mol' }, priority = priority.high}; % Quantity of ions
    end

    branches
        ndot : A.ndot -> *; % Through variable ndot
    end

    equations
        n.der == ndot;
        if n > n1
            A.mu == mu0 + log(n/(C0*m_solvent))*R*T;
        else
            A.mu == mu0 + (log(n1/(C0*m_solvent)) + n/n1 - 1)*R*T;
        end
    end

end

This component introduces two Simscape language features not yet used in the blocks looked at so far. These are:

Use of a conditional statement in the equation section. This is required to prevent taking the logarithm of zero. Hence if the molar concentration is less than the specified level n1, then the operand of the logarithm function is limited. Without this protection, the solver could perturb the value of n to zero or less.
Definition of private parameters that can be used in the equation section. Here the Universal Gas constant (R) and the Standard Concentration (C0) are defined as private parameters. Their values could equally well be used directly in the equations, but this would reduce readability of the definition. Similarly, the lower limit on the molar concentration n1 is also defined as a private parameter, but could equally well have been exposed to the user.

Defining a Cross-Domain Interfacing Component

Cross-domain blocks allow the interchange of energy between domains. For example, the Rotational Electromechanical Converter block in the Foundation library converts between electrical and rotational mechanical energy. To relate the two sets of Through and Across variables, two equations are required. The first comes from an underlying physical law, and the second from summing the powers from the two domains into the converter, which must total zero.

As an example of an interfacing component, consider the electrochemical half-cell. The chemical molar flow rate and the electrical current are related by Faraday’s law, which requires that:

$ν = \frac{i}{z F}$

where ν is the molar flow rate, i is the current, z is the number of electrons per ion, and F is the Faraday constant. The second equation comes from equating the electrical and chemical powers:

$(V_{2} - V_{1}) i = (μ_{2} - μ_{1}) ν$

which can be rewritten as:

$(V_{2} - V_{1}) = (μ_{2} - μ_{1}) \frac{ν}{i} = \frac{μ_{2} - μ_{1}}{z F}$

This is the Nernst equation written in terms of chemical potential difference, (μ₂ – μ₁). These chemical-electrical converter equations are implemented by the component file +MyElectroChem/+Elements/Chem2Elec.ssc.

component Chem2Elec
% Chemical to Electrical Converter
% Converts chemical energy into electrical energy (and vice-versa)
% assuming no losses. The electrical current flow i is related to the
% molar flow of electrons ndot by i = -ndot*z*F where F is the Faraday
% constant and z is the number of exchanged electrons.

% Copyright 2008-2017 The MathWorks, Inc.

    nodes
        p = foundation.electrical.electrical; % +:top
        n = foundation.electrical.electrical; % -:top
        A = ElectroChem.ElectroChem;  % A:bottom
        B = ElectroChem.ElectroChem;  % B:bottom
    end

    parameters
        z = {1, '1'};                % Number of exchanged electrons
    end

    parameters(Access=private)
        F = {9.6485309e4, 'C/mol'};  % Faraday constant
    end

    variables
        i = { 0, 'A'   };      % Current
        ndot = { 0, 'mol/s' }; % Molar flow rate
    end

    branches
        i   : p.i    -> n.i;    % Through variable i from node p to node n
        ndot: A.ndot -> B.ndot; % Through variable ndot from node A to node B
    end

    equations
        let
            k = 1/(z*F);
            v = p.v - n.v;    % Across variable v from p to n
            mu = A.mu - B.mu; % Across variable mu from A to B
        in
            v == k*mu;    % From equating power
            ndot == -k*i; % Balance electrons (Faraday's Law)
        end
    end

end

Note the use of the let-in-end construction in the component equations. An intermediate term k is declared as

$k = \frac{1}{z F}$

It is then used in both equations in the expression clause that follows.

This component has four ports but only two equations. This is because the component interfaces two different physical networks. Each of the networks has two ports and one equation, thus satisfying the requirement for n–1 equations, where n is the number of ports. In the case of a cross-domain component, the two equations are coupled, thereby defining the interaction between the two physical domains.

The Faraday constant is a hidden parameter, because it is a physical constant that block users would not need to change. Therefore, it will not appear in the block dialog box generated from the component file.

Customizing the Appearance of the Library

The library can be customized using lib.m files. A lib.m file located in the top-level namespace folder can be used to add annotations. The name of the top-level library model is constructed automatically during the build process based on the top-level namespace name, as namespace_lib, but you can add a more descriptive name to the top-level library as an annotation. For example, open +MyElectroChem/lib.m in the MATLAB Editor. The following line annotates the top-level library with its name:

libInfo.Annotation = sprintf('Example Electrochemical Library')

In the electrochemical library example, lib.m files are also placed in each subnamespace folder to customize the name and appearance of respective sublibraries. For example, open +MyElectroChem/+Sensors/lib.m in the MATLAB Editor. The following line causes the sublibrary to be named Electrochemical Sensors:

libInfo.Name = 'Electrochemical Sensors';

In the absence of the lib.m file, the library would be named after the subnamespace name, that is, Sensors. For more information, see Library Configuration Files.

Using the Custom Components to Build a Model

The Battery Cell with Custom Electrochemical Domain example uses the electrochemical library to model a lead-iron battery. See the example help for further information.

References

[1] Pêcheux, F., B. Allard, C. Lallement, A. Vachoux, and H. Morel. “Modeling and Simulation of Multi-Discipline Systems using Bond Graphs and VHDL-AMS.” International Conference on Bond Graph Modeling and Simulation (ICBGM). New Orleans, USA, 23–27 Jan. 2005.