\chapter{Numerical computation of the per-unit-length MTL parameters using a FEM solver} \label{Laplace}
Multiconductor Transmission Line (MTL) modelling is based on a set of coupled first order phasor MTL equations as discussed in chapter \ref{MTLtheory}. 
In the set of coupled first-order MTL equations the impedance matrix $\left[\tilde{Z}\left( \omega  \right)\right]$ and admittance matrix $\left[\tilde{Y}\left( \omega  \right)\right]$ are given by\\
\\

\begin{equation} \label{eq:Admittance}
\left[\tilde{Z}\left( \omega  \right)\right] = \left[\tilde{R}\left( \omega  \right)\right] + j\omega {\kern 1pt} \left[\tilde{L}\left( \omega  \right)\right]
\end{equation}

\begin{equation}
\left[\tilde{Y}\left( \omega  \right)\right] = \left[\tilde{G}\left( \omega  \right)\right] + j\omega {\kern 1pt} \left[\tilde{C}\left( \omega  \right)\right]
\end{equation}
%
with $\left[\tilde{R}\left( \omega  \right)\right]$ the per unit length resistance matrix, $\left[\tilde{L}\left( \omega  \right)\right]$ the per unit length inductance matrix, $\left[\tilde{G}\left( \omega  \right)\right]$ the per unit length conductance matrix and $\left[\tilde{C}\left( \omega  \right)\right]$the per unit length capacitance matrix. In order to be able to perform a propagation or crosstalk analysis using MTL modelling it is important to determine good approximations of these per-unit length parameter matrices. \\

It is assumed in the initial discussion that the inductance and capacitance matrix are independent of frequency and that dielectrics are loss-less and frequency independent. In section \ref{conductance_matrix} frequency dependent capacitance is introduced by means of a complex relative permittivity. In this case the capacitance and conductance matrix may be evaluated at specific frequencies.
\\
In the book by Paul,\cite{PaulMTL}, chapter 5, a definition of the $\left[L\right]$, $\left[G\right]$ and $\left[C\right]$ matrices is given and it is explained that these matrices can be determined from static field solutions in the transverse plane for perfect line conductors. Closed-form solutions for these static field problems however only exist for a small set of configurations and several approximate solutions have limited applicability. As a consequence approximate numerical methods must be used for the computation of the entries in the per-unit-length matrices. In the same book by Paul, chapter 5, it is also shown that in order to numerically determine the entries of the $\left[L\right]$, $\left[G\right]$ and $\left[C\right]$ matrices we only need a capacitance solver, that is, we only need a numerical solver that is able to compute a capacitance matrix $\left[C\right]$. The inductance matrix $\left[L\right]$ can be derived from the capacitance matrix solution ${\left[C_0\right]}$ for the static field problem in which we have replaced all dielectrics in the cross -section with free space, the per-unit-length inductance matrix can be obtained from\\
\\
\begin{equation}\left[L\right] = {\mu _0}{\varepsilon _0}\left[C_0\right]^{ - 1}\end{equation}
\\
with ${\mu _0}$ and ${\varepsilon _0}$ the permeability and permittivity of free space. The conductance matrix $\left[G\right]$ can be obtained by introducing a complex permittivity and solving for a complex valued capacitance matrix $\left[C\right]$. \\
\\
Two dimensional electrostatic field problems are governed by the two-dimensional Poisson equation\\
\\
\begin{equation} - \frac{\partial }{{\partial x}}\left( {{\varepsilon _{{\rm{r,x}}}}\frac{{\partial \phi }}{{\partial x}}} \right) - \frac{\partial }{{\partial y}}\left( {{\varepsilon _{{\rm{r,y}}}}\frac{{\partial \phi }}{{\partial y}}} \right) = \frac{{{\rho _{\rm{c}}}}}{{{\varepsilon _{\rm{0}}}}}\end{equation}
\\
with $\phi $ the potential, ${\rho _{\rm{c}}}$ the charge density, ${\varepsilon _{\rm{0}}}$ the permittivity of free space, ${\varepsilon _{{\rm{r,x}}}}$ and ${\varepsilon _{{\rm{r,y}}}}$ the diagonal components of the relative permittivity tensor (we have assumed all off-diagonal terms are zero) of the materials in the configuration. \\
\\
We can solve the capacitance matrix, $\left[C\right]$, by discretising the cross-section of the cable (MTL) and solving the Laplace equation for this cross section numerically.\\
\\
The two-dimensional Poisson equation  follows from Gauss's law, \\
\\
\begin{equation}\frac{\partial }{{\partial x}}{D_x} + \frac{\partial }{{\partial y}}{D_y} = {\rho _{\rm{c}}}\end{equation}
\\
the constitutive relation between the electric flux density and the electric field intensity for a linear axial medium \\
\\
\begin{equation}\left[ {\begin{array}{*{20}{c}}\\
   {{D_x}}  \\\\
   {{D_y}}  \\\\
\end{array}} \right] = {\varepsilon _0}\left[ {\begin{array}{*{20}{c}}\\
   {{\varepsilon _{{\rm{r}},x}}} & 0  \\\\
   0 & {{\varepsilon _{{\rm{r}},y}}}  \\\\
\end{array}} \right]\left[ {\begin{array}{*{20}{c}}\\
   {{E_x}}  \\\\
   {{E_y}}  \\\\
\end{array}} \right]\end{equation}
\\
and the definition of the scalar electric potential:\\
\\
\begin{equation}\left[ {\begin{array}{*{20}{c}}\\
   {{E_x}}  \\\\
   {{E_y}}  \\\\
\end{array}} \right] =  - \left[ {\begin{array}{*{20}{c}}\\
   {\frac{{\partial \phi }}{{\partial x}}}  \\\\
   {\frac{{\partial \phi }}{{\partial y}}}  \\\\
\end{array}} \right]\end{equation}
\\
In the static field problems that need to be solved in order to determine the entries of the capacitance matrix of a MTL, there is no free charge in the region between the conductors and the potential is prescribed on the conductors. As a result we have to solve the two dimensional field problem:\\
\\
\begin{equation} - \frac{\partial }{{\partial x}}\left( {{\varepsilon _{\rm{r}}}\frac{{\partial \phi }}{{\partial x}}} \right) - \frac{\partial }{{\partial y}}\left( {{\varepsilon _{\rm{r}}}\frac{{\partial \phi }}{{\partial y}}} \right) = 0\end{equation}
\\

with an appropriate boundary condition applied on the outer boundary of the domain.

In the SACAMOS project the capacitance (and conductance) matrices are calculated numerically using a Finite Element Method \cite{Jin} as outlined in the rest of this chapter. Details of the implementation are not discussed here however the software closely follows the process described in chapter 4 of \cite{Jin}.
\\
\section{Solving the capacitance matrix based on energy}
The capacitance matrix can be solved based on energy or based on charges. We will use the method based on energy.\\
The total electrostatic energy stored in the multi-conductor configuration can be computed from the static field solution as\\
\\
\begin{equation} \label{eq:EnergyStored}
W=\frac{1}{2}\iint_S\mathbf{D}\mathbf{E}dA		\\
\end{equation} 

For this method we set the pth and qth conductor to 1V, that is, ${V_{\rm{p}}} = {V_{\rm{q}}} = 1$ and connect all other conductors to reference, that is, all other conductors are set to 0V which implies${V_{\rm{i}}} = 0$ for all ${\rm{i}} \ne {\rm{q~or~p}}$, the electrostatic energy ${W_{{\rm{pq}}}}$stored in the system equals\\
\\
\begin{equation}{W_{{\rm{pq}}}} = {W_{\rm{p}}} + {W_{\rm{q}}} - {c_{{\rm{pq}}}}\end{equation}
\\
where $c_{\rm{pq}}$ is the coupling capacitance. 
As a result, in case we have determined ${W_{\rm{p}}}$ for all $p = 1, \cdots ,n$ and ${W_{\rm{pq}}}$ for all combinations of $p = 1, \cdots ,n$ and  $q = 1, \cdots ,n$ for which $p \ne q$ we can first compute all coupling capacitances ${c_{{\rm{pq}}}}$, namely\\
\\
\begin{equation}{c_{{\rm{pq}}}} = {W_{\rm{p}}} + {W_{\rm{q}}} - {W_{{\rm{pq}}}}\end{equation}
\\
Then we may compute the coupling capacitances ${c_{{\rm{p0}}}}$\\
\\
\begin{equation}{c_{{\rm{p0}}}} = 2{W_{\rm{p}}} - \sum\limits_{\scriptstyle j = 1 \hfill \atop \\
  \scriptstyle j \ne p \hfill}^n {\,{c_{{\rm{pj}}}}} \end{equation}
\\
We can also immediately compute the coefficients of capacitance ${C_{{\rm{pq}}}}$ for all combinations of $p = 1, \cdots ,n$ and  $q = 1, \cdots ,n$ as we note that:\\
\\
\begin{equation} \label{eq:Cpq} {C_{{\rm{pq}}}} =  - {c_{{\rm{pq}}}} = {W_{{\rm{pq}}}} - {W_{\rm{p}}} - {W_{\rm{q}}}\end{equation}
\begin{equation} \label{eq:Cpp} {C_{{\rm{pp}}}} = \sum\limits_{\scriptstyle j = 0 \hfill \atop \\
  \scriptstyle j \ne p \hfill}^n {{c_{{\rm{pj}}}}}  = 2{W_{\rm{p}}}\end{equation}
\\
So, with a static field solver available we are able to compute the entries in the capacitance matrix.\\
\\

\section{Procedure capacitance matrix computation based on energy analysis}\label{LC_energy}

The procedure for determining the capacitance matrix via computation of the energy storage is:\\
\begin{enumerate}
\item First select all $n$ non-reference conductors one after another and each time compute the static field solution with the voltage of the conductor under consideration set to 1V and all other conductors set to the reference 0V. With the pth conductor set to 1V (${V_{\rm{p}}} = 1$ V) and all other conductors to reference (0 V), compute the electrostatic energy ${E_{\rm{p}}}$stored in the configuration by using the surface integral \eqref{eq:EnergyStored}. This requires that the static field problem and energy integral has to be determined a total of $n$ times. \\

\item Second select all combinations of two different conductors from the $n$ non-reference conductors one after another and each time compute the static field solution with the voltage of the two conductors under consideration set to 1V and all other conductors set to the reference 0V. With the p$^{\mbox{th}}$ and q$^{\mbox{th}}$ conductor set to 1V (${V_{\rm{p}}} = {V_{\rm{q}}} = 1$) and all other conductors set to reference (0V), compute the electrostatic energy ${E_{{\rm{pq}}}}$stored in the configuration by using the surface integral \eqref{eq:EnergyStored}. This requires that the static field solution and energy integral has to be determined a total of $n\left( {n - 1} \right)/2$ times. \\

\item Determine the entries of the capacitance matrix from the computed stored energies ${E_{\rm{p}}}$and ${E_{{\rm{pq}}}}$ using expressions \eqref{eq:Cpq} and \eqref{eq:Cpp}. \\
\end{enumerate}
The reference conductor can be either a ground plane, a shield or a circular or rectangular conductor.\\

\subsection{FEM Solver for the Open Boundary Electrostatic Field Problem}

Without an outer shield the electrostatic field is no longer confined to a bounded domain while the FEM can only be used on a finite computational domain. In the course of time a number of techniques have been developed to circumvent this problem. A review of the most important finite element open boundary techniques is given by Chen and Konrad \cite{Chen}. 

%Of all available techniques we have chosen to truncate the domain by introducing a fictitious contour that encloses the domain of primary %interest and prescribe on it an asymptotic boundary condition which mimics the behaviour of the electrostatic potential in the region %outside the contour. Not only has it been shown that asymptotic boundaries give good results, the technique is relatively easy to implement %and preserves the sparsity of the finite element coefficient matrix which makes it a very efficient method as-well.\\

In unshielded conductor configurations a suitable boundary condition must be applied on the outer boundary of the problems space. The imposition of a boundary condition will necessarily have an effect on the fields within the Finite Element domain and hence the capacitance matrix to be calculated. Three possible boundary conditions are outlined here, the Dirichelet boundary condition ($\phi=0$), the Neumann boundary condition $\frac{\partial \phi(r) }{\partial r}=0$ and an asymptotic boundary condition,  

$\frac{\partial \phi(r) }{\partial r}=\frac{1}{r ln(r)}\phi(r)$. 

\subsubsection{Dirichelet boundary conditions}

The Dirichelet boundary condition imposes a $\phi=0$ condition on the outer boundary. This implies that the outer boundary is another conductor which may store charge. In order to calcualte the capacitance matrix with the Dirichelet boundary condition, the outer boundary must be included in the analysis as an additional conductor and a 'generalised capacitance' matrix calculated. The capacitance matrix for the interior conductors only may then be calculated from the generalised capacitance matrix via the procedure described in \cite{PaulMTL}.

\subsubsection{Neumann boundary conditions}

The Neumann boundary condition is expressed as
\begin{equation} 
\frac{\partial \phi(r) }{\partial r}=0
\end{equation}

This boundary condition enforces a condition that the normal derivative of the potential is zero and hence the charge inside the boundary is zero (by application of the integral form of Gauss' law). The zero charge propertiy is required by the capacitance matrix calculation and using the Neumann boundary condition leads to convergence of the capacitance matrix elements to the correct values as the boundary distance and the mesh density are increased.

\subsubsection{Asymptotic boundary conditions}

For 2D electrostatic problems, the potential in a point ${\bf{r}}$, can be written as\\

\\
\begin{equation}\phi\left(\mathbf{r}\right)=\iint_{\Omega}\rho\left(\mathbf{r}'\right) G\left(\mathbf{r},\mathbf{r}'\right)d\mathbf{r}'\end{equation}
\\
Where $\Omega$ is the charged region, $\rho \left( {{\bf{r}}'} \right)$ is the charge distribution and $G\left( {{\bf{r}},{\bf{r}}'} \right)$ is the free space Green’s function:\\
\\
\begin{equation}G\left( {{\bf{r}},{\bf{r}}'} \right) = \frac{1}{{2\pi {\varepsilon _0}}}\ln \left| {{\bf{r}} - {\bf{r}}'} \right|\end{equation}
\\
In the expression for the potential , the effect of dielectrics is taken into account by incorporating the polarization charge density in the charge distribution. Using polar coordinates $(r,\theta )$ the position vectors may be written as ${\bf{r}} = \left( {r\cos \theta ,r\sin \theta } \right)$ and ${\bf{r}}' = \left( {r'\cos \theta ',r'\sin \theta '} \right)$ and the free space Green function equals\\
\\
\begin{equation}G\left( {r,\theta ,r',\theta '} \right) = \frac{1}{{4\pi {\varepsilon _0}}}\ln \left( {{r^2} + r{'^2} - 2rr'\cos \left( {\theta  - \theta '} \right)} \right)\end{equation}
\\
and hence\\
\\
\begin{equation}\frac{\partial }{{\partial r}}G\left( {r,\theta ,r',\theta '} \right) = \frac{1}{{2\pi {\varepsilon _0}}}\frac{{r - r'\cos \left( {\theta  - \theta '} \right)}}{{{r^2} + r{'^2} - 2rr'\cos \left( {\theta  - \theta '} \right)}}\end{equation}
\\
Expression  immediately implies that the potential for $r \gg r'$, with $r'$ the distances from the origin of the charges in the charged region, satisfies \\
\\
\begin{equation}\frac{\partial }{{\partial r}}\phi \left( {r,\theta } \right) = \frac{1}{{r\ln \left( r \right)}}\phi \left( {r,\theta } \right)\end{equation}
\\

It has been observed that the use of the asymptotic boundary condition described here can have an effect on the convergence of the capacitance matrix calculation due to the total charge within the boundary not being null, which impacts the capacitance matrix calculation described in section \ref{LC_energy}. Currently it is not known whether a 'generalised capacitance matrix' type solution can be used along with this boundary condition as the outer boundary is not an equipotential (as it is in the Dirichelet case).

As a consequence of this the Neumann boundary condition is recommended and used in the SACAMOS software as the default open boundary condition.

\subsection{Computing the conductance matrix} \label{conductance_matrix}
The conductance matrix $\left[G\right]$ can be obtained by introducing a complex permittivity and solving for a complex valued capacitance matrix $\left[C\right]$. ( reference \cite{PaulMTL} equation 3.6.1)\\
\\
If the surrounding medium is homogeneous, G can be obtained from C or L by\\
\\
\begin{equation}G = \frac{\sigma }{\varepsilon }C = \mu \sigma {L_e}^{ - 1}\end{equation}
${L_e}$ is the external inductance matrix assuming perfect conductors. \\
Next we want to obtain G for inhomogeneous media. This can be the case for a conductor or shield that has a dielectric jacket (insulator). Actually, this is a simple modification of the calculation of the per-unit-length capacitance matrix C (also computed assuming perfect conductors). Losses in the medium are due to:\\
\begin{itemize}
\item Conductive losses (due to the conductance parameter $\sigma $).\\
\item Polarization losses\\
\end{itemize}
Both of the losses can be represented by a complex permittivity.\\
\\
	\begin{equation}\hat \varepsilon  = \varepsilon ' - j\varepsilon ''\end{equation}
\\
	\begin{equation}\varepsilon ' = \varepsilon \end{equation}
\\
	\begin{equation}\varepsilon '' = {\varepsilon _b} + \frac{\sigma }{\omega } = \varepsilon '\tan (\delta )\end{equation}
\\
where $\tan (\delta )$ is the loss tangent and ${\varepsilon _b}$ are the polarization losses due to bound charges. $\sigma $ represents the conduction current losses.\\
We can include losses in the medium by solving for the capacitance matrix for a medium having a complex permittivity 
\begin{equation}\hat \varepsilon  = \varepsilon'\left(1 - j\tan \left(\delta\right)\right) 
\end{equation}
and from that result directly obtain C and G.\\
	\begin{equation}\hat Y = G + j\omega C\end{equation}
	\begin{equation}\hat Y = j\omega \left(C - \frac{j}{\omega }G\right) = j\omega \hat C\end{equation}
\\
From that we can obtain the capacitance matrix $\hat C$\\
	\begin{equation}C = {\mathop{\rm Re}\nolimits} (\hat C)\end{equation}
\\
And the conductance matrix G\\
	\begin{equation}G =  - \varpi {\mathop{\rm Im}\nolimits} (\hat C)\end{equation}
\\

\subsection{Iterative solution of the Finite Element matrix equations} \label{cg_solution}

The Finite Element method used to solver the Laplace equation for multi-conductor cable cross sections results in a matrix equation which must be solved for the unknown node potentials. The matrix is sparse however the V3 version of the SACAMOS software did not take advantage of this sparsity and a direct matrix inverse was used to solve the matrix equation. This approach is wasteful in that the sparsity of the matrix equation is not used to our advantage and the size of the Finite Element mesh which can be used is severely limited to systems of the order of 1000's of unknowns. Iterative matrix solutions have the potential to improve the efficiency of solution in terms of both storage (a full matrix never needs to be stored) and runtime for a solution. The method implemented here is described in reference \cite{numerical_recipes}. A discussion of the use of conjugate gradient methods in Finite Element solutions may be found in reference \cite{Jin}.

The V4 version of the SACAMOS sooftware addresses this issue by implementing the biconjugate gradient method, an iterative approach to the solution to the Finite Element matrix equations which is suitable for real and complex symmetric matrix systems as in equation \ref{eq:matrix_eqn}.

\begin{equation} \label{eq:matrix_eqn}
\left[ A \right] \left( x \right)=\left( b \right)
\end{equation}

% Include references to Jin's book and Numerical Recipes

The method is derived by writing a functional whose minimum corresponds to the solution of \ref{eq:matrix_eqn}. Here we minimise the function

\begin{equation} \label{eq:matrix_eqn}
f(x)=\frac{1}{2}\left( x \right) \left[ A \right] \left( x \right)-\left( b \right)\left( x \right)
\end{equation}

This function is minimised when it's gradient is zero i.e.

\begin{equation} \label{eq:matrix_eqn}
\nabla f=\left[ A \right] \left( x \right)-\left( b \right)=0
\end{equation}


The inner product is defined as follows: if we have two column vectors $(a)$ and $(b)$ each with n elements then

\begin{equation} \label{eq:inner_product}
\left< a,b \right>=\sum_{i=1}^{n} a_i b_i
\end{equation}

The bicongugate gradient algorithm is described by the following process:

given an initial guess of the solution vector $(x_1)$, calculate the residual

\begin{equation} \label{eq:BCG_1}
\left( r_1 \right)=\left( b \right)-\left[ A \right] \left( x_1 \right)
\end{equation}
then set 
\begin{equation} \label{eq:BCG_2}
\left( r'_1 \right)=\left[ A \right] \left( r_1 \right)
\end{equation}
This choice of $\left( r'_1 \right)$ gives the so called 'minimum residual' algorithm which is suitable for
non-positive definite matrices. 

Also set the initial values of $\left( p_1 \right)$ and $\left( p'_1 \right)$ as
\begin{equation} \label{eq:BCG_3}
\left( p_1 \right)=\left( r_1 \right)
\end{equation}
\begin{equation} \label{eq:BCG_4}
\left( p'_1 \right)=\left( r'_1 \right)
\end{equation}

Then iterate the following process for $k=1,2,3 \hdots$ until the size of the residual is less than a specified tolerance $(10^{-12})$

\begin{equation} \label{eq:BCG_5}
a_k=\frac{\left< \left(r'_k\right),\left(r_k\right) \right>}{\left< \left(p'_k\right)\left[ A \right] \left(p_k\right) \right> }
\end{equation}

The improved estimate of the unknown solution vector is now
\begin{equation} \label{eq:BCG_6}
\left(x_{k+1}\right)=\left(x_k\right)+a_k \left(p_k\right)
\end{equation}

The residual is given by

\begin{equation} \label{eq:BCG_7}
\left( r_k \right)=\left( b \right)-\left[ A \right] \left( x_k \right)
\end{equation}

Calculate the size (norm) of the residual
\begin{equation} \label{eq:BCG_8}
norm=\left< r_k,r^*_k \right>=\sum_{i=1}^{n} r_{k,i} r^*_{k,i}
\end{equation}

If norm is less than the specified tolerance then finish, otherwise calculate 

\begin{equation} \label{eq:BCG_9}
\left( r_{k+1} \right)=\left( r_k \right)-a_k\left[ A \right] \left( p_k \right)
\end{equation}
   

\begin{equation} \label{eq:BCG_10}
\left( r'_{k+1} \right)=\left( r'_k \right)-a_k\left[ A^T \right] \left( p'_k \right)
\end{equation}

\begin{equation} \label{eq:BCG_11}
b_k=\frac{\left< \left(r'_{k+1}\right),\left(r_{k+1}\right) \right>}{\left< \left(r'_k\right),\left(r_k\right) \right> }
\end{equation}
   

\begin{equation} \label{eq:BCG_12}
\left( p_{k+1} \right)=\left( r_{k+1} \right)+b_k\left( p_k \right)
\end{equation}

\begin{equation} \label{eq:BCG_12}
\left( p'_{k+1} \right)=\left( r'_{k+1} \right)+b_k\left( p'_k \right)
\end{equation}
        
Then return to the start of the iteration loop  

Two implementations of the above iterative process have been developed, one for real matrices and one for complex matrices. 
The complex solution is only required if lossy dielectric materials are present in the cross section.  



\clearpage