ECE1254H Modeling of Multiphysics Systems. Lecture 4: Modified nodal analysis and solutions of large systems. Taught by Prof. Piero Triverio

September 26, 2014 ece1254 , , , ,

Disclaimer

Peeter’s lecture notes from class. These may be incoherent and rough.

Modified nodal analysis

• branch currents for voltage sources
• all elements for which it is impossible or inconvenient to write $$i = f(v_1, v_2)$$.

Imagine, for example, that we have a component illustrated in fig. 1.

fig. 1. Variable voltage device

\label{eqn:multiphysicsL4:20}
v_1 – v_2 = 3 i^2

• any current which is controlling dependent sources, as in fig. 2.

fig. 2. Current controlled device

• Inductors
\label{eqn:multiphysicsL4:40}
v_1 – v_2 = L \frac{di}{dt}.

Solving large systems

We are interested in solving linear systems of the form

\label{eqn:multiphysicsL4:60}
M \overline{{x}} = \overline{{b}},

possibly with thousands of elements.

Gaussian elimination

\label{eqn:multiphysicsL4:80}
\begin{bmatrix}
M_{11} & M_{12} & M_{13} \\
M_{21} & M_{22} & M_{23} \\
M_{31} & M_{32} & M_{33} \\
\end{bmatrix}
\begin{bmatrix}
x_1 \\
x_2 \\
x_3 \\
\end{bmatrix}
=
\begin{bmatrix}
b_1 \\
b_2 \\
b_3 \\
\end{bmatrix}

It’s claimed for now, to be seen later, that back substitution is the fastest way to arrive at the solution, less computationally complex than completion the diagonalization.

Steps

\label{eqn:multiphysicsL4:100}
(1) \cdot \frac{M_{21}}{M_{11}} \Rightarrow
\begin{bmatrix}
M_{21} & \frac{M_{21}}{M_{11}} M_{12} & \frac{M_{21}}{M_{11}} M_{13} \\
\end{bmatrix}

\label{eqn:multiphysicsL4:120}
(2) \cdot \frac{M_{31}}{M_{11}} \Rightarrow
\begin{bmatrix}
M_{31} & \frac{M_{31}}{M_{11}} M_{32} & \frac{M_{31}}{M_{11}} M_{33} \\
\end{bmatrix}

This gives

\label{eqn:multiphysicsL4:140}
\begin{bmatrix}
M_{11} & M_{12} & M_{13} \\
0 & M_{22} – \frac{M_{21}}{ M_{11} } M_{12} & M_{23} – \frac{M_{21}}{M_{11}} M_{13} \\
0 & M_{32} – \frac{M_{31}}{M_{11}} M_{32} & M_{33} – \frac{M_{31}}{M_{11}} M_{33} \\
\end{bmatrix}
\begin{bmatrix}
x_1 \\
x_2 \\
x_3 \\
\end{bmatrix}
=
\begin{bmatrix}
b_1 \\
b_2 – \frac{M_{21}}{M_{11}} b_1 \\
b_3 – \frac{M_{31}}{M_{11}} b_1
\end{bmatrix}.

Here the $$M_{11}$$ element is called the {pivot}. Each of the $$M_{j1}/M_{11}$$ elements is called a {multiplier}. This operation can be written as

\label{eqn:multiphysicsL4:160}
\begin{bmatrix}
M_{11} & M_{12} & M_{13} \\
0 & M_{22}^{(2)} & M_{23}^{(3)} \\
0 & M_{32}^{(2)} & M_{33}^{(3)} \\
\end{bmatrix}
\begin{bmatrix}
x_1 \\
x_2 \\
x_3 \\
\end{bmatrix}
=
\begin{bmatrix}
b_1 \\
b_2^{(2)} \\
b_3^{(2)} \\
\end{bmatrix}.

Using $$M_{22}^{(2)}$$ as the pivot this time, we form
\label{eqn:multiphysicsL4:180}
\begin{bmatrix}
M_{11} & M_{12} & M_{13} \\
0 & M_{22}^{(2)} & M_{23}^{(3)} \\
0 & 0 & M_{33}^{(3)} – \frac{M_{32}^{(2)}}{M_{22}^{(2)}} M_{23}^{(2)} \\
\end{bmatrix}
\begin{bmatrix}
x_1 \\
x_2 \\
x_3 \\
\end{bmatrix}
=
\begin{bmatrix}
b_1 \\
b_2 – \frac{M_{21}}{M_{11}} b_1 \\
b_3 – \frac{M_{31}}{M_{11}} b_1
– \frac{M_{32}^{(2)}}{M_{22}^{(2)}} b_{2}^{(2)} \\
\end{bmatrix}.

LU decomposition

Through Gaussian elimination, we have transformed the system from

\label{eqn:multiphysicsL4:200}
M x = b

to
\label{eqn:multiphysicsL4:220}
U x = y.

Writing out our Gaussian transformation in the form $$U \overline{{x}} = b$$ we have

\label{eqn:multiphysicsL4:240}
U \overline{{x}} =
\begin{bmatrix}
1 & 0 & 0 \\
-\frac{M_{21}}{M_{11}} & 1 & 0 \\
\frac{M_{32}^{(2)}}{M_{22}^{(2)}}
\frac{M_{21}}{M_{11}}

\frac{M_{31}}{M_{11}}
&

\frac{M_{32}^{(2)}}{M_{22}^{(2)}}
& 1
\end{bmatrix}
\begin{bmatrix}
b_1 \\
b_2 \\
b_3
\end{bmatrix}.

We can verify that the operation matrix $$K^{-1}$$, where $$K^{-1} U = M$$ that takes us to this form is

\label{eqn:multiphysicsL4:260}
\begin{bmatrix}
1 & 0 & 0 \\
\frac{M_{21}}{M_{11}} & 1 & 0 \\
\frac{M_{31}}{M_{11}} &
\frac{M_{32}^{(2)}}{M_{22}^{(2)}} &
1
\end{bmatrix}
\begin{bmatrix}
U_{11} & U_{12} & U_{13} \\
0 & U_{22} & U_{23} \\
0 & 0 & U_{33} \\
\end{bmatrix}
\overline{{x}} = \overline{{b}}

Using this {LU} decomposition is generally superior to standard Gaussian elimination, since we can use this for many different $$\overline{{b}}$$ vectors using the same amount of work.

Our steps are

\label{eqn:multiphysicsL4:280}
b
= M x
= L \lr{ U x}
\equiv L y.

We can now solve $$L y = b$$, using substitution for $$y_1$$, then $$y_2$$, and finally $$y_3$$. This is called {forward substitution}.

Finally, we can now solve

\label{eqn:multiphysicsL4:300}
U x = y,

using {back substitution}.

Note that we produced the vector $$y$$ as a side effect of performing the Gaussian elimination process.

It is already time for remembrance day propaganda

September 24, 2014 Incoherent ramblings , , , ,

The propaganda for celebration of war, destruction, killing, and rape is already starting.  I was greeted today on facebook by the serene and smiling image of a world war II vet, along with his seemingly innocent request for “likes”

This is probably supposed to make people feel patriotism, the worship of an artificial boundary on a map and its associated label system.  Patriotism and national pride worship are designed to ensure that no rational thought opposes the psychopathic figureheads that run our modern “democracies”.

When we see vet worship propaganda like this, you have to think about the interests that encouraged the war in the first place.  These are power brokers like the Carnegies that requested that Wilson not end the war (WWI) too soon (and who had concluded in their minutes before that there was no better institution for profits than the business of war).  These are the armaments and munitions  manufacturers.  These are the companies that built the tanks and the planes that bombed uncounted civilians in Germany and the UK.  These are the makers of weapons that use white phosphorus and depleted uranium, and the makers of cluster bombs that shred children into hamburger meat before and after they are deployed.

I would describe us as being in world war III now.  This is the series of undeclared and covert wars that the USA and NATO have waged since world war II.  This is a war that happens with our tacit approval.  This approval that is assumed because there is no widespread disapproval.  Unfortunately, so much of this modern and current war happens without our even knowing about it, that nobody thinks to object.  This war is disconnected from our immediate attention, and most of the time we don’t even realize that we are funding it.  As always, I don’t want my tax dollars, money stolen from me at the point of a gun, to be used to fund the slaughter of people in Yugoslavia or Libya or <enter your favorite target in some little known corner of the world>, but I can’t control this.  It happens again and again, and is truly disgusting.

I am happy to see that this poor fellow survived.  However, if he had fought a real war, it would have been the war against the propaganda and brainwashing that involved him in the combat that got most of his compatriots killed, and against the same combat that killed uncountable numbers of civilians and “enemies”, and against the combat that encouraged these “enemies” to do the same.

I will not be liking this picture.

ECE1254H Modeling of Multiphysics Systems. Lecture 3: Nodal Analysis. Taught by Prof. Piero Triverio

September 24, 2014 ece1254 , , ,

Disclaimer

Peeter’s lecture notes from class. These may be incoherent and rough.

Nodal Analysis

Avoiding branch currents can reduce the scope of the computational problem. Consider the same circuit fig. 1, this time introducing only node voltages as unknowns

fig. 1. Resistive circuit with current sources

Unknowns: node voltages: $$V_1, V_2, \cdots V_4$$

Equations are KCL at each node except $$0$$.

1. $$\frac{V_1 – 0}{R_A} + \frac{V_1 – V_2}{R_B} + i_{S,A} = 0$$

2. $$\frac{V_2 – 0}{R_E} + \frac{V_2 – V_1}{R_B} + i_{S,B} + i_{S,C} = 0$$

3. $$\frac{V_3 – V_4}{R_C} – i_{S,C} = 0$$

4. $$\frac{V_4 – 0}{R_D} +\frac{V_4 – V_3}{R_C} – i_{S,A} – i_{S,B} = 0$$

In matrix form this is

\label{eqn:multiphysicsL3:20}
\begin{bmatrix}
\inv{R_A} + \inv{R_B} & – \inv{R_B} & 0 & 0 \\
-\inv{R_B} & \inv{R_B} + \inv{R_E} & 0 & 0 \\
0 & 0 & \inv{R_C} & -\inv{R_C} \\
0 & 0 & -\inv{R_C} & \inv{R_C} + \inv{R_D}
\end{bmatrix}
\begin{bmatrix}
V_1 \\
V_2 \\
V_3 \\
V_4 \\
\end{bmatrix}
=
\begin{bmatrix}
-i_{S,A} \\
-i_{S,B} – i_{S,C} \\
i_{S,C} \\
i_{S,A} + i_{S,B}
\end{bmatrix}

Introducing the nodal matrix

\label{eqn:multiphysicsL3:40}
G \overline{{V}}_N = \overline{{I}}_S

We identify the {stamp} for a resister of value $$R$$ between nodes $$n_1$$ and $$n_2$$

stamp matrix

where we have a set of rows and columns for each of the node voltages $$n_1, n_2$$.

Note that some care is required to use this nodal analysis method since we required the invertible relationship $$i = V/R$$. We also cannot handle short circuits $$V = 0$$, or voltage sources $$V = 5$$ (say). We will also have trouble with differential terms like inductors.

Recap of node branch equations

• KCL: $$A \cdot \overline{{I}}_B = \overline{{I}}_S$$
• Constitutive: $$\overline{{I}}_B = \alpha A^\T \overline{{V}}_N$$,
• Nodal equations: $$A \alpha A^\T \overline{{V}}_N = \overline{{I}}_S$$

where $$\overline{{I}}_B$$ was the branch currents, $$A$$ was the incidence matrix, and $$\alpha = \begin{bmatrix}\inv{R_1} & & \\ & \inv{R_2} & \\ & & \ddots \end{bmatrix}$$.

The stamp can be observed in the multiplication of the contribution for a single resistor, where we see that the incidence matrix has the form $$G = A \alpha A^\T$$

stamp factor

Theoretical facts

Noting that $$\lr{ A B }^\T = B^\T A^\T$$, it is clear that the nodal matrix $$G = A \alpha A^\T$$ is symmetric

\label{eqn:multiphysicsL3:80}
G^\T
=
\lr{ A \alpha A^\T }^\T
=
\lr{ A^\T }^\T \alpha^\T A^\T
=
A \alpha A^\T
= G

Modified nodal analysis (MNA)

This is the method that we find in software such as spice.

To illustrate the method, consider the same circuit, augmented with an additional voltage sources as in fig. 4.

fig. 4. Resistive circuit with current and voltage sources

We know wish to have the following unknowns

• node voltages ($$N-1$$): $$V_1, V_2, \cdots V_5$$
• branch currents for selected components ($$K$$): $$i_{S,C}, i_{S,D}$$

We will have two less unknowns for this system than with standard nodal analysis. Our equations are

1. $$-\frac{V_5-V_1}{R_A} +\frac{V_1-V_2}{R_B} + i_{S,A} = 0$$

2. $$\frac{V_2-V_5}{R_E} +\frac{V_2-V_1}{R_B} + i_{S,B} + i_{S,C} = 0$$

3. $$-i_{S,C} + \frac{V_3-V_4}{R_C} = 0$$

4. $$\frac{V_4-0}{R_D} +\frac{V_4-V_3}{R_C} – i_{S,A} – i_{S,B} = 0$$

5. $$\frac{V_5-V_2}{R_E} +\frac{V_5-V_1}{R_A} + i_{S,D} = 0$$

Put into giant matrix form, this is

giant matrix

Call the extension to the nodal matrix $$G$$, the {voltage incidence matrix} $$A_V$$.

ECE1254H Modeling of Multiphysics Systems. Lecture 2: Assembling system equations automatically. Taught by Prof. Piero Triverio

September 23, 2014 ece1254 , ,

Disclaimer

Peeter’s lecture notes from class. These may be incoherent and rough.

Assembling system equations automatically. Node/branch method

Consider the sample circuit of fig. 1.

fig. 1. Sample resistive circuit

Step 1. Choose unknowns:

For this problem, let’s take

• node voltages: $$V_1, V_2, V_3, V_4$$

• branch currents: $$i_A, i_B, i_C, i_D, i_E$$

We do not need to introduce additional labels for the source current sources.
We always introduce a reference node and call that zero.

For a circuit with $$N$$ nodes, and $$B$$ resistors, there will be $$N-1$$ unknown node voltages and $$B$$ unknown branch currents, for a total number of $$N – 1 + B$$ unknowns.

Step 2. Conservation equations:

KCL

• 0: $$i_A + i_E – i_D = 0$$
• 1: $$-i_A + i_B + i_{S,A} = 0$$
• 2: $$-i_B + i_{S,B} – i_E + i_{S,C} = 0$$
• 3: $$i_C – i_{S,C} = 0$$
• 4: $$-i_{S,A} – i_{S,B} + i_D – i_C = 0$$

Grouping unknown currents, this is

• 0: $$i_A + i_E – i_D = 0$$
• 1: $$-i_A + i_B = -i_{S,A}$$
• 2: $$-i_B -i_E = -i_{S,B} -i_{S,C}$$
• 3: $$i_C = i_{S,C}$$
• 4: $$i_D – i_C = i_{S,A} + i_{S,B}$$

Note that one of these equations is redundant (sum 1-4). In a circuit with $$N$$ nodes, we can write at most $$N-1$$ independent KCLs.

In matrix form

\label{eqn:multiphysicsL2:20}
\begin{bmatrix}
-1 & 1 & 0 & 0 & 0 \\
0 & -1 & 0 & 0 & -1 \\
0 & 0 & 1 & 0 & 0 \\
0 & 0 & -1 & 1 & 0
\end{bmatrix}
\begin{bmatrix}
i_A \\
i_B \\
i_C \\
i_D \\
i_E \\
\end{bmatrix}
=
\begin{bmatrix}
-i_{S,A} \\
-i_{S,B} -i_{S,C} \\
i_{S,C} \\
i_{S,A} + i_{S,B}
\end{bmatrix}

We call this first matrix of ones and minus ones the incidence matrix $$A$$. This matrix has $$B$$ columns and $$N-1$$ rows. We call the known current matrix $$\overline{I}_S$$, and the branch currents $$\overline{I}_B$$. That is

\label{eqn:multiphysicsL2:40}
A \overline{I}_B = \overline{I}_S.

Observe that we have both a plus and minus one in all columns except for those columns impacted by our neglect of the reference node current conservation equation.

Algorithm for filling $$A$$

In the input file, to describe a resistor of fig. 2, you’ll have a line of the form

R name $$n_1$$ $$n_2$$ value

fig. 2. Resistor node convention

The algorithm to process resistor lines is

\begin{equation*}
\begin{aligned}
&A \leftarrow 0 \\
&ic \leftarrow 0 \\
&\text{forall resistor lines} \\
&ic \leftarrow ic + 1, \text{adding one column to $$A$$} \\
&\text{if}\ n_1 != 0 \\
&A(n_1, ic) \leftarrow +1 \\
&\text{endif} \\
&\text{if}\ n_2 != 0 \\
&A(n_2, ic) \leftarrow -1 \\
&\text{endif} \\
&\text{endfor} \\
\end{aligned}
\end{equation*}

Algorithm for filling $$\overline{I}_S$$

Current sources, as in fig. 3, a line will have the specification (FIXME: confirm… didn’t write this down).

I name $$n_1$$ $$n_2$$ value

fig. 3. Current source conventions

\begin{equation*}
\begin{aligned}
&\overline{I}_S = \text{zeros}(N-1, 1) \\
&\text{forall current lines} \\
&\overline{I}_S(n_1) \leftarrow \overline{I}_S(n_1) – \text{value} \\
&\overline{I}_S(n_2) \leftarrow \overline{I}_S(n_1) + \text{value} \\
&\text{endfor}
\end{aligned}
\end{equation*}

Step 3. Constitutive equations:

\label{eqn:multiphysicsL2:60}
\begin{bmatrix}
i_A \\
i_B \\
i_C \\
i_D \\
i_E
\end{bmatrix}
=
\begin{bmatrix}
1/R_A & & & & \\
& 1/R_B & & & & \\
& & 1/R_C & & & \\
& & & 1/R_D & & \\
& & & & & 1/R_E
\end{bmatrix}
\begin{bmatrix}
v_A \\
v_B \\
v_C \\
v_D \\
v_E
\end{bmatrix}

Or

\label{eqn:multiphysicsL2:80}
\overline{I}_B = \alpha \overline{V}_B,

where $$\overline{V}_B$$ are the branch voltages, not unknowns of interest directly. We can write

\label{eqn:multiphysicsL2:100}
\begin{bmatrix}
v_A \\
v_B \\
v_C \\
v_D \\
v_E
\end{bmatrix}
=
\begin{bmatrix}
-1 & & & \\
1 & -1 & & \\
& & 1 & -1 \\
& & & 1 \\
& -1 & &
\end{bmatrix}
\begin{bmatrix}
v_1 \\
v_2 \\
v_3 \\
v_4
\end{bmatrix}

Observe that this is the transpose of $$A$$, allowing us to write

\label{eqn:multiphysicsL2:120}
\overline{V}_B = A^\T \overline{V}_N.

Solving

• KCLs: $$A \overline{I}_B = \overline{I}_S$$

• constitutive: $$\overline{I}_B = \alpha \overline{V}_B \Rightarrow \overline{I}_B = \alpha A^\T \overline{V}_N$$

• branch and node voltages: $$\overline{V}_B = A^\T \overline{V}_N$$

In block matrix form, this is

\label{eqn:multiphysicsL2:140}
\begin{bmatrix}
A & 0 \\
I & -\alpha A^\T
\end{bmatrix}
\begin{bmatrix}
\overline{I}_B \\
\overline{V}_N
\end{bmatrix}
=
\begin{bmatrix}
\overline{I}_S \\
0
\end{bmatrix}.

Is it square? We see that it is after observing that we have

• $$N-1$$ rows in $$A$$ , and $$B$$ columns.
• $$B$$ rows in $$I$$.
• $$N-1$$ columns.

ECE1254H Modeling of Multiphysics Systems. Lecture 1: Analogies to circuit systems. Taught by Prof. Piero Triverio

September 22, 2014 ece1254

Disclaimer

Peeter’s lecture notes from class. May not be entirely coherent.

In slides

A review of systematic nodal analysis for a basic resistive circuit was outlined in slides, with a subsequent attempt to show how many similar linear systems can be modeled as circuits so that the same toolbox can be applied. This included blood flow through a body (and blood flow to the brain), a model of antenna interference in a portable phone, heat conduction in a one dimensional conductor under a heat lamp, and a few other systems.

This discussion reminded me of the joke where the farmer, the butcher and the physicist are all invited to talk at a beef convention. After meaningful and appropriate talks by the farmer and the butcher, the physicist gets his chance, and proceeds with “We begin by modeling the cow as a sphere, …”. The ECE equivalent of that appears to be a Kirchhoff circuit problem.

Mechanical structures example

Continuing the application of circuits like linear systems to other systems, let’s consider a truss system as illustrated in fig. 1, or in the simpler similar system of fig. 2.

fig. 1. A static loaded truss configuratio

Our unknowns are

• positions of the joints after deformation $$(x_i, y_i)$$.
• force acting on each strut $$\BF_j = (F_{j,x}, F_{j,y})$$.

The constitutive equations, assuming static conditions (steady state, no transients)

• Load force. $$\BF_L = (F_{L, x}, F_{L, y}) = (0, -m g)$$.
• Strut forces. Under static conditions the total resulting force on the strut is zero, so $$\BF’_j = -\BF_j$$. For this problem it is redundant to label forces on both ends, so we mark the labeled end of the object with an asterisk as in fig. 3.

fig 3. Strut model

Consider a simple case

One strut as in fig. 4.

fig. 4. Very simple static load

\label{eqn:multiphysicsL1:20}
\BF^\conj = – \Ba_x
\underbrace{
\epsilon
}_{\text{constant, describes the beam elasticity, given}}
\biglr{
\underbrace{
L
}_{\text{unloaded length $$L = \Abs{x^\conj – 0}$$, given}}
– L_0}

The constitutive law for a general strut as in fig. 5 is

fig 5. Strut force diagram

The force is directed along the unit vector

\label{eqn:multiphysicsL1:40}
\Be = \frac{\Br^\conj – \Br}{\Abs{\Br^\conj – \Br}},

and has the form
\label{eqn:multiphysicsL1:60}
\BF^\conj = – \Be \epsilon \lr{ L – L_0 }.

The value $$\epsilon$$ may be related to Hooks’ constant, and $$L_0$$ is given by

\label{eqn:multiphysicsL1:80}
L = \Abs{\Br^\conj – \Br} = \sqrt{(x^\conj – x)^2 + (y^\conj – y)^2}.

Observe that the relation between $$\BF^\conj$$ and position is nonlinear!

Treatment of this system will be used as the prototype for our handling of other nonlinear systems.

Returning to the simple static system, and introducing force and joint labels as in fig. 6, we can examine the \textAndIndex{conservation law}, the balance of forces.

fig 6. Strut system

• At joint 1:\label{eqn:multiphysicsL1:100}
\Bf_A + \Bf_B + \Bf_C = 0
or
\label{eqn:multiphysicsL1:120}
\begin{aligned}
\Bf_{A,x} + \Bf_{B,x} + \Bf_{C,x} &= 0 \\
\Bf_{A,y} + \Bf_{B,y} + \Bf_{C,y} &= 0
\end{aligned}
• At joint 2:\label{eqn:multiphysicsL1:140}
-\Bf_C + \Bf_D + \Bf_L = 0
or
\label{eqn:multiphysicsL1:160}
\begin{aligned}
-\Bf_{C,x} + \Bf_{D,x} + \Bf_{L,x} &= 0 \\
-\Bf_{C,y} + \Bf_{D,y} + \Bf_{L,y} &= 0
\end{aligned}

We have an equivalence

• Force $$\leftrightarrow$$ Current.
• Force balance equation $$\leftrightarrow$$ KCL