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 directory, for
which you have write permission, and add that directory 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 package
directory:

matlabroot/toolbox/physmod/simscape/simscapedemos/+ElectroChem

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

matlabroot

in the MATLAB Command Window.

After copying the files, change the directory name `+ElectroChem`

to another name, for example `+MyElectroChem`

, so that your copy of
the library builds with a unique name.

To build the library, type

ssc_build MyElectroChem

in the MATLAB Command Window. If building from within the
`+MyElectroChem`

package directory, you can omit the argument
and type just

ssc_build

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.

Simscape software comes with several Foundation domains, such as mechanical
translational, mechanical rotational, electrical, hydraulic, 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:

$$Pb\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:

$$Pb+2F{e}^{3+}\leftrightarrow P{b}^{2+}+2F{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:

$$\mu ={\mu}_{0}+RT\mathrm{ln}a$$

where μ_{0} 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')

It is good practice to structure a library by adding hierarchy. To do this, you
can subdivide the package directory into subdirectories, each subdirectory name
starting with the `+`

character. If you look at the
`+MyElectroChem`

directory, you will see that it has
subdirectories `+Elements`

, `+Sensors`

, and
`+Sources`

. Open the library by typing
`MyElectroChem_lib`

, and you will see the three corresponding
sublibraries.

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.

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.

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.

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}=\nu $$

$$a=\frac{n}{{C}_{0}M}$$

$$\mu ={\mu}_{0}+RT\mathrm{ln}a$$

where *n* is the number of moles of the ion,
*C _{0}* is the standard concentration
of 1 mol/kg, and

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.

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:

$$\nu =\frac{i}{zF}$$

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:

$$\left({V}_{2}-{V}_{1}\right)i=\left({\mu}_{2}-{\mu}_{1}\right)\nu $$

which can be rewritten as:

$$\left({V}_{2}-{V}_{1}\right)=\left({\mu}_{2}-{\mu}_{1}\right)\frac{\nu}{i}=\frac{{\mu}_{2}-{\mu}_{1}}{zF}$$

This is the Nernst equation written in terms of chemical potential difference,
(μ_{2} – μ_{1}). 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}{zF}$$

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.

The library can be customized using `lib.m`

files. A
`lib.m`

file located in the top-level package directory 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 package name, as

, but you can add a
more descriptive name to the top-level library as an annotation. For example, open
* package*_lib

`+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 subpackage directory 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 subpackage name, that is, `Sensors`

. For more
information, see Library Configuration Files.

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.

[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.