Update theory manual and user guide

Chris Smartt
1 parent 73f9f473
Showing 6 changed files with 143 additions and 7 deletions Show diff stats
DOCUMENTATION/SACAMOS_Spacecraft_cable_MOD_development.pdf
DOCUMENTATION/SACAMOS_TheoryManual.pdf
DOCUMENTATION/SACAMOS_UserGuide.pdf
DOCUMENTATION/THEORY_MANUAL/Tex/Laplace_Solver_Documentation.tex
DOCUMENTATION/THEORY_MANUAL/Tex/bibliography.tex
DOCUMENTATION/USER_GUIDE/Tex/creating_a_cable_bundle_model.tex
@@ -61,10 +61,8 @@ In the static field problems that need to be solved in order to determine the en
 \\
 \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}
-\\
+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}
@@ -97,7 +95,7 @@ We can also immediately compute the coefficients of capacitance ${C_{{\rm{pq}}}}
 \\
 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}
+\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}
@@ -108,9 +106,27 @@ The procedure for determining the capacitance matrix via computation of the ener
 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.\\
+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 condcutor 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 derivitive of the potential is zero and hence the charge on the boundary is zero everywhere. 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}
  
-\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}
@@ -131,6 +147,11 @@ Expression  immediately implies that the potential for $r \gg r&#39;$, with $r&#39;$ the
 \\
 \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 a net charge on the outer boundary which is not taken into account in 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)\\
 \\
@@ -165,4 +186,112 @@ From that we can obtain the capacitance matrix $\hat C$\\
 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 does not take advantage of this sparsity and a direct matrix inverse is 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}.
+
+CCN2 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
  
@@ -39,6 +39,11 @@ J-M ~Jin,
 \textsl{The Finite Element Method in Electromagnetics} 
 3rd Edition, John Wiley \& sons 2014.
  
+\bibitem{numerical_recipes} 
+W.H.~Press, S.A. Teukolsky, W. T. Vetterling, B. P. Flannery
+\textsl{Numerical Recipes in C} 
+Second Edition, Cambridge University Press, 1992.
+
 \bibitem{Chen}
 Q. Chen, A. Konrad, 
 \textsl{A Review of Finite Element Open Boundary Techniques for Static and Quasi-Static Electromagnetic Field Problems}, 
@@ -92,6 +92,8 @@ If the Laplace solver is used then the mesh generation is be controlled by the p
  
 \item `Max\_mesh\_edge\_length'  This (real) parameter determines the maximum finite element mesh edge length for the rectangular conductors of flex cables. It may be used to force an increase in mesh density in flex cables. \\
  
+\item `gp\_edge\_length'  This (real) parameter determines the finite element mesh edge length at the centre of the ground plane. It may be used to increase the mesh density in the ground plane below the cable bundle to improve accuracy of the per-unit-length parameter calculation. \\
+                                   
 \end{enumerate}
  
 The default parameters are a compromise between accuracy and computation time for the Laplace solution.