multiphysicsL13.tex

%
% Copyright © 2014 Peeter Joot.  All Rights Reserved.
% Licenced as described in the file LICENSE under the root directory of this GIT repository.
%
% for template copy, run:
%
% ~/bin/ct multiphysicsL1  multiphysicsLectureN tl1
%
%\input{../blogpost.tex}
%\renewcommand{\basename}{multiphysicsL13}
%\renewcommand{\dirname}{notes/ece1254/}
%\newcommand{\keywords}{ECE1254H}
%\input{../peeter_prologue_print2.tex}
%
%\usepackage{kbordermatrix}
%
%\beginArtNoToc
%\generatetitle{ECE1254H Modeling of Multiphysics Systems.  Lecture 13: Continuation parameters and Simulation of dynamical systems.  Taught by Prof.\ Piero Triverio}
%\chapter{Continuation parameters and Simulation of dynamical systems}
\label{chap:multiphysicsL13}
%
%\section{Disclaimer}
%
%Peeter's lecture notes from class.  These may be incoherent and rough.
%
%\chapter{Simulation of dynamical systems}
%\section{Simulation of dynamical systems}
\index{dynamical systems}
%
%Example high level system in \cref{fig:lecture13:lecture13Fig2}.
%
%\imageFigure{../figures/ece1254-multiphysics/lecture13Fig2}{Complex time dependent system.}{fig:lecture13:lecture13Fig2}{0.15}
%
%\section{Assembling equations automatically for dynamical systems.}
\section{Dynamical systems equation assembly.}
%
\index{modified nodal analysis!rc circuit}
\makeexample{RC circuit.}{example:multiphysicsL13:220}{
	The method can be demonstrated well by example.  Consider the RC circuit \cref{fig:lecture13:lecture13Fig3} for which the capacitors \index{capacitor} introduce a time dependence
%
%\imageFigure{../figures/ece1254-multiphysics/lecture13Fig3}{RC circuit}{fig:lecture13:lecture13Fig3}{0.2}
\imageFigure{../figures/ece1254-multiphysics/RC-circuit.pdf}{RC circuit.}{fig:lecture13:lecture13Fig3}{0.2}
%
The unknowns are \( v_1(t), v_2(t) \).
%
The equations (KCLs) at each of the nodes are
\begin{enumerate}
\item
\(
\frac{v_1(t)}{R_1}
+ C_1 \frac{dv_1}{dt}
+ \frac{v_1(t) - v_2(t)}{R_2}
+ C_2 \frac{d(v_1 - v_2)}{dt}
- i_{s,1}(t) = 0
\)
\item
\(
\frac{v_2(t) - v_1(t)}{R_2}
+ C_2 \frac{d(v_2 - v_1)}{dt}
+ \frac{v_2(t)}{R_3}
+ C_3 \frac{dv_2}{dt}
- i_{s,2}(t)
= 0
\)
\end{enumerate}
%
This has the matrix form
\begin{equation}\label{eqn:multiphysicsL13:140}
\begin{aligned}
&\begin{bmatrix}
Z_1 + Z_2 & - Z_2 \\
-Z_2 & Z_2 + Z_3
\end{bmatrix}
\begin{bmatrix}
v_1(t) \\
v_2(t)
\end{bmatrix} \\
&\qquad +
\begin{bmatrix}
C_1 + C_2 & -C_2 \\
- C_2 & C_2 + C_3
\end{bmatrix}
\begin{bmatrix}
\frac{dv_1(t)}{dt} \\
\frac{dv_2(t}{dt})
\end{bmatrix} 
=
\begin{bmatrix}
1 & 0 \\
0 & 1
\end{bmatrix}
\begin{bmatrix}
i_{s,1}(t) \\
i_{s,2}(t)
\end{bmatrix}.
\end{aligned}
\end{equation}
%
%
Observe that the capacitor between node 2 and 1 is associated with a stamp of the form
%
\begin{equation}\label{eqn:multiphysicsL13:180}
\begin{bmatrix}
C_2 & -C_2 \\
-C_2 & C_2
\end{bmatrix},
\end{equation}
%
very much like the impedance stamps of the resistor node elements.
}
%
The RC circuit problem has the abstract form
%
\begin{equation}\label{eqn:multiphysicsL13:160}
\BG \Bx(t) + \BC \frac{d\Bx(t)}{dt} = \BB \Bu(t),
\end{equation}
%
which is more general than a \textAndIndex{state space} equation of the form
%
\begin{equation}\label{eqn:multiphysicsL13:200}
\frac{d\Bx(t)}{dt} = \BA \Bx(t) + \BB \Bu(t).
\end{equation}
%
Such a system may be represented diagrammatically as in \cref{fig:lecture13:lecture13Fig4}.
%
\imageFigure{../figures/ece1254-multiphysics/lecture13Fig4}{State space system.}{fig:lecture13:lecture13Fig4}{0.1}
%
The \( C \) factor in this capacitance system, is generally not invertible.  For example, if consider a 10 node system with only one capacitor, for which \( C \) will be mostly zeros.
In a state space system, in all equations have a derivative.  All equations are dynamical.
%
The time dependent MNA \index{modified nodal analysis} system for the RC circuit above, contains a mix of dynamical and algebraic equations.  This could, for example, be a pair of equations like
%
\begin{subequations}
\begin{equation}\label{eqn:multiphysicsL13:240}
   \frac{dx_1}{dt} + x_2 + 3 = 0
\end{equation}
\begin{equation}\label{eqn:multiphysicsL13:260}
   x_1 + x_2 + 3 = 0
\end{equation}
\end{subequations}
%
\paragraph{How to handle inductors}
%\index{Modified nodal analysis!inductor}
\index{modified nodal analysis!inductor}
%
A pair of nodes that contains an inductor element, as in \cref{fig:lecture13:lecture13Fig5}, has to be handled specially.
%
\imageFigure{../figures/ece1254-multiphysics/lecture13Fig5}{Inductor configuration.}{fig:lecture13:lecture13Fig5}{0.1}
%
The KCL at node 1 has the form
%
\begin{equation}\label{eqn:multiphysicsL13:280}
\cdots + i_{\textrm{L}}(t) + \cdots = 0,
\end{equation}
%
where
%
\begin{equation}\label{eqn:multiphysicsL13:300}
v_{n_1}(t) - v_{n_2}(t) = L \frac{d i_{\textrm{L}}}{dt}.
\end{equation}
%
It is possible to express this in terms of \( i_{\textrm{L}} \), the variable of interest
%
\begin{equation}\label{eqn:multiphysicsL13:320}
i_{\textrm{L}}(t) = \inv{L} \int_0^t \lr{ v_{n_1}(\tau) - v_{n_2}(\tau) } d\tau
+ i_{\textrm{L}}(0).
\end{equation}
%
Expressing the problem directly in terms of such integrals makes the problem harder to solve, since the usual differential equation toolbox cannot be used directly.  An integro-differential toolbox would have to be developed.  What can be done instead is to introduce an additional unknown for each inductor current derivative \( di_{\textrm{L}}/dt \), for which an additional MNA row is introduced for that inductor scaled voltage difference.
%
\section{Numerical solution of differential equations.}
%
Consider the one variable system
%
\begin{equation}\label{eqn:multiphysicsL13:340}
G x(t) + C \frac{dx}{dt} = B u(t),
\end{equation}
%
given an initial condition \( x(0) = x_0 \).  Imagine that this system has the solution sketched in \cref{fig:lecture13:lecture13Fig6}.
%
\imageFigure{../figures/ece1254-multiphysics/lecture13Fig6}{Discrete time sampling.}{fig:lecture13:lecture13Fig6}{0.2}
%
Very roughly, the steps for solution are of the form
%
\begin{enumerate}
\item Discretize time
\item Aim to find the solution at \( t_1, t_2, t_3, \cdots \)
\item Use a finite difference formula to approximate the derivative.
\end{enumerate}
%
There are various schemes that can be used to discretize, and compute the finite differences.
%
\section{Forward Euler method.}
\index{Euler!forward method}
%
One such scheme is to use the forward differences, as in \cref{fig:lecture13:lecture13Fig7}, to approximate the derivative
%
\begin{equation}\label{eqn:multiphysicsL13:360}
\dot{x}(t_n) \approx \frac{x_{n+1} - x_n}{\Delta t}.
\end{equation}
%
\imageFigure{../figures/ece1254-multiphysics/lecture13Fig7}{Forward difference derivative approximation.}{fig:lecture13:lecture13Fig7}{0.3}
%
Introducing this into \cref{eqn:multiphysicsL13:340} gives
%
\begin{equation}\label{eqn:multiphysicsL13:350}
G x_n + C \frac{x_{n+1} - x_n}{\Delta t} = B u(t_n),
\end{equation}
%
or
%
\begin{equation}\label{eqn:multiphysicsL13:380}
C x_{n+1} = \Delta t B u(t_n) - \Delta t G x_n + C x_n.
\end{equation}
%
The coefficient \( C \) must be invertible, and the next point follows immediately
%
%\begin{equation}\label{eqn:multiphysicsL13:400}
\boxedEquation{eqn:multiphysicsL13:400}{
x_{n+1} = \frac{\Delta t B}{C} u(t_n) + x_n \lr{ 1  - \frac{\Delta t G}{C} }
}
%\end{equation}
%
%\EndArticle
%\EndNoBibArticle