Continuous Electrode Models and Application of Exact Schemes in Modeling of Electrical Impedance Measurements

Abstract

The crucial issue in electrical impedance (EI) measurements lies in the galvanic interaction between the electrodes and the investigated material. This paper brings together the basic and applied research experience and combines their results with excellent properties. Consequently, innovative precise methodologies have emerged, enabling the direct modeling of EI measurements, free from the inaccuracies often associated with numerical approaches. As an outcome of the efficiency and robustness of the applied method, the conductivity of the material and the electrodes are represented by a common piecewise function, which is used to solve the differential equation modeling of the EI measurement. Moreover, this allows the possibility for modeling the conductivity of electrodes with continuous functions, providing an important generalization of the Complete Electrode Model (CEM), which has been widely used so far. The effectiveness of the novel approach was showcased through two distinct case studies. In the first case study, potential functions within both the material and the electrodes were computed using the CEM. In the second case study, calculations were performed utilizing the newly introduced continuous electrode model. The simulation results suggest that the new method is a powerful tool for biological research, from in vitro experiments to animal studies and human applications.

1. Introduction

The main goal of this paper is to combine the measurement and theoretical research results in a joint study in order to introduce a completely new electrical impedance (EI) modeling approach. This new modeling method is now able to entirely represent the measurement circuit, including electrodes and electrode–material interactions, constructed during the EI measurements. Moreover, due to the use of a unique mathematical method, the resulting procedure is able to calculate the potentials in the electrodes and in the material without the errors that are common in numerical computation methods. This provides a completely unique and new basis for EI measurement modeling, both for understanding the behavior of electrodes and for reconstructing the impedance of the material.

The relevance of EI measurements is increasing [1,2]. This non-destructive technique detects the physicochemical properties of the measured material and its variations [3,4]. Several approaches for the implementation of EI measurements are currently being developed: discrete frequency EI, EI spectrum (EIS) measurement, and EI tomography (EIT) [3]. The combination of all these methods offers interesting and useful implementation possibilities. For example, combining EIS and EIT can contribute to overcoming difficulties in the image reconstruction with the EIT method [4,5]. An example of such a technology is the multi-frequency EIT (mfEIT), where the application-specific prototypes have been developed beyond the basic research [5,6]. The motivation of this research is to develop a completely new mathematical method in order to implement the mfEIT. Since this is a strongly ill-conditioned, non-linear, and unstable inverse problem, even very small perturbations of the voltage and/or potential values used in the model can cause significant variations in the results [3,4]. Therefore, it is important to apply models that accurately describe the physics of the measurement method. Thus, errors caused by model inaccuracies in the reconstructed impedance images can be minimized.

Regardless of which type of EI measurement is used, the implementation is always based on similar principles [3,4]. Electrodes are placed on the surface of the investigated material and an excitation signal is applied to the selected electrodes. The parameters of the generated electric field are measured and, based on the obtained values of the measured parameters, conclusions are drawn about the complex electrical impedance of the material. Further, based on the results, the physicochemical properties of the material can be analyzed [2,3]. The generator used in the measurement may be either a current generator or a voltage generator [3]. The most commonly used signal type for multi-frequency approaches is the monochromatic sine signal [7,8].

Based on these considerations, it is evident that the interaction between the electrode and the investigated material is crucial for the reliability of the analysis, since the galvanic contact between the electrode and the surface of the material is required to close the circuit utilized for the measurement [3]. The role of the electrode is also essential from another point of view, since electrodes are made of a conductive material and, due to the presence of galvanic connections, they significantly influence the electric field in the material [9]. Incorrect recognition and modeling of the electrodes (and the interaction between electrode and material) significantly reduces the reliability of the EI measurement since the resulting distortions and artifacts are propagated into the result generation, causing errors [10,11]. The state-of-the-art electrode modeling technique is the Complete Electrode Model (CEM), which is the basis for various EIT imaging algorithms [12,13,14].

The basic idea behind the CEM is that the electrical properties of the electrode, as a function of its material, are considered as an impedance in a series connection with the material under measurement [11]. Consequently, the voltage measured on an electrode is modeled as the sum of two components—the voltages across the electrode and the measured material. Naturally, this effect is significant if the current passing through the electrode is not zero [11]. The CEM is a convenient model for the implementation of EIT numerical methods; however, the CEM approach can be easily applied to EIS and more general EI methods [4]. The implementation of the CEM can be described as follows. First, the elliptic differential equation modeling of the electric field is discretized with the finite element method over the given domain, and then the potential values representing the boundary conditions are substituted with the values expressed from the voltages across the electrodes [11]. In this way, the system of linear equations obtained with the domain discretization under the consideration directly represents the electrode impedance (in the form of a concentrated parameter) and the measured potential value. The potential value arises from the electrode–material interaction and thus it is eliminated from the model [3,4,11].

In contrast, the approach presented in this paper, beyond the possibilities exploited in the case of the CEM, is able to model the electrodes and the complete measurement circuit used in EI measurements without any errors caused by the numerical modeling. This is explained by the possibility that the method presented in this paper applies the same differential equation not only to the material, but also to the electrodes and even to the material–electrode interaction. This raises the possibility that not only the material, but also the electrodes are represented by continuous functions in the model. As a consequence, the values of the analytical solution of the differential equation, which is the physical representation of the measurement, are calculated at any point in the measurement assembly, regardless of whether the electrode, material, or interaction between these two is being investigated. Therefore, the main contribution of this research is the improvement of the efficacy and robustness of EI measurements—from spectroscopy to tomography—or even in hybrid technologies.

This publication is structured as follows. The Section 1 contains the introduction, the Section 2 describes the related works, the Section 3 introduces a completely new modeling approach, the continuous electrode model. The Section 4 describes the case studies comparing the results obtained using the CEM and the continuous model, while conclusions are drawn in the Section 5.

2. Related Works

The behavior of a wide variety of electrodes used to perform EI measurements often cannot be represented by simple RC elements. Electrodes may also have properties that change over time; therefore, special techniques are needed to compensate for this behavior [15]. Considering the relevant literature, it can be concluded that, in general, the multi-frequency approach is the one that represents the most electrode artefacts, and therefore most studies are related to this issue [16]. The most common application of multifrequency measurements is in the field of human body composition [4,17]. In this context, a large number of studies have been published on electrode–skin interaction properties and their potential errors [18,19]. Since in EI measurements, electrodes with large surfaces are usually used, electrode contact effects were investigated in a separate study. Significant work on understanding the failures has been undertaken by Hwang et al. [20], who investigated the effects of temperature variation, changes in oxygen partial pressure, and other mechanical effects on the contact impedance of the electrode. The failure effects of imperfect contact on solid electrolytes have also been investigated by Fleig et al. [21]. In addition, studies have been carried out to understand and eliminate the error phenomena caused by electrode mismatches in EI measurements [22,23]. The detection and elimination of these errors is crucial for the success of EIS measurements [24].

In the previous study, Vizvari et al. [8] presented a new EI prototype and a specific data collection procedure on simple RC elements, which provided impressive results even at an early stage. In order to implement the technology in praxis, further research is needed on the degree of uncertainty, and additional steps may be required to create a technology suited for practical utilization. In the next step of the development, the use of electrodes complicates the implementation of the method. The electrodes naturally involve the previously described risk of errors; therefore, it is essential to describe and model the electrode–matter interaction in detail.

Moreover, great attention is paid to the application of the most efficient computational procedures for modeling EI measurements. Vizvari et al. [5,25] have established a completely new basis for the mathematical modeling procedures required for their interpretation. One of the achievements is the absolutely new mathematical approach, where Vizvari et al. [25] introduced the exact scheme for second-order ordinary differential equations (ODEs) using arbitrary spatial discretization. Exact schemes were also first introduced by Vizvari et al. [25] for second-order ODEs with a self-adjoint differential operator. The exact schemes are characterized by the property that they always provide the values of the analytic solution of the ODEs in the grid points, independently of the spatial discretization. The efficiency and robustness of this outstanding ability and the mathematical structures have been demonstrated in detail by Vizvari et al. [25].

Moreover, the previous studies have also introduced a specific property of exact schemes. Since they are based on Local Green’s functions defined by solutions to the homogeneous ODE of the original problem, the exact schemes are perfectly suitable for difficult mathematical problems where the functions in the ODE have discontinuities [25]. All these outstanding properties are exploited in this study, where a completely new modeling aspect of EIS measurements is presented, based on the improvement of the CEM and on the use of exact schemes. In this completely new approach, the electrodes and the material under investigation are modeled using a piecewise conductivity function, which has the advantage that the electrode–material interaction can be modeled more directly and accurately. Moreover, the electrical behavior of the electrodes may be modeled using continuous functions in addition to RC elements. Utilizing all this, it is possible to model and visualize the frequency-dependent potential profiles in the material under investigation and in the electrodes with the least possible error.

3. The Continuous Electrode Model Using the Exact Schemes

The basic equation for the physical modeling of EI measurements is the ODE with the following self-adjoint operator [3]:

$$Du\left(x\right)=\frac{d}{dx}\left(\kappa \left(x\right)\frac{d}{dx}u\left(x\right)\right)=f\left(x\right)x\in [0,L],$$

(1)

where

–

D denotes the second order, self-adjoint differential operator;

–

$\kappa \left(x\right)\ge {\kappa }_0>0$ denotes a positive, isotropic conductivity function ($\mathrm{S}/\mathrm{m}$);

–

$u\left(x\right)$ denotes the electric potential function (V);

–

$f\left(x\right)$ denotes the source function ($\mathrm{A}/{\mathrm{m}}^2$).

It is assumed that the classical solution $u\left(x\right)$ exists on $[0,L]$ with the appropriate boundary conditions.

Vizvari et al. have implemented EI measurements by specifying the following mixed boundary conditions [8]:

$${\mathfrak{B}}_{DN}=\left\{\kappa \left(0\right)u\left(0\right)=i,u\left(L\right)=0\right\}.$$

(2)

The exact scheme for Equation (1) is defined by the following result, specifying the boundary conditions in Equation (2).

Theorem 1 (Exact scheme for Dirichlet and Neumann boundaries).

Let

$$x_0=0<x_1<\cdots <x_{i-1}<x_i<x_{i+1}<\cdots <x_{n-1}<x_n<x_{n+1}=L$$

(3)

be an arbitrary discretization of the interval $[0,L]$ into $(n+1)$ subintervals. Let ${\psi }_{i-1}\left(x\right)={\int }_{x_{i-1}}^x\frac{1}{\kappa \left(s\right)}ds$ and ${\phi }_i\left(x\right)={\int }_x^{x_{i+1}}\frac{1}{\kappa \left(s\right)}ds$ be the test functions obtained from the following initial value problems:

$$\{\begin{array}{ccc}D{\psi }_{i-1}\left(x\right)=0,\hfill & & \hfill \mathrm{(4a)}\\ {\psi }_{i-1}\left(x_{i-1}\right)=0,\hfill & & \hfill \mathrm{(4b)}\\ \kappa \left(x_{i-1}\right)\frac{d}{dx}{\psi }_{i-1}\left(x_{i-1}\right)=1,\hfill & & \hfill \mathrm{(4c)}\end{array}$$

$$\{\begin{array}{ccc}D{\psi }_i\left(x\right)=0,\hfill & & \hfill \mathrm{(5a)}\\ {\psi }_i\left(x_{i+1}\right)=0,\hfill & & \hfill \mathrm{(5b)}\\ \kappa \left(x_{i+1}\right)\frac{d}{dx}{\psi }_i\left(x_{i+1}\right)=-1,\hfill & & \hfill \mathrm{(5c)}\end{array}$$

where $i=1,2,\cdots ,n$. By using these structures and considerations, the following system of $n\ge 3$ linear equations is constructed:

$$\left\{\begin{array}{cc}\left(a_0+a_1\right)u_1-a_1{u}_2& =a_{0}u\left(0\right)+a_0{G}_0+a_1{H}_1,\hfill \\ -a_{i-1}{u}_{i-1}+\left(a_{i-1}+a_i\right)u_i-a_i{u}_{i+1}& =a_{i-1}{G}_{i-1}+a_i{H}_i,\hfill \\ -a_{n-1}{u}_{n-1}+\left(a_{n-1}+a_n\right)u_n& =a_{n}u\left(L\right)+a_{n-1}{G}_{n-1}+a_n{H}_n,\hfill \end{array}\right.$$

(6)

with indexes $i=2,3,\cdots ,(n-1)$. The coefficients in Equation (6) are defined as

$$a_{i-1}=\frac{1}{{\psi }_{i-1}\left(x_i\right)}=\frac{1}{{\phi }_i\left(x_{i-1}\right)},$$

(7)

$$G_{i-1}={\int }_{x_{i-1}}^{x_i}f\left(t\right){\psi }_{i-1}\left(t\right)dt,$$

(8)

$$H_i={\int }_{x_i}^{x_{i+1}}f\left(t\right){\phi }_i\left(t\right)dt,$$

(9)

where $i=1,2,3,\cdots ,n$. The $u\left(L\right)$ value can be substituted directly from Equation (2), while the $u\left(0\right)$ value is calculated using the Neumann-to-Dirichlet transformation:

$${\displaystyle u\left(0\right)=u\left(L\right)+{\phi }_n\left(0\right)\kappa \left(0\right)u\left(0\right)+{\int }_0^{L}f\left(t\right){\phi }_n\left(t\right)dt.}$$

(10)

The solution vector of Equation (6) $U={\left[u_1,u_2,\cdots ,u_n\right]}^T$ leads to the same values as the solution $u\left(x\right)$ of the second-order ODE (1) with boundary conditions (Equation (2)) at the interior grid points (Equation (3)) without any error; that is, $u\left(x_i\right)=u_i$.

A detailed proof of Theorem 1 can be found in Vizvari et al. [25] for further reading.

Theorem 1 is applied in cases where the material sample and the measurement assemblies are represented by a one-dimensional model. This is correct in all cases where the electric current density propagates in such geometry, where the perpendicular cross-section is constant along the entire length of the sample. Then, the electrodes are fixed to the two sides of the sample, where the surface area is equal to the whole perpendicular cross-section. In metrology practice, because of easier adaptation to sampling methods, the cylindrical geometry is preferred, for example in geophysical [26,27] or medical applications [28] of EIS. Based on this, a schematic illustration of the cylindrical geometry material sample and the associated measurement setup is shown in Figure 1.

The one-dimensional nature of the model problem depicted in Figure 1 arises from the specific geometry of the electrodes and the material, which possess a constant cross-sectional area denoted with A. As can be seen in Figure 1, the current generator (i) is connected to Electrode 1 with width a, the ground point ($u_{n+1}=0$) of the generator is connected to Electrode 2 with width b. As a consequence, it is easy to see from Theorem 1 that the integral functions, defined by Equations (4) and (5), represent the concentrated parameter derived from the impedance of the material. The potential values ($[u_0,u_1,u_2,u_3,\cdots ,u_i,\cdots ,u_n]$), which can be calculated using Theorem 1, are consequently the values of $u_i$ pertaining exclusively to the cross-sectional position $x_i$:

$$u\left(x_i\right)=u_i,i=1,2,\cdots ,n.$$

(11)

Based on Figure 1, the measurement method is modeled using the following statements defined in Equation (2). The application of these constraints, and the substitution of Equation (10), simplifies the corresponding system of linear equations in Equation (6) in the following, easily applicable, form:

$$\left\{\begin{array}{cc}\left(a_0+a_1\right)u_1-a_1{u}_2& =a_0{\phi }_n\left(0\right)i,\hfill \\ -a_{i-1}{u}_{i-1}+\left(a_{i-1}+a_i\right)u_i-a_i{u}_{i+1}& =0,\hfill \\ -a_{n-1}{u}_{n-1}+\left(a_{n-1}+a_n\right)u_n& =0,\hfill \end{array}\right.$$

(12)

with indexes $i=2,3,\cdots ,(n-1)$. The matrix in Equation (12) is symmetric and tridiagonal, and the sum of row elements is equal to zero, except the first and last rows. As a consequence of these advantageous properties, this reduced Laplacian matrix [5] is always invertible and Equation (12) has a unique $U={\left[u_1,u_2,u_3,\cdots ,u_n\right]}^T$ solution vector.

Corresponding to the model concept, since the domain is extended with the electrode lengths a and b, the $\kappa \left(x\right)$ function now describes the conductivity function on the whole $[0,L]$ interval.

$$\kappa (j\omega ,x)=\left\{\begin{array}{cc}{\kappa }_{el,1}(j\omega ,x),\hfill & ifx<a,\hfill \\ {\kappa }_m(j\omega ,x),\hfill & ifa\le x\le (L-b),\hfill \\ {\kappa }_{el,2}(j\omega ,x),\hfill & if(L-b)<x,\hfill \end{array}\right.$$

(13)

where

• $\kappa (j\omega ,x)$ is the complex admittance of the complete measured setup;
• ${\kappa }_{el,1}(j\omega ,x)$ is the complex admittance of the electrode placed on the left side of the investigated material sample;
• ${\kappa }_m(j\omega ,x)$ is the complex admittance of the investigated material sample;
• ${\kappa }_{el,2}(j\omega ,x)$ is the complex admittance of the electrode placed on the right side of the investigated material sample;
• $\omega$ is the angular frequency, $\omega =2\pi ·f$;
• f is the frequency; and
• $j=\sqrt{-1}$.

The impedance of the investigated material and the electrodes is defined in the model as the function $\kappa (j\omega ,x)$ in Equation (13). The functions ${\kappa }_{el,1}(j\omega ,x)$ and ${\kappa }_{el,2}(j\omega ,x)$, interpreted on $[0,a)$ and $\left(\right(L-b),L]$, respectively, represent the impedance of the electrodes applied to the measurement and they are placed on the surface of the material.

The handling of piecewise functions defined in Equation (13) is highly efficient when using the exact scheme defined in Theorem 1. In a related result by Vizvari et al. [25], it is demonstrated that, even for the piecewise function $\kappa \left(x\right)$, the values of the analytic solution in the mesh points can be calculated (if the function $\kappa \left(x\right)$ is integrable on the interval $[0,L]$). Moreover, the application of the exact scheme allows for the calculation of potential values at the discontinuities of the $\kappa \left(x\right)$ function. In a case study by Vizvari et al. [25], it has already been shown that the exact scheme can provide the analytical solution (in addition to the values recorded in the mesh points) in the simple case. The resulting analytic solution ($u\left(x\right)$) is continuous, but is not derivable at the discontinuity points of $\kappa \left(x\right)$ [25]. Naturally, if the differential operator D is applied to $u\left(x\right)$, the resulting discontinuities are eliminated, since $u\left(x\right)$ satisfies Equation (1). Therefore, for the piecewise conductivity function $\kappa \left(x\right)$, defined in Equation (13), these properties are expected using the exact scheme; however, in the case studies presented in Section 4, the symbolic construction of the analytical solution is omitted.

4. Brief Case Studies

The aim of the case studies is to investigate how the composition of the test material and the electrodes applied for the measurement is applicable to the modeling of the EIS measurement using the exact scheme described in Theorem 1. In this brief case study, two electrode modeling approaches are described:

1.

The electrodes are modeled with concentrated parameters, as is usual for the CEM;

2.

Modeling the impedance of electrodes with functions using continuous electrode models to include anomalies in the model more accurately.

For each of these case studies (based on Figure 1), the spatial discretization parameters are

$$x_{el,1}=[0,0.5,0.625,0.75,0.875,1],$$

(14)

$$x_m=[1,2,3,4,5,6,7,8,9,10],$$

(15)

$$x_{el,2}=[10,10.125,10.25,10.375,10.5,11],$$

(16)

which implies

$$a=1,b=1,L=11.$$

(17)

The parameters of the EIS measurement model are as follows:

$$f=[0.1,10,25.119,63.096,158.49,398.107,{10}^3,{10}^5],$$

(18)

and

$$\kappa \left(0\right)u\left(0\right)=i={10}^{-3} \mathrm{A},\mathrm{and}u\left(L\right)=0 \mathrm{V}.$$

(19)

Now consider the mathematical method used to model the electrical behavior of the material under investigation.

4.1. The Model of Measured Sample

The basic and well-known method for electrical modeling of the measured sample, the Cole–Cole model, can be used to describe the frequency-dependent behavior of the material [29]. The flexibility of the Cole–Cole model lies in describing the relaxation of the materials under investigation, whether it is described in terms of concentrated parameters or material properties (conductivity, resistivity, or dielectric constant) [28,29]. In the case studies, however, in addition to the frequency domain behavior of the sample under study, inhomogeneities along the x coordinate were considered in the model, since the resistivity of the material sample is used to build the model. Thus, the Cole–Cole parameters represent the physical properties of the material with continuous functions, which leads to the longitudinal behavior of the model. These functions are substituted into the Cole–Cole equation, whose variable is, by definition, the angular frequency. Finally, all these together yield a common model to represent the spatial and frequency domain properties of the material’s impedance. The extended mathematical model can be represented as follows:

$${\rho }_m(j\omega ,x)=\frac{1}{{\kappa }_m(j\omega ,x)}={\rho }_{\infty }\left(x\right)+\frac{{\rho }_0\left(x\right)-{\rho }_{\infty }\left(x\right)}{1+{\left(j\omega \tau \left(x\right)\right)}^{\alpha \left(x\right)}},$$

(20)

where

• ${\rho }_m(j\omega ,x)=\frac{1}{{\kappa }_m(j\omega ,x)}$ is the complex impedance of the material sample $\left(\mathsf{\Omega }\mathrm{m}\right)$;
• ${\rho }_{\infty }\left(x\right)=R_{\infty }\left(x\right)·\frac{A}{L}$ is the resistivity corresponding to the ∞ frequency $\left(\mathsf{\Omega }\mathrm{m}\right)$;
• ${\rho }_0\left(x\right)=R_0\left(x\right)·\frac{A}{L}$ is the resistivity corresponding to 0 Hz frequency $\left(\mathsf{\Omega }\mathrm{m}\right)$;
• $\tau \left(x\right)$ is the time constant $\left(\mathrm{s}\right)$;
• $\alpha \left(x\right)$ is the exponent parameter ($0<\alpha \left(x\right)\le 1$).

In the case of the location-dependent model parameters, inhomogeneities are modeled with the following functions:

$${\rho }_{\infty }\left(x\right)=10+10·e^{-{(x-2.5)}^2},$$

(21)

$${\rho }_0\left(x\right)=110+10·e^{-\frac{{(x-7.5)}^2}{4}},$$

(22)

$$\tau \left(x\right)={10}^{-3}+{10}^{-3}·e^{-\frac{{(x-4.5)}^2}{8}},$$

(23)

$$\alpha \left(x\right)=0.75+5·{10}^{-2}·e^{-\frac{{(x-7.5)}^2}{8}}.$$

(24)

Graphs of the Cole–Cole parameters defined in this way are illustrated in Figure 2.

In EIS measurements, the material sample cannot be accessed directly, hence the electrical behavior of the electrodes used to perform the measurements must be added to the mathematical model. Using the modeling approach provided by Equation (13), the reciprocal of the function in Equation (20) appears in the middle term of the piecewise function, i.e.,

$${\kappa }_m\left(x\right)=\frac{1}{{\rho }_m(j\omega ,x)}x\in [1,10].$$

(25)

The functions ${\kappa }_{el,1}\left(x\right)$ and ${\kappa }_{el,2}\left(x\right)$ are defined in the case of the complete and continuous electrode models in the next two subsections separately. The use of two different electrode modeling techniques leads to two different mathematical problems, which can be solved uniformly by applying Theorem 1 and Equation (12) respectively.

4.2. Complete Electrode Model Approach

The CEM considers the electrodes used for EIS measurement concentrated parameters. These concentrated parameters are included in the mathematical model in the form of resistance and capacitance, from which RC components can be obtained. Therefore, in case of the CEM approach, the electrodes in Figure 1 can be modeled by RC circuits. Figure 3 shows the changed electrodes in parallel RC circuits.

Figure 3 shows that the electrodes from Figure 1 are replaced by RC elements where the admittance of the new circuits can be obtained as follows:

$${\kappa }_{el,1}(j\omega ,x)=\frac{1}{R_{in}}+j\omega ·C_{in}={10}^{-2}+j\omega ·{10}^{-4}x\in [0,1],$$

(26)

and

$${\kappa }_{el,2}(j\omega ,x)=\frac{1}{R_{out}}+j\omega ·C_{out}=1.25\times {10}^{-2}+j\omega ·1.2\times {10}^{-4}x\in [10,11],$$

(27)

where

• $R_{in}$ is the resistance of the electrode on the input in Figure 3;
• $C_{in}$ is the capacitance of the electrode on the input in Figure 3;
• $R_{out}$ is the resistance of the electrode on the output in Figure 3;
• $C_{out}$ is the capacitance of the electrode on the output in Figure 3.

It is easy to observe from Equations (26) and (27) that, in the case of the CEM, the admittance of the electrodes depends only on the frequency, while the concentrated parameter values are constant along the whole length of the electrode.

In this case, the function $\kappa (j\omega ,x)$ is constructed using Equations (26) and (27). The reciprocal of the constructed piecewise function (impedance function) is defined in Equation (28) and is shown in Figure 4.

$${\rho }_{CEM}(j\omega ,x)=\frac{1}{{\kappa }_{CEM}(j\omega ,x)}=\left\{\begin{array}{cc}\frac{1}{{\kappa }_{el,1}(j\omega ,x)},\hfill & \mathrm{if}0<x<1,\hfill \\ \frac{1}{{\kappa }_m(j\omega ,x)},\hfill & \mathrm{if}1\le x\le 10,\hfill \\ \frac{1}{{\kappa }_{el,2}(j\omega ,x)},\hfill & \mathrm{if}10<x<11.\hfill \end{array}\right.$$

(28)

Naturally, the main properties of the impedance function ${\rho }_{CEM}(j\omega ,x)$ (Equation (28)) can be deduced from Equations (26) and (27); however, the graphical representation of the most important aspects that influence the measurement is even more illustrative in Figure 4. Figure 4 clearly shows that, based on the impedance values of the same color and considered at the same frequency, the electrodes always produce different values to the material under investigation. Further, it can be observed that, while the impedance of the material is strongly spatially dependent, the electrodes are represented as a constant function. Focusing on the frequency dependence, it is noticed that the impedance of the electrodes decreases significantly as the frequency increases, hence the influence of the electrodes is relevant only at low frequencies. It is also remarkable that the difference between the electrodes placed on the input and the output is also highlighted.

The next step in solving the case study using the CEM is to calculate the potentials based on Theorem 1 and on the derived system of the linear equation in Equation (12). The first step is to calculate the values of $a_{i-1}$ (based on Equation (7)), which is easily carried out using the discretization (Equations (14)–(16)) and the function ${\rho }_{CEM}$ (Equation (28)). Therefore, the matrix in Equation (12) can be constructed from the values $a_{i-1}$, while the right-hand side can be constructed using the measurement properties defined in Equation (19). The solution of the resulting system of linear equations provides the potential values in each x obtained in the spatial discretization and at each frequency. The calculated potential values are depicted in Figure 5.

Figure 5 clearly shows that the sharp potential changes are in the electrodes; in fact, they are the most important factor that determines the change in potential values. It is especially noticeable, that although the electrodes are now modeled with constant functions (Equations (26) and (27)), the steepest changes in potential values occur at the electrode–matter transition, since the length distribution of the impedance is described by a discontinuity function (Figure 4) and the potential functions are continuous (but not differentiable at the discontinuity points). It is interesting to observe that, as a consequence of the spectral behavior of the electrodes (Equations (26) and (27)), at low frequencies, there is a high potential difference between the electrodes and the material under investigation, which decreases significantly with the increasing frequency.

To allow for an easier comparison of the results obtained using the CEM with the continuous case, the magnitude of the potential values is not only presented graphically (Figure 5), but also in tabular form.

Table 1. Typical values of the potential functions $u\left(x\right)$ calculated with the CEM (shown in Figure 5) at mesh points $x=0$, $x=1$, $x=5$, $x=10$, and $x=11$.

$\mathit{u}\mathbf{\left(}\mathit{x}\mathbf{\right)} \mathbf{\left(}\mathit{V}\mathbf{\right)}$	0.1 Hz	10 Hz	25 Hz	63 Hz	158 Hz	398 Hz	1 kHz	10 kHz
$x=0$	1.2026	1.1065	0.9691	0.7878	0.5692	0.3635	0.2269	0.1091
$x=1$	0.7183	0.6516	0.5524	0.4382	0.3138	0.1996	0.1243	0.0595
$x=5$	0.6666	0.6103	0.5308	0.4300	0.3105	0.1984	0.1238	0.0595
$x=10$	0.6054	0.5632	0.5073	0.4211	0.3067	0.1967	0.1232	0.0595
$x=11$	0	0	0	0	0	0	0	0

summarizes the functions shown in Figure 5 with a representation of the most typical values of the analytical solution of Equation (1), calculated with mixed boundary conditions in Equation (19) using Theorem 1.

The potential values shown in Table 1, based on the spatial discretization defined in Equations (14)–(16) were chosen at the following grid points: the Neumann boundary condition ($x=0$), the input electrode and material boundary ($x=1$), the midpoint of the material length ($x=5$), the output electrode and material boundary ($x=10$), and the Dirichlet boundary condition ($x=11$). In Table 1, essentially the same trend is observed as in Figure 5 and, in addition, it is noted that, with increasing frequency, the individual potential profiles are better described with constant functions.

From the modeling of the measurement setup, these main properties can be concluded; however, it is also clear from Figure 5 that the modeling of the electrode–material interface is a sensitive area in the field of EIS. Now, let us consider an extended version of the case study presented so far.

4.3. Continuous Electrode Approach

In the following, consider a case where the measurement setup (in Figure 1) is modeled with the Cole–Cole model defined in Equation (20) and the electrodes are modeled with continuous functions. These functions are able to consider the properties of the electrodes that the CEM is not capable of. This study details the case where the conductivity at the electrode–material contact points is degraded, for example, due to erroneous electrode contact setup.

Now, let the complex conductivity of the electrodes be defined with the following functions:

$${\kappa }_{el,1}^{★}(j\omega ,x)=\left(1-e^{10(x-1.01)}\right){\kappa }_{el,1}(j\omega ,x)x\in [0,1],$$

(29)

and

$${\kappa }_{el,2}^{★}(j\omega ,x)=\left(1-e^{-10(x-9.99)}\right){\kappa }_{el,1}(j\omega ,x)x\in [10,11].$$

(30)

From Equations (29) and (30), it can be seen that the functions ${\kappa }_{el,1}^{★}(j\omega ,x)$ and ${\kappa }_{el,2}^{★}(j\omega ,x)$ are generated from ${\kappa }_{el,1}(j\omega ,x)$ and ${\kappa }_{el,2}(j\omega ,x)$ such that, respectively, by multiplying them with a special space-dependent function, it results in a significant decrease in conductivity in the electrode’s material as it approaches the material under investigation. The piecewise impedance function now takes the following form:

$${\rho }^{★}(j\omega ,x)=\frac{1}{{\kappa }^{★}(j\omega ,x)}=\left\{\begin{array}{cc}\frac{1}{{\kappa }_{el,1}^{★}(j\omega ,x)},\hfill & \mathrm{if}0<x<1,\hfill \\ \frac{1}{{\kappa }_m(j\omega ,x)},\hfill & \mathrm{if}1\le x\le 10,\hfill \\ \frac{1}{{\kappa }_{el,2}^{★}(j\omega ,x)},\hfill & \mathrm{if}10<x<11.\hfill \end{array}\right.$$

(31)

Figure 6 shows the function ${\kappa }^{★}(j\omega ,x)$ at only one frequency, 0.1 Hz, due to the high dynamics of the frequency dependence.

Figure 6 clearly shows that the electrodes represented by Equations (29) and (30) model an extreme measurement artefact. Compared to the previous case in Figure 4, where the electrodes were represented by the concentrated parameter used by the CEM, in the case of Figure 6, the electrodes near the points $x=0$ and $x=L$ also behave as concentrated parameters, since their conductivity is constant. However, now, approaching the investigated material, i.e., $x\to 1^{-}$ and $x\to {10}^{+}$, the conductivity starts to decrease very drastically, hence the applicability of the electrode for EI measurements is seriously reduced. (Naturally, the piecewise function used to model conductivity is always positive on the interval $[0,11]$.) A further remarkable feature of the piecewise functions shown in Figure 6 is that there is a discontinuity at $x=1$ and $x=10$, since the impedance of the electrodes is drastically degraded at the boundary of the material, but this property of the electrodes does not affect the conductivity of the material, hence the conductivity function defined in Section 4.1 is still obtained. Figure 6 shows the piecewise function defined in Equation (31) at a single frequency (0.1 Hz) for clarity. Obviously, varying the frequency generates different functions; however, from Equations (29) and (30), regardless of the frequency at which the electrode conductivity function is considered, a very drastic decrease at the material–electrode interface is always apparent.

The impedance functions $\frac{1}{{\kappa }_{el,1}^{★}}$ and $\frac{1}{{\kappa }_{el,2}^{★}}$ are illustrated in Figure 7.

Compared to Figure 4, while most of the findings are still valid, Figure 7 clearly shows the most important difference: the impedance of the electrodes starts to increase significantly as they approach the material under investigation. This tendency remains true regardless of frequency; however, at lower frequencies, where the electrode impedance becomes dominant, the impedance of the electrodes increases more drastically. Importantly, the curves with the same color in Figure 7 are still the piecewise impedance functions for the same frequency.

Based on Theorem 1, using the exact scheme adapted specifically for the measurement (Equation (12)) and the measurement parameters in Equation (19), the potential values along the length of the material under investigation can still be calculated. The absolute values of the potentials are shown in Figure 8 as a function of the frequency and the x coordinate.

Typical $u\left(x\right)$ potential values in this case can also be presented in tabular form, similar to Table 1.

Table 2. Typical values of the potential functions $u\left(x\right)$ calculated with continuous electrode approach defined in Equation (31) (shown in Figure 8) at mesh points $x=0$, $x=1$, $x=5$, $x=10$, and $x=11$.

$\mathit{u}\mathbf{\left(}\mathit{x}\mathbf{\right)} \mathbf{\left(}\mathit{V}\mathbf{\right)}$	0.1 Hz	10 Hz	25 Hz	63 Hz	158 Hz	398 Hz	1 kHz	10 kHz
$x=0$	1.2449	1.1402	0.9869	0.7947	0.5720	0.3647	0.2273	0.1091
$x=1$	0.7547	0.6806	0.5677	0.4440	0.3162	0.2006	0.1246	0.0595
$x=5$	0.6897	0.6286	0.5404	0.4336	0.3120	0.1990	0.1241	0.0595
$x=10$	0.6128	0.5694	0.5109	0.4225	0.3073	0.1970	0.1233	0.0595
$x=11$	0	0	0	0	0	0	0	0

summarizes the functions shown in Figure 8 with a representation of the $u\left(0\right)$, $u\left(1\right)$, $u\left(5\right)$, $u\left(10\right)$, and $u\left(11\right)$ values of the analytical solution calculated using the continuous electrode approach defined in Equation (31).

Table 2 illustrates clearly (in addition to the properties shown in Figure 8) that, when the frequency increases, the resulting potential profiles are increasingly represented by constant functions. The modeling of the electrodes with a continuous function have provided interesting results; as can be seen in Figure 8, the decrease in conductivity considered for the continuous model visibly modified the potential values measured in the measurement setup. The properties of the functions shown in Figure 8 are almost identical to those in Figure 5; however, the shape of the curves is visibly different from the results calculated using the CEM. The differences between the results are further highlighted in the following discussion by summarizing the results.

4.4. Discussion

The case studies developed in this study present the mathematical modeling of EIS measurements using two different approaches to analyze the modeling of electrodes used in measurements. The investigated material is represented by a location- and frequency-dependent Cole–Cole model, while the electrodes are modeled with two approaches. In the first case study, the well-known and widely used CEM approach is utilized, while in the second, a continuous electrode approach was formulated to exploit exact schemes. The same material was analyzed in all case studies. While for the CEM approach, a location-independent function with constant concentrated parameters was applied, the behavior of the functions used for the continuous electrode approach modelled the conductivity loss of the electrodes.

The calculated potential values show that different potential profiles were obtained for the two electrode approaches. Figure 9 compares the potentials calculated at 0.1 Hz (Figure 9a) and 1 kHz (Figure 9b), respectively.

Figure 8 shows, systematically, that, despite a minimal mathematical change in the applied functions (Equations (26) and (27) vs. Equations (29) and (30)), where the conductivity of the electrodes decreases in the direction of the material, the calculated potential values differ significantly at 0.1 Hz (Figure 9a). The potential profile indicated by the dashed line takes higher values as it approaches the left endpoint. It is interesting to note that, at higher frequencies (1 kHz, Figure 9b), the difference disappears as the electrode impedance becomes insignificant.

Comparing Table 1 and Table 2, it can be seen that the continuous modeling of the electrodes (although only the impedance of the electrode–matter interface was modified) affected the calculated potential values. As a consequence of the properties of the functions used to model the electrodes, the largest difference between the electrode models is observed at the lowest frequency (0.1 Hz). This difference disappears when the frequency increases. This suggests that the differences arising from electrode modeling are significant at frequencies where the electrode properties dominate over the material properties. As a consequence, the inappropriate modeling of the electrodes causes errors in the calculated potential values, which distort the reconstructed impedance profiles in the case of possible impedance measurements.

These results certainly bring a new perspective to the completed research for developing increasingly effective EI measurement and data evaluation procedures. Namely, the methods presented in this paper can be easily implemented in curricular EI-based biological research such as the study of non-alcoholic fatty liver disease [30], the rapid detection of breast cancer in mouse models [31,32], and the development of a new method for in vitro biological studies [33,34].

5. Conclusions

This paper presents a more complete and extended version of combined EI measurements modeling and its advantages. To achieve this, the introduced exact schemes were applied by combining the conductivity functions of the material under investigation and the electrodes into a general piecewise conductivity function. Consequently, the exact scheme of the ODE, defined in this way, produces the values of its analytical solution in the grid points independently of the spatial discretization. Accordingly, a general modeling procedure has been developed, which is capable of implementing a conventional CEM approach; however, it can be efficiently applied to modeling electrodes with continuous functions.

The longitudinal change in the electrode conductivity is dominant in the frequency range where the electrode impedance is comparable with the impedance of the tested material. This justifies the use of continuous functions to model the electrodes used for measurements in terms of space and frequency. The accuracy of the electrode modeling affects the potential values obtained since, in both case studies, the same material is modeled; however, the resulting potentials are different (obviously due to differences in the electrode models). Consequently, when the purpose of the measurement is to reconstruct the impedance profile of the material, two different reconstructed impedance profiles may be obtained for the same material, which is a significant discrepancy. This emphasizes the need for correct and accurate consideration of the electrode artifacts, which minimizes the detailed error effects.

6. Future Works

Consequently, the intention is to apply the results presented in this study to the modeling of EI measurement solutions that have been developed, and to identify several directions for further development.

The first, and perhaps the most important, task is to transfer the findings presented in this paper into practical applications, whether in in vitro or human EI measurements, tomography, or spectroscopy, for each developed prototype. Additionally, the modeling of electrode anomalies is also a major research topic, where the aim is to investigate how these phenomena can be represented as continuous functions and how they can be incorporated into the presented method. Moreover, the aim is to generalize the obtained results to two- and/or three-dimensional cases, in order to extend the relevance of the method to the whole range of EI measurements. Further, research on exact schemes for solving model equations used to describe non-stationary cases for electric field modeling and charge distribution during electrode–material interaction will play a prominent role in future work.