\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 Dirichlet boundary condition\\
\\
\begin{equation}\phi  = p\end{equation}
\\
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}

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 (0 V), 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.\\

\subsection{Asymptotic Boundary Conditions}
For 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}
\\
\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}
\\