
\documentclass[12pt,ragged2e,oneside,a4paper]{book}

% For embeddeing the source file
\usepackage{attachfile}
\usepackage{caption}

\renewcommand{\captionfont}{\small}


\usepackage[english]{babel}
\usepackage{tocloft}
\usepackage{titlesec}
\usepackage{enumitem}
\usepackage{marginnote}
\usepackage{graphicx}
% for different enumerate styles
\usepackage{multicol}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{amsmath}
\usepackage{amsthm}
\usepackage{bookmark}
\usepackage{pgfplots}
\usepackage{import}
\usepackage{xcolor}
\usepackage{soul}
\newcommand{\mathcolorbox}[2]{\colorbox{#1}{$\displaystyle #2$}}

\usepgfplotslibrary{external}
\tikzexternalize
\pgfplotsset{width=8cm,compat=1.9}

\usepackage[pdf]{graphviz}
% for unnumbered theorems

\usepackage{mathtools}
% for :=

\usepackage{amssymb}
% for mathbb 



\usepackage[top=2cm, bottom=3cm, outer=7cm, inner=1cm, heightrounded,
marginparwidth=5cm, marginparsep=0.3cm]{geometry}

\title{MATH3102 Lecture Notes}
\author{Alistair Michael}

\titleformat{\chapter}{\normalfont\huge\bf}{\thechapter.}{20pt}{\huge\bf}

\graphicspath{{images/}}

\newcommand{\sforall}{\enspace\forall}

\renewcommand{\labelenumi}{\roman{enumi})}
\newcommand{\dsum}{\displaystyle\sum}
\newcommand{\ds}{\displaystyle}
\newcommand{\lsum}{\underline{\int_a^b}}
\newcommand{\usum}{\overline{\int_a^b}}
\newcommand{\rieman}{\int_a^b}
\newcommand{\Arg}{\text{Arg}\,}
\newcommand{\Log}{\text{Log}\,}
\newcommand{\Int}{\text{Int}\,}
\newcommand{\Ext}{\text{Ext}\,}
\newcommand{\p}{\partial}
\newcommand{\pd}{\partial}

\newtheoremstyle{theoremstyle} % name
    {\topsep}                    % Space above
    {\topsep}                    % Space below
    {}                   % Body font
    {}                           % Indent amount
    {\itshape}                   % Theorem head font
    {}                          % Punctuation after theorem head
    {.5em}                       % Space after theorem head
    {}  % Theorem head spec (can be left empty, meaning ‘normal’)

\theoremstyle{theoremstyle}
\newtheorem{thm}{Theorem}[subsection]
\newtheorem{lemma}[thm]{Lemma}
\newtheorem{corr}{Corrolary}[thm]
\newtheorem{claim}{Claim}[thm]

\newtheorem{definition}{Definition}[subsection]
\newtheorem{prop}{Proposition}[subsection]
\newtheorem{ex}{Example}[subsection]

\newtheoremstyle{ittheoremstyle} % name
    {\topsep}                    % Spaceabove
    {\topsep}                    % Space below
    {\itshape}                   % Body font
    {}                           % Indent amount
    {\scshape}                   % Theorem head font
    {.}                          % Punctuation after theorem head
    {.5em}                       % Space after theorem head
    {}  % Theorem head spec (can be left empty, meaning ‘normal’)

\theoremstyle{ittheoremstyle}


\newcommand{\incfig}[1]{%
    \def\svgwidth{\columnwidth}
    \import{./img/}{#1.pdf_tex}
}

\newcommand{\incfigw}[2]{%
    \begin{center}
    \def\svgwidth{#2\columnwidth}
    \import{./img/}{#1.pdf_tex}
    \end{center}
}


\newcommand{\N}{\mathbb N}
\newcommand{\R}{\mathbb R}
\newcommand{\Z}{\mathbb Z}
\newcommand{\C}{\mathbb C}
\newcommand{\Q}{\mathbb Q}
\newcommand{\contra}{\Rightarrow\Leftarrow}

\setlength{\parskip}{0.5em}

\setlength{\parindent}{0em}

\newcommand{\nombreindice}{Contents By Week}
\newlistof{lecture}{lecn}{\nombreindice}

\newcounter{weekno}
\newcounter{lecno}

% definition of the commands used for the Spanish ToC;
% \captce for chapters, \sectce for sections and
% \ssectce for subsections
\newcommand\captce[1]{%
  \addcontentsline{lecn}{chapter}{\protect\makebox[1.3em][l]{\thechapter}#1}}
  \newcommand\week[1]{\stepcounter{weekno}\setcounter{lecno}{0}
  \addcontentsline{lecn}{section}{\protect\makebox[4.8em][l]{Week \theweekno}#1}}
  \newcommand\lecture[1]{\stepcounter{lecno} 
\addcontentsline{lecn}{subsection}{\protect\makebox[6em][l]{Lecture \thelecno}#1}
\marginnote{(Lecture \theweekno.\thelecno)}
}

\renewcommand{\cfttoctitlefont}{\normalfont\bfseries\Large}
\renewcommand{\cftlecntitlefont}{\normalfont\bfseries\Large}

\begin{document}



\maketitle

This document can be found at
\href{https://alistairmichael.com/notes/applied.pdf}{alistairmichael.com/notes/applied.pdf} with source available at
\href{https://git.topost.net/alistair/applied-notes}{git.topost.net/alistair/applied-notes}.
It contains many typos and errors.

\tableofcontents
\listoflecture

\pagebreak


\week{}
\chapter{Introduction}
\lecture{}

Applied mathematics is about creating and analyzing mathematical models of 
phenomena and processes in the real world. 

The textbook is M. Holmes \emph{Introduction to Foundations of Applied
Mathematics, Springer.}


\includegraphics[height=0.5\linewidth]{ab.pdf}


\textbf{Modelling:} Creating models from words to equations, identifying key
variables and processes. Models are as simple as possible.

\emph{Every model is wrong, ... some are useful.}

Limited insight and predictive power in highly complex systems or systems that
are very sensitive to the initial conditions. 

\section{Topics}

\begin{itemize}
    \item Dimensional analysis
    \item Perturbation methods
    \item Traffic flow models
    \item Continuum mechanics in 1D
    \item Elastic and viscoelastic materials
    \item Continuum mechanics in 3D
    \item Fluid flows
\end{itemize}


\chapter{Dimensional Analysis}

\textbf{Dimensional Estimate:} Qualitative prediction of the order of magnitude
of key quantities without using a model.

\textbf{Non-dimensionalisation:} Reducing the numbner of effective indpendent
parameters.

So to create a dimensional estimate for the period of a pendulum $P$, we
understand it may depend on time, distance ($L$), velocity, mass and gravity.

$P = f(L, v_0, M, g)$

How can we combine these paramaters to get the correct dimension as a result?

$$\frac L  g = \frac {m} {m / s^2}$$

So as an estimate we can say $P = const \times \sqrt{ \frac L g }$. 

\lecture{}

\begin{ex}
Ex. Energy of the Bomb

We want to determine the energy released by a bomb using dimensional analysis. 

We know the radius of the shockwave increases with time, at a rate which depends on the
air density and the energy released. 

\textbf{Radius $R$:} $d$

\textbf{Energy:} $m \frac {d^2}{t^2}$

\textbf{Air Density ($\rho$):} $\frac {m}{d^3}$

So $\frac \rho E = \frac {t^2} {d^5} \implies t^2 \frac E \rho = d^5 $

So $R(t) = C \cdot \frac E \rho ^{\frac 15} \cdot t^{\frac 25}$ 
\end{ex}

\subsubsection{Non-dimensionalization}

\begin{enumerate}
    \item[1.] Introduce new units to each independent quantity to reduce the number
        of parameters.
    \item[2.] Identify relative magnitude of processes and identify dominant
        terms and small perturbations. 
    \item[3.] Choose the non-dimensional form of the equation so that all terms
        are either order 1 or small (perturbation).
\end{enumerate}

\begin{ex} Pendulum

    $\phi$ is angle to overtical, $l$ is length of the pendulum.



    \begin{align*}
        ml \frac {d^2 \phi}{dt^2} &= -mg \sin(\phi) \\
         \frac {d^2 \phi}{dt^2} &= -\frac gl \sin(\phi) \\
    \end{align*}

    $\frac gl$ can be eliminited by choosing a new time unit.

    Call $t' = \frac t {t_0}$ the nondimensional time unit. 

    $g$ = distance / time$^2$. $l$ = distance, $\phi = $  1. $\sqrt{{\frac lg}}
    = $ time $\implies $ choose $t_0 = \sqrt{\frac lg} $ 

    \marginnote{$\phi$ is a ratio of two distances it doesnt have a dimension.}
    $$\frac {d^2 \phi}{d{t'}^2} = -\sin (\phi)$$
    
    Which has no parameters.
\end{ex}

\begin{ex} Damped Pendulum

    \begin{align*}
        ml \frac {d^2\phi}{dt^2} &= - mg\sin(\phi) - kl\frac {d\phi}{dt} \\
        \frac {d^2\phi}{dt^2} &= - \frac gl\sin(\phi) - \frac km\frac {d\phi}{dt} \\
    \end{align*}

    So there are two parameters $\frac km$ and $\frac gl$.

    So choose a new time unit $t = t' t_0$  to make one of the coefficients
    disappear. And make $t_0 = \sqrt{\frac lg} $

    \begin{align*}
    \frac {d^2\phi}{dt^2} &= - \frac gl\sin(\phi) {t_0}^2 - \frac km\frac
        {d\phi}{dt} t_0 \\
        \frac {d^2\phi}{dt^2} &= - \sin(\phi)- \frac km\frac
        {d\phi}{dt} \sqrt{\frac lg}  \\
        \frac {d^2\phi}{dt^2} &= - \sin(\phi)- \alpha\frac
        {d\phi}{dt} 
    \end{align*}

    Natural time unit proportional to the period. 
    
    The solution depends on a single parameter $\alpha = (\frac km)\sqrt{\frac l g}$ 
\end{ex}

\begin{ex}
    Rabbits and Foxes.

    \begin{align*}
        \frac {dR}{dt} &= rR \left(1- \frac RK\right) - pF \frac R{R + S} \\
        \frac {dF}{dt} &= \gamma p F \frac {R}{R + S} - \delta F \\
    \end{align*}

    We want to reduce the parameters: it is a six-dimensional parameter space
    for a model of two species.

    \marginnote{So the name of the game is choosing units to absorb all the coefficients.}
    Choose natural units for the variables $R$, $F$ and $t$: which look like the
    most concrete ones. 

    \begin{align*}
        R &= R_0 R' \\
        F &= F_0 F' \\
        t &= t_0 t' \\
    \end{align*}

    Exercise.
\end{ex}

\lecture{}

\subsubsection{Non-Dimensional Variables and Parameters}

Where we have a variable $x$ and introduce $x = x_0 x'$ such that $x_0$
minimises the number of parameters in the model, we call $x'$ a non-dimensional
variable.

\begin{ex} Cancer Growth (PDE)

    $$\frac {\partial C}{\partial t} = rC \left(1 - \frac CK\right) + D \Delta C$$

    \begin{itemize}
        \item First term is growth from cell division.
        \item Second term is diffusion due to random cell movements.
    \item $\Delta$ is the laplacian. $\sum \frac {\partial^2}{\partial {x_i}^2}$
    \end{itemize}


    What are the dimensions of $r$ and $D$?

    \begin{itemize}
        \item $r = \frac 1 {time}$ 
        \item $D = \frac {dist^2}{time}$
    \end{itemize}

    We want to get something with the dimension $\frac {dist}{time}$.

    So $v \sim \sqrt{{rD}}$.  

    So use $t_0 = \frac 1r$ $x_0 = \sqrt{{\frac Dr}}$ and $C_0 = K$. 

    \begin{align*}
        C &= C'K \\
        t &= \frac {t'}r \\
        x &= x'\sqrt{\frac Dr} \\
        \\
        \frac {\partial C'}{\partial t'} &= C' (1 - C') + \Delta' C'
    \end{align*}

\end{ex}

\chapter{Perturbation Methods}

    Transform the original model into non-dimensional form and identify a small
    non-dimensional parameter $\varepsilon < 1$.

    Assume that in the special case $\varepsilon = 0$, that the model simplifies
    into a problem that is fully solvable analytically. 

    This method gives good approximations for small epsilons. 

    \begin{ex}
        $m \frac {d^2 y}{dt^2} = -mg $, $y(0) = 0$, $\frac {dy}{dt}|_{t = 0} = v_0$.

        Introduce a nondimensional variables for $t$ and height $y$. 

        $y = y_0 t'$ and $t = t_0 t'$ so choose $t_0 = v_0 / g$ and $y_0 =
        {v_0}^2 / g$ which leads to the non-dimensional problem

        $$\frac {d^2 y'}{dt'^2} = -1,\quad\quad y'(0)=0,\quad\quad \frac
        {dy'}{dt'}|_{t' = 0} = 1$$

        Its still fully solvable: $y'(t') = - \frac 12 {t'}^2 + t'$. If we want
        to generalise it, for example saying gravitational acceleration is no
        longer constant. Where $R$ is the radius of the Earth then

        $$\frac {d^2y}{dt^2} = - g \frac {R^2}{(R + y)^2}$$

        Using the same non-dimensionalisation as before we get:

        $$\frac {d^2 y'}{dt'^2} = -\left(1 + \frac{{v_0}^2}{Rg}y'\right)^{-2},\quad y'(0)=0,\quad \frac
        {dy'}{dt'}|_{t' = 0} = 1$$

        So the small non-dimensional parameter is $\varepsilon = \frac {{v_0}^2}{Rg}$

        So can find approximate solutions when $\varepsilon \ll 1$ is close to $0$.


    \end{ex}



    \section{Perturbation Methods}

    \week{}
    \lecture{}

    \begin{ex}
        Solve $x^2 + 2\varepsilon x - 1 = 0$ with perturbation methods. 

        The equation has a small non-dimensional parameter $\varepsilon$.

        It can be solved exactly: $x_{1,2} = -\varepsilon \pm
        \sqrt{{\varepsilon^2} + 1}$ but we will use it as an example. 
        See it is easily solvable when $\varepsilon = 0$.

        When $\varepsilon \ll 1$ we can use a taylor series expansion (which is
        only valid around $0$ when $\varepsilon \ne 0$).

        {
        $$x_{1,2}(\varepsilon) \approx \pm 1 - \varepsilon \pm \frac 12 \varepsilon^2 \mp \frac 18 \varepsilon^4 + \cdots$$
        \captionof{figure}{Taylor Series Expansion}}

        Look for an approximation of the form: $x(\varepsilon) = c_0 + x_1
        \varepsilon + x_2 \varepsilon^2$

        Assume $\varepsilon \ll 1$ and $x_0, x_1, x_2,\dots$ are unknown
        coefficients of magnitude $\sim 1$. We want to be able to order the
        terms by the magnitude according only to the power of $\varepsilon$.

        In $(x_0 + x_1\varepsilon + \dots)^2 + 2\varepsilon(x_0 + x_1\varepsilon
        + \dots) - 1 = 0$ it is clear that the $x_0$ terms dominate, and we can
        take out the dominant terms and group them by the powers of epsilon. 

        The largest terms where $\varepsilon ^ 0$

        $${x_0}^2 - 1 = 0 \quad \quad x_0 = \pm 1$$

        The next largest terms, terms proportional to $\varepsilon^1$: 

        $$2x_0 x_1 \varepsilon + 2\varepsilon x_0 = 0 \quad \quad x_1 = -1$$

        Terms proportional to $\varepsilon^2$: 

        $${x_1}^2 \varepsilon^2 + 2 x_0 x_2 \varepsilon^2 + 2x_1 \varepsilon^2 =
        0 \quad \quad x_2 = \pm \frac 12$$

        We can neglect the rest assuming $\varepsilon^3$ is insignificantly
        small.

        So we have the approximate solution: 

        \begin{align*}
            x(\varepsilon) &= x_0 +  x_1 \varepsilon + x_2 \varepsilon^2 + \cdots \\
            x(\varepsilon) &= \pm 1 - \varepsilon \pm \frac 12 \varepsilon^2 + \cdots 
        \end{align*}
        \marginnote{Which matches the taylor series.}


        All we need to be able to solve these kinds of problems is to be able to
        solve for the coefficients along the way and use the approximation of
        $\varepsilon$ close to $0$: we do not need the explicit form that we
        have in this example. 


    \end{ex}

    \textbf{Perturbation method: Solve $F(y; \varepsilon) = 0$ for $\varepsilon \ll 1$ assuming
    that $F(y; 0) = 0$ has an exact solution $y_0$ to find an approximate
    analytical solution.}

    \begin{definition}{Perturbation Method}{}

        Look for a solution in the form of an asymptotic expansion

        $$y(\varepsilon) = y_0 + y_1 f_1(\varepsilon) + y_2f_2(\varepsilon) +
        \cdots $$

        Where $1 \gg f_1(\varepsilon) \gg f_2(\varepsilon) \gg \cdots$ when
        $\varepsilon \ll 1$ (for example a power series like in Example 2.0.1).
        
    \end{definition}
    
    \begin{ex}
        Find an asymptotic solution for $x^3 + \varepsilon x - 1 = 0$ when
        $\varepsilon \ll 1$

        \marginnote{These problems are artificially simple but it is fine for the
        coefficients to contain parameters.}

        Special case: $\varepsilon = 0,\; x^3 - 1 = 0$ has the exact solution $x =
        1$.

        $(x_0 + x_1 \varepsilon + x_2\varepsilon^2 \cdots) ^3 + \varepsilon(x_0
        + x_1\varepsilon + x_2 \varepsilon^2 \cdots) - 1 = 0$

        \begin{align*}
            e^0 &: x_0^3 - 1 = 0 &\implies x_0 = 1 \\
            e^1 &: 3x_0^2x_1 + x_0 = 0, \; 3x_1 + 1 = 0 &\implies x_1 = \frac 13 \\
            e^2 &: 3x_0^2x_2 +  3x_1^2x_0 + x_1 = 0,\; 3x_2 + \frac 13 - \frac 13 =
            0 &\implies x_2 = 0\\
            e^3 &: {...} &\implies x_3 = \frac 1 {81}
        \end{align*}


        {
        $$x = 1 - \frac 13 \varepsilon + \frac 1 {81}\varepsilon^3 + \cdots$$
        \captionof{figure}{The Asymptotic Series Solution}
        }

    \end{ex}

    \begin{ex}
        Point mass in a non-constant gravitational field:

        \begin{equation}
        \frac {d^2y'}{dt'^2} = - \left(1 + \frac {{v_0}^2}{Rg}y'\right)^{-2}, y'(0)
            = 0, \frac {dy'}{dt'}(t = 0) = 1
        \end{equation}

        Small non-dimensional parameter $\varepsilon = \frac {{v_0}^2}{Rg}$. It
        is not solvable for nonzero $\varepsilon$ but can be approximated when
        $\varepsilon \ll 1$.

        \marginnote{You could drop the primes after the change of variables to
        make the notation less confusing.}

        $$\varepsilon = 0 \implies \frac {d^2y'}{dt'^2} = -1 $$

        $y'(t') = - \frac 1 2 t'^2 + t'$

        $$\frac {d^2y'}{dt'^2}\left(1 + \frac {{v_0}^2}{Rg}y'\right)^{2} = -1 $$

        \marginnote{Rearranging (1)}

        \begin{align*}
            \varepsilon^0:& \frac {d^2y_0}{dt'^2} = -1, y_0(0) = 0, \frac
            {dy_0}{dt'}(t'=0) = 1 && y_0(t') = - \frac 12 t'^2 + t'  \\
            \varepsilon^1:& \frac {d^2y_1}{dt'^2} + 2y_0 \frac {d^2y_0}{dt'^2} =
            0 &&\implies \frac{d^2 y_1}{2t'^2} = -t'^2 + 2t' \\ 
              &\text{And solve with initial condition} 
              && \implies y_1(t') = - \frac 13 t'^3 + \frac 1 {12}t'^4 
        \end{align*}

        $$y'(t') = - \frac 12 t'^2 + t' - \varepsilon \frac {t'^3(4 - t')}{12} +
        \cdots$$

    \end{ex}

    \subsection{Singular Perturbation}

    \lecture{}

    There are some cases where $\varepsilon = 0$ and small $\varepsilon \ne 0$
    are actually qualitatively different (according to the model). 

    A typical case is when setting $\varepsilon = 0$ results in a change of
    order of the equation (i.e. when the epsilon is on the highest order term).  
    These are called \emph{singular perturbation} problems, and the standard
    method does not work. 

    \begin{ex}  $\varepsilon ^2 + x - 1 = 0, \; \varepsilon << 1$

        If we use the regular perturbation method: 

        \[
        x(\varepsilon) = x_0 + x_1 \varepsilon + x_2 \varepsilon^2 + \dots
        \quad\quad\qquad (*)
         \]

        \marginnote{Substitute $x(\varepsilon)$ into the equation.}
        $\varepsilon(x_0 + x_1 \varepsilon + \dots) ^2 + (x_0 + x_1\varepsilon +
        \dots) - 1 = 0$


        Largest $\varepsilon^0$: $x_0 - 1 = 0 \implies x_0 = 1$.

        \marginnote{Using solution $x_0 = 1$}
        $\varepsilon^1: $  $x_0^2 + x_1 = 0 \implies x_1 = -1$
        
        So $x = 1 - \varepsilon + 2\varepsilon^2 + \dots$

        What is the problem with this approach? This is a quadratic equation but
        we only got \textbf{one solution.} This is exactly solvable, so we know
        the solution should be $x_{1,2} = \frac 1 {2\varepsilon}(-1 \pm
        \sqrt{{1-4\varepsilon}})$. 

        \begin{center}
            \begin{tikzpicture}{testp}
        \begin{axis}[xlabel = $\varepsilon$, ylabel=$x$, xtick={0,0.1}]
                \addplot[domain=0:0.1, samples=90,color=blue,]{ 1/(2 * x) * (-1 + sqrt(1 -4 * x))};
                \addplot[domain=0:0.1, samples=90,color=blue]{ 1/(2 *x) * (-1 - sqrt(1 - 4* x))};
            \end{axis}
        \end{tikzpicture}
        \end{center}

        Graphing this solution we can see we have large negative numbers in $x$ as
        $\varepsilon$ approaches zero so the assumption about the form of the
        solution $(*)$ is not correct: all out coefficients are not order $1$.
        
        When we set $\varepsilon = 0$ it reduces the order of the equation to
        $1$. 
    \end{ex}

        \subsubsection{Solving Singular Perturbation Problems}

        Use a transformation of variables $x = e^\alpha z$, where
        $\alpha$ is determined based on the condition that the highest order
        term participates in the dominant part of the equation when $\varepsilon
        \to 0$. Then solve it as a regular perturbation problem.

        \begin{ex} Solving $\varepsilon x^2 + x - 1 = 0$.

            \begin{align*}
                x = \varepsilon^\alpha z &\implies \varepsilon^{2\alpha + 1}z^2 +
            \varepsilon^\alpha z - 1 = 0 
            \end{align*}

            We need to match terms of the same order, and there always needs to
            be at least two of them so they can cancel out to zero. So we need
            to decide which powers of epsilon match up. 

            Either: 

            \begin{itemize}
                \item $2 \alpha + 1 = 0$
                \item $\alpha = 0$, or
                \item $2 \alpha + 1 = \alpha$
            \end{itemize}


            Case 1: 
            \begin{align*}
                2 \alpha + 1 &= 0 \implies \alpha = \frac 12 \\
                \varepsilon^0 z^2 + \varepsilon ^{-\frac 12} - 1 &= 0\\
                \varepsilon^0 z^2 - 1 &= 0 
            \end{align*}

            When $\varepsilon \to 0$ the first and last terms are of the same
            magnitude but the middle is larger: even though we chose the first
            and last to be the dominant terms. This doesnt work so we have to
            drop this case. 

            Case 2: $\alpha = 0$

            \begin{align*}
                \varepsilon z^2 + \varepsilon^0z - 1 &= 0 \\
                \varepsilon z^2 + z - 1 &= 0 \\
            \end{align*}

            This is not even a transformation and it is stilla singular
            perturbation problem. 

            Case 3: $2\alpha + 1 = \alpha \implies \alpha = -1$

            \begin{align*}
                \varepsilon^{-1} z^2 + \varepsilon^{-1} - 1 &= 0 \\
                \implies z^2 + z - \varepsilon &= 0 \\
            \end{align*}

            So now the third term $-1$ is small compared to the very large
            $\varepsilon^{-1}$ when $\varepsilon \to 0$. If we set $\varepsilon
            = 0$ we still have two roots. 
            \begin{itemize}
                \item The highset order term is part of the dominant balance.
                \item We don't loose one of the solutions when $\varepsilon =
                    0$.
            \end{itemize}


            Substitute $z = z_0 + z_1 \varepsilon + \dots$ to get
            \begin{align*}
                (z_0 + z_1\varepsilon + z_2\varepsilon^2+\dots)^2 + (z_0 + z_1
                \varepsilon + z_2 \varepsilon^2) - \varepsilon &= 0 \\
                z_0^2 + z_0 = 0 \implies z_0 &= 0, -1 \\
                2z_0 z_1 + z_1 - 1 = 0 \implies z_1 &= 1,-1 \\
                \dots
            \end{align*}

            Two solutions: $z = \varepsilon - \varepsilon^2 + \dots$ and $z = -1
            - \varepsilon + \varepsilon^2 + \dots$.

            Now do the backwards transformation $x = e^{-1}z$: 
            \begin{align*}
                x &= 1 - \varepsilon + \cdots \\
                x &= -\varepsilon^{-1} - 1 + \varepsilon + \cdots
            \end{align*}
        \end{ex}

        \hrule
        \lecture{}

        Generally considering an $f(x;\varepsilon)$, $x = ?$ when $\varepsilon
        \ll 1$ epsilon allowing us to approximate a solution.  
        \begin{itemize}
            \item Substitute a series approximation: $x = x_0 + x_1 \varepsilon
                + x_2 \varepsilon^2 + \cdots$
            \item Use larger numbers of terms and eventually truncate when there
                is enough accuracy. 
        \end{itemize}

        Boundary conditions:

        Eg $\varepsilon y'' 2y' + 2y = 0$ want to find $y$ in $[0,1]$ with
        boundary conditions for example $y(0) = 0$ and $y(1) = 1$.

        Note $\varepsilon$ multiplies the term with the highset derivative. If
        $\varepsilon = 0$ then it becomes a first order ode and cannot satisfy
        both boundary conditions, so it is a singular perturbation problem.

        In this case we can directly solve it for reference: 

        \begin{center}
            \begin{tikzpicture}{testp}
        \begin{axis}[width=0.9\linewidth,height=2in,xlabel = $x$, ylabel=$y$, xtick={0,0.5}]
                \addplot[domain=0:0.5, samples=200,color=blue]{-2.71964267074 *
                    (exp(-18.994 * x) -  exp(-1.0553 * x ) };
            \addlegendentry{$\varepsilon= 0.1$}
                \addplot[domain=0:0.5, samples=200,color=red]{-2.73205 *
                    (exp(-198.995 * x) -  exp(-1.00505 * x ) };
            \addlegendentry{$\varepsilon= 0.1$}
                \addplot[domain=0:0.5, samples=200,color=green]{-2.71964267074 *
                    (exp(-1999 * x) -  exp(-1.0005 * x ) };
            \addlegendentry{$\varepsilon= 0.1$}
            \end{axis}
        \end{tikzpicture}
        \end{center}

        You can see there is a sharp transition point in the solution
        called `boundary layers'

        \subsubsection{Boundary Layers, Matched Asymptotic Expansions}
        
        For a singular perterbation in ODE boundary value problems, for example
        $\varepsilon y '' + 2y' + 2y = 0 $ where $y(0) = 0, y(1) = 1$, obviously
        this is singular since $\varepsilon = 0$ makes the highest order term
        disappear. Standard perterbation method is only good far from the sharp
        transition point shown above---the `boundary layers'.

        Matched asymptotic expansions makes an approximation that is valid in
        the entire domain. 

        \begin{itemize}
            \item Use rescaling $x = \varepsilon^{\alpha} z$. We want to create
                a transformation that absorbs the singularity so we can solve
                the transformed equation using the normal perterbation method
                and then transform it back.  


                So if we were to set $\varepsilon = 0$ the order should not
                change. 
            \item Chain rule transforms derivatives:

                $Y'(z) = \frac{dY}{dx} \frac {dx}{dz} = y'(x) \varepsilon^a$

                $Y''(z) = y''(x) \varepsilon^{2a} $
            \item After change of variables to the rescaled equation do a
                regular perturbation so we still get a nonzero second derivative
                term. 

                 $\varepsilon y''(x) 2y'(x) + 2y(x) = 0 \to$

               Rescaled equation: $\varepsilon^{1-2\alpha}Y''(z) + \varepsilon^{-\alpha}2Y'(z) +
                2Y(z) = 0$

                \week{}\lecture{}
            \item Determine $\alpha$

                Choose an alpha so the term with the highest derivative is part
                of the dominant balance. 
                \begin{enumerate}
                    \item $1 - 2 \alpha = -\alpha \implies \alpha = 1$
                    \item $1 - 2\alpha = 0 \implies \alpha = \frac 12$
                    \item $-\alpha = 0 \implies \alpha = 0$ (no rescaling)
                \end{enumerate}

                \begin{enumerate}
                    \item $1 - 2 \alpha = -\alpha \implies \alpha = 1$
                        \begin{align*}
                            \varepsilon ^{-1} Y''(z) + \varepsilon^{-1}2Y'(z) +
                            2Y(z) &= 0 \implies  \\
                            Y''(z) + 2Y'(zz) + \varepsilon2Y(z) &= 0
                        \end{align*}
                            The two largest terms (with
                            $\varepsilon^{-1}$) can cancel out and the other is a
                        small correction

                        So this one is okay.
                            
                    \item $1 - 2\alpha = 0 \implies \alpha = \frac 12$

                        $Y''(z ) + \varepsilon^{-\frac 12} 2Y'(z) + 2Y(z)$

                        The largest term (the middle term) does not have another
                        term of the same order which can balance it out to reach
                        a possible solution, so this cannot produce a solution. 
                \end{enumerate}

            \item  So we choose $\alpha = 1$, and can now sovle with
                perterbation.

                \emph{Boundary condition:} We only expect this to give a good
                solution near $x=$, not for the whole domain, so we can only use the boundary
                condition at $x = 0$.

                $\alpha = 1 \implies Y''(z) + 2Y'(z) + \varepsilon2Y(z) =
                0,\quad Y(0) = 0$ 

                \begin{align*}
                    Y(z) &= Y_0(z) + \varepsilon Y_1(z) + \cdots \\
                    0 &=  (Y_0(z) + \varepsilon Y_1(z) + \dots)'' + 2(Y_0(z) +
                    \varepsilon Y_1(z) + \dots)' + \varepsilon(Y_0(z) +
                    \varepsilon Y_1(z) + \dots) \\
                \end{align*}
                \begin{itemize}
                    \item [$e^0$] We get $Y_0'' + 2Y_0'(z) = 0$, which is
                        solvable:

                        $Y_0(z) = e^{\lambda z} \implies \lambda^2 + 2\lambda =
                        0 \implies \lambda_1 = 0, \lambda_2 = -2$ so

                        $y_0(z) = c(1 - e^{-2z})$, with an unknown constant $c$.

                        Transforming back: $x = \varepsilon z \implies$

                        \marginnote{The boundary layer solution, which is only
                        valid around $x = 0$.}
                        $$y_0(x) =  c(1 - e^{\frac{-2x} \varepsilon})$$
                \end{itemize}

            \item Matching the two solutions: 

                Boundary layer (sometimes called inner solution), valid near $x = 0$. $y_{BL}(x) =  c(1 - e^{\frac{-2x}
                \varepsilon})$. 

                \marginnote{$y_{out}$ obtained before using normal perterbation}
                Outer solution: $y_{out} = e^{1 - x}$. Valid away from $x = 0$.

                We want to combine them to crate an approximation that is valid
                for the whole domain. 

                Find the \textbf{common limit}:

                $y_{BL} (x, \varepsilon \to 0) = c$

                $Y_{out}(x \to 0)  = e$.

                This implies $e^1 = c$, $y_{BL} = e(1 - e^{\frac{-2x}{\varepsilon}})$

            \item Combine to get composite solution:
                \begin{align*}
                    y(x) &= y(x)_{OUT} + y(x)_{BL} - \text{common limit} \\  
                        &=  y(x) e^{1-x} + e(1 - e^{-2x/\varepsilon}) - e \\
                        &= e^{1-x} - e^{1 - \frac {2x}e}
                \end{align*}


        \end{itemize}

        By solving directly we know the exact solution for $\varepsilon = 0.1$ 
        is 

        $-2.874 \left(e^{-18.944 x}-e^{-1.056 x}\right)$

        \begin{center}
            \begin{tikzpicture}{testp}
        \begin{axis}[width=0.9\linewidth,height=2in,xlabel = $x$, ylabel=$y$, xtick={0,0.5}]
            \addplot[domain=0:0.5, samples=200,color=red]{exp(1-x) };
            \addlegendentry{$y_{out}(x)$}
                \addplot[domain=0:0.5, samples=200,color=blue]{ e * (1 -
                        exp((-2 * x) /
                0.1))};
            \addlegendentry{$y_{BL}(x)$}
            \addplot[domain=0:0.5, samples=200,color=green]{exp(1-x) -  exp(1 -
                2 * x / 0.1)};
            \addlegendentry{Combined}
                \addplot[domain=0:0.5, samples=200,color=orange]{-2.71964267074 *
                    (exp(-18.994 * x) -  exp(-1.0553 * x ) };
            \addlegendentry{Exact}
            \end{axis}
        \end{tikzpicture}
        \end{center}


        \lecture{}

        Mathematica.

        \lecture{}

        \subsubsection{Methods of Multiple Scales} 

        \begin{ex}
            Nonlinear pendulum without damping: 

            $ml\theta'' = -mg \sin(\theta)$

            In nondimensiona form: 
            $\theta '' = \sin(\theta),\quad\theta(t) = ?$

            To do perterbation we add an assumption, $\theta(0) = \varepsilon$,
            $\theta'(0) = 0$, and a small initial angle $\epsilon \ll 1$

            Trying regular perterbation: $\theta(t) = \theta_0 (t) + \varepsilon
            \theta_1(t) + \cdots. $
            Use taylor expansion of $\sin(\theta)$ around $0$ for small
            oscillations. $\sin(\theta) \approx \theta - \frac {\theta^3} 6 +
            \cdots $ So when oscillation is small you can use the taylor series to
            represent $\sin(\theta)$.

            Sub the expansion in the initial conditions: $\theta_0(0) +
            \varepsilon \theta_1(0) + ... = \varepsilon. $ so $\theta_0(0) = 0$,
            $\theta_1(0) = 1$. And 

            $\theta_0'(0) + \varepsilon\theta_1'(0) + ...  = 0 \implies
            \theta'_0(0) = 0, \theta_1'(0) = 0$.

            \begin{align*}
                \varepsilon^0 :&& \theta_0'' = -\theta_0' + \frac 16
                (\theta_0^3) + \cdots = -\sin(\theta_0) \\
                              &&\theta_0(0) = 0\quad \theta_0'(0) =  0 \\
                              &&\theta_0(t) = 0 \\
                \varepsilon^1:&& \theta_1'' = -\theta_1  \\
                              && \theta_1(0) = 1, \theta_1'(0) = 0 \\
                              \theta_1(t) = \cos(t) \implies \theta(t) =
                              \varepsilon \cos(t)
            \end{align*}

            The approximation is not uniform in time. 

            Solving for the next power of $\varepsilon$...

            $$\theta(t) = \varepsilon\cos(t) - \frac 1 {16} \varepsilon^3 t \sin(t) +
            \cdots$$

            For very large $t$ the second term may be larger than the first
            which is a problem since we expect the order to be based only on the
            power of $\varepsilon$.

        \begin{center}
            \begin{tikzpicture}{testp}
        \begin{axis}[width=0.9\linewidth,height=2in,xlabel = $t$,
            xtick={180,190,200}]
            \addplot[domain=0:50, samples=200,color=blue]{  exp(-0.1/(2*x)) *
                sin(1/2 * sqrt(4 - 0.1^2) * deg(x)) };
            \addlegendentry{Exact}
            \addplot[domain=0:50, samples=200,color=red]{  (1 - 0.1*x/2) *
                sin(deg(x)) };
            \addlegendentry{Approx}
            \end{axis}
        \end{tikzpicture}
        \end{center}

        There are two timescales in the problem the oscillation period and the
        slow drift in the phase of the oscillation which is not described  by
        the first term so we can add parameters to the time to split it into
        multiple timescales. 

        \end{ex}

        \week{}
        \lecture{}

        The ordering of the terms, their relative magnitude, may change in a
        way which violates our assumptions. This is what the method of multiple
        scales addresses. 

       There are two characteristic timescales: oscillation period,
                the time of decay.

        \subsubsection{Method of multiple scales}

        \begin{itemize}
            \marginnote{(the equation is already non-dimensional)}
            \item Introduce two new nondimensional time variables to rescale.

                $t_1 = t$ and $t_2 = \varepsilon^\alpha t$ so $y(t) \implies Y(t_1,t_2)$
                
            \item Use the chain rule:

                $$\frac d{dt} = \frac {\partial}{\partial t_1} 
                + \varepsilon^\alpha \frac \partial{\partial t_2}$$

                $$\frac{d^2}{dt^2} = \frac{\partial^2}{\partial {t_1}^2} +
                2\varepsilon^\alpha \frac {\partial^2}{\partial t_1\partial t_2} +
            \varepsilon^{2\alpha} \frac {\partial^2}{\partial {t_2}^2}$$

        \item Initial conditions:

            $Y(0,0) = 0, \quad \frac{\partial Y}{\partial t_1}(0,0) +
            \varepsilon^\alpha \frac {\partial Y}{\partial t_2}(0,0) = 1$
        \item Substitute

            $$Y(t_1,t_2) = Y_0(t_1,t_2) + \varepsilon Y_1(t_1,t_2) +
            O(\varepsilon^2)$$
        \item

            Largest:$\varepsilon^0: \frac{\partial^2 Y_0}{\partial t_1^2} + Y_0 =
            0$

            $Y_0(0,0) = 0, \quad \frac{\partial Y_0} {\partial t_1}(0,0) = 1$

            We know there is a general solution: 

            $$Y_0(t_1, t_2) = A(t_2) \sin(t_1) + B(t_2)\cos(t_2)$$

            Where $B(0) = 0$ and $A(0) = 1$ is the initial condition, by
            susbstituting the initial conditions in the general solution. 

        We will choose $A(t)$ and $B(t)$ in such a way as to avoid
        resonancnes. 


            In order for the friction terms to be balanced by the mixed
            derivatives choose $\alpha = 1$, (slow time: $t_2 = \varepsilon^1
            t$) 


            Terms $\varepsilon^1$: 

            $\frac {\partial^2 Y_1} {\partial t_1^2} + Y_1 = - \frac {\partial
            y_0} {\partial t_1} - 2 \frac {\partial ^2 Y_0} {\partial t_1
        \partial t_2}$ 

        Substitute

        $\frac {\partial^2 Y_1} {\partial t_1^2} + Y_1 = - A(t_2)\cos(t_1) +
        B(t_2) \sin(t_1) - 2A'\cos(t_1) + 2B'(t_2) \sin(t_1)$

        Group together forcing terms.

        $\frac {\partial^2 Y_1} {\partial t_1^2} + Y_1 = [- A(t_2) + 2A'(t_2)]\cos(t_1) +
        [B(t_2) \sin(t_1) +  2B'(t_2)] \sin(t_1)$

        To avoid resonant forcing which leads to secular terms we need $A + 2A'
        = 0$ and $B + 2B' = 0$. So we need to solve these differential
        equations:

        $$A(t_2) = A(0) e^{\frac{-t_2}2} \quad \quad B(t_2) = B(0)
        e^{\frac{-t_2}2}$$

    \item Transform back and substitute

        $y(t) = e^{\frac {-\varepsilon t}2 }   \sin(t)$

        \end{itemize}


\chapter{Traffic Flow Models}
        \lecture{}

        \textbf{Agent based models:} rules for individual cars and how they
        interact can be stochastic. can be used for simulations, not suitable
        for mathematical analysis.

        \textbf{Conntinuum model:} car density function and velocity
        distribution etc. Treated as a continuous variable.

        \begin{definition}{}{} Density function (1D)

            $\rho(x,t) = \frac {\text{The number of cars in } [x - \Delta x, x + \Delta x]} {
        \text{distance } 2\Delta x}$

        The $\Delta x$ needs to be large compared to the size of the cars but
        small compared to the region.
        \end{definition}

        \begin{definition}{}{} Flux $J(x,t)$

            The number of cars that passing a position $x$ per unit time
        $\frac {[t-\Delta t, t + \Delta t]}  {2\Delta t}$. 

        For constant velocity of cars with distantce $d$ $J = \frac v {l + d} =
        v \rho$.
        \end{definition}

        \begin{definition}{}{} Balance Law

            Conservation law for cars along the region $[x - \Delta x, x +
            \Delta x]$ of the road. 

            Initial cars + cars going in - cars going out. 

            $$J(x - \Delta x, t) \Delta t - J(x + \Delta x , t) \Delta t =
            \frac{J(x - \Delta x, t) - J( x + \Delta x, t)}{2\Delta x}$$
            
            Taking the limit $\Delta t \to 0$ and $\Delta x \to 0$

            $$\frac {\partial \rho}{\partial t} = - \frac {\partial J}{\partial
            x}$$

            The flux can also be written as $J = \rho v$ so we have the equation

            $$\frac {\partial \rho}{\partial t} + \frac {\partial (\rho
            v)}{\partial x} = 0$$
        \end{definition}

        In order to obtain a closed system we need a relationship between the
        velocity and density of traffic.

        \begin{definition}{}{}
            Constitutive law

        $v(\rho)$ is known as the constitutive law

        It is natural to expect $v(\rho)$  is monotonously decreasing: as the
        density of cars increases the traffic slows down. And $v(0) = v_M$ when
        there are no other cars around traffic moves at the maximum allowable
        velocity. 

        This can be based on observations from real traffic. 
        \end{definition}

        The simplest constitutive law is greenshield's law assumes a linear
        function:

        $$v(\rho) = v_M\left(1 - \frac \rho{\rho_M}\right)$$

        substituting into the balance law this gives 

        $$\frac {\partial \rho}{\partial t} + c(\rho) \frac {\partial (\rho
            )}{\partial x} = 0 \quad\quad \text{where} $$


            $$c(\rho) = v_M\left(1 - \frac {2\rho}{\rho_M}\right)$$

        If you assume an almost uniform density and they propogate like waves
        you can get wave velocity, which is the velocity of propogation of 
        perterbations in the density of traffic:

        $$c(\rho) = v(\rho) + \rho v'(\rho)'$$

        $$\frac {\partial \rho} { \partial t} + c(\rho) \frac {\partial \rho}{\partial x} =
        0$$

        \subsubsection{Solving}

        Consider the simplest case, which is consntant velocity: $v(\rho) = a$.
        But we don't assume density is uniform: $v$ is independent of the density, 
        which is not realistic for traffic it
        is just a simpler maths problem; advection equation for fluid flow. 

        The solution is $\rho(x,t) = f(x - at)$, can check this with
        substitution. 

        $$\frac {\partial \rho}{\partial t} + a \frac {\partial \rho }{\partial x}
        = 0$$

        Initial: $\rho (x,0) = f(x) $

        Boundary: $\rho (x = 0, t) = g(t)$.

        Don't need boundary conditions on both sides of the domain, because it
        only goes in one direction. 

        Solution : $$\rho(x,t)  = \begin{cases}f(x - at) & x > at \\ g(t - \frac
        )xa & x < at\end{cases}$$

        Nonconstant velocity:

        $$\frac {\partial \rho}{\partial t} + c (\rho) \frac {\partial
        \rho}{\partial x} = 0, \quad \rho(x,t=0) = f(x)$$ 

        The density changes according to the velocity. 

\week{}

\lecture{}

\begin{ex}
\end{ex}

\begin{itemize}
    \item uniform positive velocity
    \item  density linearly linked to velocity according to greenfields law
\end{itemize}



Perterbations propagate along the roads according to wave equation:
areformation of hte balance equation:

$$\frac {\partial \rho}{\partial t} + c (\rho) \frac {\partial \rho}{\partial x}
= 0$$

Where $c(\rho) = J' (\rho) = v + \rho v'(\rho)$.

\textbf{Special case with constatnt velocity} gives the travelling wave
solution: $\rho(x,t) = f(x - ct)$ the density distribution is transported along hte x axis with
        constatnt velocity without any distortion. The shape is exactly the
        same. This is when the cars dont care about the density. 

        \textbf{Method of characteristic trajectories}: $x(t) = x_0 + ct$. 
find th trajectory ending at $(x,t) \implies x_0 = x(t) - ct$; 
  The density is conserved along the
        lines $x(t) = x_0 + c(\rho(x_0))t$.  (A transformation into an arbitrary
        space such that there are lies where $\rho$ is constant to make it easy to
        solve.)                            


The same method works for non-constant velocity, we know the waves propogate at
constant speed (because of our assumption) but we don't know the slope. 

For any $(x,t)$ find an $x_0$  such that $x = x_0 + c(f(x_0))t$ then $\rho(x,t)
= \rho(x_0, 0) = f(x_0)$. 

\emph{characteristics:} points $x_0$ that move to the point $(x,t)$.

\begin{itemize}
    \item There may be regions with no characteristics
    \item there may be regions with multipl eintersecing characteristics
\end{itemize}

\section{Rankine-Hugoniot condition:}

\begin{itemize}
    \item If there are discontinities in the initial conditions of density
        ($\rho$), you can
        end up with diverging or overlapping characteristic lines and there can
        be no $\rho_l < \rho_r$ or no $\rho_l > \rho_r$ unique solutions for a particular region. 
\end{itemize}

        \begin{center}
            \begin{tikzpicture}{testp}
                \begin{axis}[xmin=0, xmax=10, ymin=0, ymax=10, ticks=none, xlabel=x, ylabel=t]
                    \addplot[domain=0:7.5,samples=90,color=red]{ 2.5 * (x - 5)};
                \addplot[domain=0:20,samples=90,color=black,very thick,dotted]{ 2 * x};
                \addplot[domain=0:20,samples=90,color=black,very thick,dotted]{ 2 * (x - 1)};
                \addplot[domain=0:20,samples=90,color=black,very thick,dotted]{ 2 * (x - 2)};
                \addplot[domain=0:20,samples=90,color=black,very thick,dotted]{ 2 * (x - 3)};
                \addplot[domain=0:20,samples=90,color=black,very thick,dotted]{ 2 * (x - 4)};
                \addplot[domain=0:20,samples=90,color=black,very thick,dotted]{ 2 * (x - 5)};
                \addplot[domain=0:20,samples=90,color=black,very thick,dotted]{ 3 * (x - 5)};
                \addplot[domain=0:20,samples=90,color=black,very thick,dotted]{ 3 * (x - 6)};
                \addplot[domain=0:20,samples=90,color=black,very thick,dotted]{ 3 * (x - 7)};
                \addplot[domain=0:20,samples=90,color=black,very thick,dotted]{ 3 * (x - 8)};
                \addplot[domain=0:20,samples=90,color=black,very thick,dotted]{ 3 * (x - 9)};
            \end{axis}
        \end{tikzpicture}
        \end{center}


You can integrate the balance equation  to get information about the
discontinuity. Integrate $x$ in $[s(t) - \varepsilon, s(t) + \varepsilon]$,
givin

$$\int_{s - \varepsilon}^{s + \varepsilon} \frac {\partial \rho}{\partial t} dx
+ J(\rho_r) - J(\rho_l) = 0$$

Can differentiate this again using \textbf{Leibniz rule}: the moving integrating bounds
changes the value of the integral. 

If they are constant:

$$\frac d {dt} \int_\alpha ^\beta f(x,t) dx = \int_\alpha^\beta \frac {\partial
f} {\partial t} dx$$

Otherwise 

$$\cdots$$

Specifically for the non-constant functions $\alpha(t) = s - \varepsilon$ and $\beta(t) = s +
\varepsilon$

$$\frac d {dt} \int_{x - \varepsilon}^{s + \varepsilon} \rho dx - s'(t) \rho_r +
s'(t) + \rho_l + J(\rho_r) - J(\rho_l) = 0$$


$$\frac d {dt} \int_{x - \varepsilon}^{s + \varepsilon} \rho dx = \frac d {dt}
(\epsilon \rho_l + \epsilon \rho _r) = 0$$

\marginnote{cool}
$$s'(t) = \frac {J(\rho_r) - J(\rho l)} {\rho_r - \rho_L} \quad\quad\quad(*)$$

From this we can also show that the discontinuity moves with the velocity 

$$s'(t) = \frac 1 {\rho_r - \rho _l } \int_{\rho_l} ^{\rho_r} c(\rho) d\rho$$


\begin{itemize}
    \item This solution is not limited to greenshield's law it gives
        caharacteristic lines for intersections of the characteristic lines and
        they are related to the function $c(\rho)$ and the way these boundaries
        move through traffic (because of hte moving $s$ in the integral before). 
\end{itemize}

Under \textbf{Greenshield's law} it becomes the average speed of the wave velocity:

$$s'(t) = \frac 12 (c_R + c_L)$$


We want to solve these kinds of initial conditions cause the models often
produce steps (singularities in the solution) or 'shocks' between fast and slow
moving regions (because the faster cars `run into' the slow cars and teh
slowdown region gets smaller until it is a jump). 

note that $c(\rho)$ is not strictly positive: perturbation waves can move
backwards through the traffic, eg if a car suddenly stops on a highway. 


        \lecture{}

        \textbf{Initial Conditions $\rho_l > \rho_r$}

        \begin{center}
            \begin{tikzpicture}{testp}
                \begin{axis}[xmin=0, xmax=10, ymin=0, ymax=10, ticks=none, xlabel=x, ylabel=t]
                \addplot[domain=0:20,samples=90,color=black,very thick,dotted]{ 3 * x};
                \addplot[domain=0:20,samples=90,color=black,very thick,dotted]{ 3 * (x - 1)};
                \addplot[domain=0:20,samples=90,color=black,very thick,dotted]{ 3 * (x - 2)};
                \addplot[domain=0:20,samples=90,color=black,very thick,dotted]{ 3 * (x - 3)};
                \addplot[domain=0:20,samples=90,color=black,very thick,dotted]{ 3 * (x - 4)};
                \addplot[domain=0:20,samples=90,color=black,very thick,dotted]{ 3 * (x - 5)};
                \addplot[domain=0:20,samples=90,color=black,very thick,dotted]{ 2 * (x - 5)};
                \addplot[domain=0:20,samples=90,color=black,very thick,dotted]{ 2 * (x - 6)};
                \addplot[domain=0:20,samples=90,color=black,very thick,dotted]{ 2 * (x - 7)};
                \addplot[domain=0:20,samples=90,color=black,very thick,dotted]{ 2 * (x - 8)};
                \addplot[domain=0:20,samples=90,color=black,very thick,dotted]{ 2 * (x - 9)};
            \end{axis}
        \end{tikzpicture}
        \end{center}

        The initial jump is dissolved into a transition zone that gets wider
        over time: ``expansion fan''. In this case the Rankine-Huginot
        condition cant be used, since the jump discontinuity doesnt persist for
        greater $t$. 

        \begin{ex} Left side maximum density and zero velocity, right side
            maximum velocity and zero density. 

            \begin{itemize}
                \item Smooth out the initial step function with a linear interpolation 

                    $$f(x) = \begin{cases} \rho_M & x < -\varepsilon \\ 
                        \rho_M \frac{\varepsilon - x}{2\varepsilon} & -\varepsilon < x <
                \varepsilon \\
                        0 & \varepsilon < x
                \end{cases}$$

                In the transition zone you have $-\varepsilon < x < \varepsilon$
                and greenshields law to calculate wave velocity $c(f(x)) = V_M \frac x
                \varepsilon$. 

                solve stuff
            \end{itemize}
            
        \end{ex}

        \week{}
        \lecture{}
        \chapter{Continuum Mechanics}

        Describing the movement and deformation dymaics of materials represented
        as a continuous medium (rather than representations as connected
        particles). For example elastic solids, flowing gases and liquids. 


        \incfig{ref}


        \subsection{Coordinates}

        Consider two different coordinate systems: 
\begin{itemize} \item [] \textbf{Material coordinates:} label a pooint within the
                continous material and follow its movement along hte
                trajectories $X(A,t)$ where $A$ is the intial location: $X(A,0)
                = A$. (If you painted a dot on a piece of jelly).
            \item [] \textbf{Spatial coordinates:} Fixed coordinate system $x$
                of an external observer. Does not follow the movement of the
                material.
        \end{itemize}

        \textbf{Material Coordinates:}

        \begin{itemize}
            \item Trajectory $X(A,t)$ where $X(A,0) = A$

            \item Displacement: $U(A,t) = X(A,t) - A$

            \item Velocity: $V(A,t) = \frac {\partial U(A,t)} {\partial t}$
        \end{itemize}

        \textbf{Spatial Coordinates:}

        \begin{itemize}
            \item Displacement: $u(x,t)$ of the material $x$ at time $t$.
                To relate it to $U(A,t) $ find $A$ at $t = 0$ by solving $X(a,t) = x$ for
                $A = a(x,t)$  (a in the spatial coordinates).

            \item $u (x,t) = U(a(x,t) , t)$
            \item $v(x,t) =  \frac {\partial U(a(x,t),t)} {\partial t} =  V(a(x,t),t)$ which is not equal to $\frac {\partial
                u}{\partial t}$ (there is an additional dependence on time, we
                evaluate by traversing the trajectory function $A$).
        \end{itemize}      

        Example: 

        Consider a material described by the trajectories $X(A,t) = A \frac{1 +
        2t}{1 + t}$. Initial domain: $0 < A < l_0$.

        The displacement from $A$ of a material point at time $t$ is

        $$U(A,t) = X (A,t) - A = A \frac t {1 + t}$$

        \marginnote{Can check this is different from $\frac{\partial u}{\partial
        t}$. However it is the same as $\frac {\partial X}{\partial t}$ since $A$
        is a constant.}

        The velocity: 
        $$V(A,t) = \frac {\partial U(A,t)}{\partial t} = A \frac 1{(1 + t)^2}$$

        The domain in material coordinates is $0 < A < l_0$, in spatial
        coordinates $X(0,t) = 0 < x < X(l_0,t) = l_0 \frac {1 + 2t}{1 + t}$

        If final coordinate is $x$ at $t$ find the initial coordinate $A$

        $$x = X(A,t) \implies A = x \frac {1 + t}{1 + 2t}$$

        The displacement of a material point at $x$ at $t$ relative to its
        initial position: 


        $$u(x,t) = U(x\frac{1 + 2t}{1 + t}, t) = x \frac t {1 + 2t}$$

        $$v(x,t) = V(x \frac {1 + 2t}{1 + t}, t) = x \frac 1 {(1 + t)(1 + 2t)}$$

        \subsection{Material Derivative}

        Consider a scalar quantity that describes a propertiy transported by the
        material. In spatial coordinates the distribution is described by
        $f(x,t) $ and in material coordinates is $F(A,t)$. They must satisfy the
        relationship $f(X(A,t),t) = F(A,t)$.

        \emph{What is the rate of change of this property in a material
        coordinate system (along the trajectories of movement of the material
        points).}

        \begin{align*}
            \frac {\partial F}{\partial t} &= \frac {\partial f}{\partial t} +
        \frac {\partial f} {\partial x} \frac {\partial X(A,t)} {\partial t} \\
                                           &=
        \frac {\partial f}{\partial t} + V(A,t) \frac {\partial f}{\partial x}
                                        \\ &=
        \frac {\partial f}{\partial t} + v(x,t) \frac {\partial f} {\partial x}
        \end{align*}

        Material Derivative: $$\frac D {Dt} \equiv \frac \partial{\partial t} +
        v \frac \partial {\p x} \equiv \frac {\partial F}{\partial t}$$

        \lecture{}

        \subsection{Continuity Equation}

        The law of conservation of mass 

        $$\frac {\partial \rho}{\partial t} + \frac {\partial(v \rho)}{\partial
        x} = 0$$

        This is not sufficient to solve, we need another equation to describe
        the velocity field $v(x,t)$. We use a continuum mechanics version of
        newtons law: conservation of momentum.

        \subsection{Momentum Equation}

        The motion of a point particle under forces: $ma = \sum F_i$.
        The rate of change of momentum:

        $$\frac d {dt} mv = F$$

        Typoes of forces in continuum mechanics

        \begin{itemize}
            \item external surface forces
                \begin{itemize}
                    \item Act along the boundary of the domain
                \end{itemize}
            \item external body forces
                \begin{itemize}
                    \item Act on every material point inside the domain, for
                        example gravity. 
                \end{itemize}
            \item internal forces
                \begin{itemize}
                    \item Stress ($\tau(x,t) = $ force / area) generated by different parts of
                        the material acting on eachother. For example restoring
                        force in elastic material.
                \end{itemize}
        \end{itemize}


        Consider a segment of material between $A_L$ and $A_R$ with trajectories
        $X(A_L,t) = \alpha (t) $ and $X(A_R,t) = \beta(t)$.

        The rate of change of the momentum of this segment is the sum of the
        forces acting on it, external body force and internal surface force. 
        \marginnote{$S$ is the surface force. }

        $$\frac d {dt} \int_{\alpha(t)}^{\beta (t)} S \rho v \; dx =
        \int_{\alpha(t)}^{\beta(t)} S\rho f \; dx + S\tau(\beta,t) - S
        \tau(\alpha, t)$$

        $$\frac d {dt}\int_{\alpha(t)}^{\beta(t)} \rho v \; dx =
        \int_{\alpha(t)}^{\beta(t)} \rho f \; dx + \tau(\beta,t) - \tau(\alpha,
        t)$$

        $$\beta'(t)\rho(\beta9t))v(\beta(t)) -
        \alpha'(t)\rho(\alpha(t))v(\alpha(t)) = \int_{\alpha(t)}^{\beta(t)} \rho
        f \; dx +\int_{\alpha(t)}^{\beta(t)} \frac {\partial \tau}{\partial x}\;
        dx$$

        $\cdots$

        $$\frac {\partial(\rho v)}{\partial t} + \frac {\partial (\rho
        v^2)}{\partial x} = \rho f + \frac {\partial \tau}{\partial x}$$

        $$\rho v_t + v \rho_t + \rho v_x  + v \frac {\partial \rho v}{\partial
        x} = \rho f +_ \frac {\partial \tau} {\partial x}$$

        $\cdots$

        Then you get the momentum equation as:

        $$\rho \frac {Dv}{Dt} = \rho f + \frac {\partial \tau}{\partial x}$$

        \begin{center}
        Material Derivative $=$ body force $+$ internal stress
        \end{center}

        You can think of this as a generalisation of Newton's law with internal
        forces in addition to the external forces. 

        \subsubsection{In Material Coordinates}

        \begin{itemize}
            \item Density $R(A,t) = R_0(A) \left(1 + \frac {\partial U}{\p
                A}\right)^{-1} \quad (*)$
            \item Velocity $V(A,t)$
            \item Stress $T(A,t)$
            \item Body force $F(A,t)$ 
        \end{itemize}

        All considered along the trajectories $X(A,t)$.

        Using the transformation 

        $$\frac {\partial \tau}{\partial x} = \frac{\partial T}{\partial
        A}\left(\frac {\p U}{\p A} + 1\right)^{-1}$$

        We need to apply this to the stress, and momentum equation.

        \begin{align*}
        R(A,t) \frac {\p V}{\p t} &= RF + \frac {\p T}{\p A}\left(\frac {\p
            U}{\p A} + 1\right)^{-1} \\
        R \left(\frac {\p
                U}{\p A} + 1\right)\frac {\p V}{\p t} &= R\left(\frac {\p
            U}{\p A} + 1\right)F + \frac {\p T}{\p A} 
            \\  
            \marginnote{By the continuum equation in material coordinates $(*)$}
            \\
            R_0 \frac {\p V}{\p t} &= R_0 F + \frac {\p T}{\p A}
        \end{align*}

        \lecture{}

        \subsection{Constsitutive Laws}

        \subsubsection{Constitutive Law for Elastic Materials}

        Here, stress ($T$) is generated as a result of the displacement of the
        material elemeents $U(A,t)$. However $T \ne T(U)$ since then a constant
        displacement does not generate zero stress. Stress is generated by the
        relative elongation of material elements: \emph{strain}.

        \subsubsection*{Strain in 1D}

        Consider a small section $A, A + \Delta$ at $t = 0$. At time $t$, 

        \begin{align*}
            A &\to A + U(A, t) \\
            A + \Delta &\to A + \Delta  + U(A + \Delta,t)
        \end{align*}

        After substituting the relative elongation (local strain) becomes

        $$\frac {l(t) - l_0}{l_0} = \frac {U(A + \Delta, t) - U(A,t)}{\Delta}$$

        Then since $\Delta \ll 1$ we can use ethe taylor series approximation
        for $U(A + \Delta, t)$,  $U(A + \Delta, t) \approx U(A) + \frac {\p
        U}{\p A} \Delta + \cdots$

        $$\frac {l(t) - l_0}{l_0} = \frac {U(A + \Delta, t) - U(A,t)}{\Delta} =
        \frac {\p U}{\p A}$$

        Hence in elastic materials the stress is a function of strain

        $$T = f(\frac {\p U}{\p A})$$

        $f(0) = 0$, but it is otherwise determined by the properties of the
        material.

        Assuming a \textbf{linear elastic material}: $T = E\frac {\p U}{\p A}$, where $E$ is a
            constant.
            \marginnote{We can often choose a region around which a linear
            approximation is a good solution}

        Then we get the momentum equation in this form;

        $$R_0 \frac {\p^2 U}{\p t^2} = R_0 F + E \frac {\p^2 U}{\p A^2}$$

        \emph{Note:} we still don't know about the external forces $F$. 

        \begin{ex}
            Consider an elastic bungie cord with length $l_0$: how long is it
            when stretched by its own weight hanging free in vertical position.

            Initial position $0 < A < l_0$ at $t = 0$.

            Boundary conditions: $U(0,t) = 0$ at the fixed end.

            $T(l_0, t) = 0$ at the free end, $\implies$ $E \frac {\p U}{\p A}(A = l_0, t) = 0$
            
            $U(0) = 0$, $U'(l_0) = 0$. 

            $$U(A) = - \frac {R_0 g}{2E} + C_1 A$$

            $$U'(A) = - \frac {R_0 g}E A + C_1$$

            $$C_1 = \frac {R_0 g}E l_0$$

            $$U(A) = \frac {R_0 g}{2 E} A(2 l_0 - A)$$

        final length is $l = l_0 + U(l_0) = l_0 + \frac {R_0 g} {2E} {l_0}^2$

        We can also calculate the stress distribution along the cord:

        $$T(A) = E U'(A) = R_0 g(l_0 - A)$$ 

        The density distribution (non-constant because it is now stretyched.) 

        $$R(A) = \frac {R_0} {1 + \frac {R_0 g}E (l_0 - A)}$$

        Note we are working in material coordinates, if we want them as a
        function of spatial coordinates you need to invert $x = A + U(A)$.

        \end{ex}


        \week{}
        \lecture{}

\section{Review of Continuum Mechanics in 1D}

        \incfig{ref}


        \subsubsection*{Material coordinate}

        \begin{itemize}
            \item $t$ time
            \item $A$ cross section
            \item $X(A,t)$ a spatial position
            \item 
                $U(A,t) = X(A,t) - X(A,0)$ a displacement at a time from the
            \item 
                $V = \p_t U(A,t) $ velocity

            \item $R(A,t) = \rho(X(A,t), t)$ density
            \item $T(A,t) = \tau(X(A,t), t)$ stress
            \item $F(A,t) = f(X(A,t), t)$ body force
        \end{itemize}

        \marginnote{bad handwriting idk if these are right}
        \subsubsection*{Spatial coordinate}
        \begin{itemize}
            \item $t$ time
            \item $a(x,t) = X^{-1}$ 
            \item $x$ spatial position
            \item $u(X,t) = U(a(x,t) ,t)$ displacement
            \item $v = D_t u = \frac {\partial_t u}{1 - u_X}$ 

            \item $\rho(x,t)$ density ($\frac {mass}{volume}$)
            \item $\tau(x,t)$ stress ($\frac {force}{area}$)
            \item $f(x,t)$ body force ($\frac {force}{mass}$)
        \end{itemize}

        \subsubsection*{Governing equations}


        $$\p_t f + \p _x(v \rho) = v \quad{\text{continuity equation}}$$

        $$\rho(\p_t v + v v_x) = \rho f + \p _X \tau = D_t v
        \quad{\text{momentum eqn}}$$

        Have the unknowns $R, V,T$ and only two equations so need a
        constitutive equation for $\tau$ or $T$. This can be found using a creep
        experiment: stretch a material an amount and see how much force is
        requried, or apply an amount of force and see how much it streteches. 

        $$\frac{\ell} {\ell_0} = \lambda \quad\quad\text{Extension ratio}$$

        \begin{itemize}
            \item Typically extension grows with stress.
        \end{itemize}

        At steady state $\p_{t t} U = 0 \implies \p _A T = 0$ stress $T$ is
        constant along the material (from the momentum equation) so we cannot have $T = T(U)$ stress is a
        function of displacement.


        $$
\begin{array}{|l|l|l|}
\hline \text { Name } & \text { Ratio } & \text { Definition } \\
\hline \text { Lagrangian Strain } & \left(\ell-\ell_{0}\right) / \ell_{0} & \epsilon=U_{A} \\
\hline \text { Eulerian Strain } & \left(\ell-\ell_{0}\right) / \ell & \epsilon_{e}=u_{x} \\
\hline \text { Green Strain } & \left(\ell^{2}-\ell_{0}^{2}\right) /\left(2 \ell_{0}^{2}\right) & \epsilon_{g}=U_{A}+\frac{1}{2} U_{A}^{2} \\
\hline \text { Almansi Strain } & \left(\ell^{2}-\ell_{0}^{2}\right) /\left(2 \ell^{2}\right) & \epsilon_{a}=u_{x}-\frac{1}{2} u_{x}^{2} \\
\hline \text { Midpoint Strain } & 2\left(\ell-\ell_{0}\right) /\left(\ell+\ell_{0}\right) & \epsilon_{m}=U_{A} /\left(1+\frac{1}{2} U_{A}\right) \\
\hline \text { Hencky Strain } & \ln \left(\ell / \ell_{0}\right) & \epsilon_{h}=\ln \left(1+U_{A}\right) \\
\hline
\end{array}
$$


\lecture{}
Eg For lagrangian strain: $\epsilon \sim \frac {\ell - \ell_0}{\ell_0}$

\begin{align*}
    \frac{l - l_0}{l_0} &= \frac {X(A + \Delta A, t) - X(A - \Delta t, t) -
    2\Delta A - A \times A}{2 \Delta A} \\
                        &= \frac {X(A + \Delta A, t) - (A + \Delta A) - (X(A -
                        \Delta t, t) - (A - \Delta A))}{2\Delta A} \\
                        & \quad\quad \text{Taylor expansion} \\
                        &= \frac {U(A,t) + \Delta A U_A(A,t) + O(\Delta A^2) - (U(A,t)
                    - \Delta A(U_A(A,t) + O(\Delta A)^2)) } {2\Delta A}\\
                        &= U_A(A,t) + O(\Delta A) \\
                        &\implies \lim_{\Delta A \to 0} \frac {\ell -
                        \ell_0}{\ell_0} = U_A(A,t) =: \epsilon \\
\end{align*}

$\implies$ that  for elastic materials 

                        $$T = T(\epsilon) = T(U_A) $$

Now the momentum equation is given by 
\marginnote{(non-linear) wave equation} 

$$R_0 U_{tt} = R_0 F + T'(U_A) U_{AA}$$

$$\lambda = \frac {\ell }{\ell_0} = \frac {\ell - \ell_0 + \ell_0} {\ell_0} \sim
\epsilon + 1 $$

$\implies$ the reference state has $\lambda = 1 \Leftrightarrow \epsilon = 0$

For a small enough $\epsilon$ any material will exhibit a linear stress--strain
relation. 

$$T(\epsilon) = E\epsilon - E\p_A U$$

Wher $E$ is young's modulus also known as the elastic modulus. 

Physical dimensions:

\begin{itemize}
    \item[] $[T] = \frac {force}{area}$ 
    \item[] $[\epsilon] = [\frac {\ell - \ell_0}{\ell_0}] = 1$
    \item[] $[E] = \frac{force}{area}$
\end{itemize}

\begin{center}
\begin{tabular}{l|l|l}
\hline Material & Young's modulus $(\mathrm{GPa})$ & Density $\left(\mathrm{kg} / \mathrm{m}^{3}\right)$ \\
\hline Diamond & 1000 & 3500 \\
\hline Stainless steel & 200 & 8030 \\
\hline Glass & 65 & 2600 \\
\hline White oak & 12 & 770 \\
\hline Beeswax & $0.2$ & 960 \\
\hline Rubber & $0.007$ & 1200 \\
\hline Silica aerogel & $0.001$ & 100 \\
\hline
\end{tabular}
\end{center}

The momentum equation: $R_0 U_{tt} = R_0 F + EU_{AA}$ multiplied by $\frac
1{R_0}$ gives  $U{tt} = F + \frac E {R_0} U_{AA}$, which is the linear constant
coefficient wave equation. So we can understand that the characteristic frequencies 
are proportional to wave speed $c = \sqrt{\frac E{R_0}}$. For example a high
elastic modulus in the body of a brass instrument has a high wave speed so
expresses those frequencies more. 

Another way of imagining this is if the linear elastic material is imagined as a
chain of beads connected by Hookean springs (at reference length $l_0$), then $F = K(\ell - \ell_0)$ where
$K$ is the spring constant. 

\incfig{linearelasticspring}


\begin{align*}
    T &=E \frac {\ell - \ell_0}  {\ell_0} \\
    &\text{multiply by $\sigma$, cross sectional area} \\
    \sigma T &= \frac { \sigma E}{ \ell_0 } (\ell - \ell_0) \\
    &= K 
\end{align*}

Spring constant is $\frac {\sigma E}{\ell_0}$

So we can use \textbf{the creep test} to find young's modulus. 

We need to solve the system of equations 

$$\begin{cases} 0 = EU_{AA} \\ 0 = U(0) \\ U(\ell_0) = \ell - \ell_0 \end{cases}$$

\begin{align*}
    U &= \alpha A = B \\
    U(0) &= 0 = \beta \\
    U(\ell_0) &= \ell - \ell_0 = \alpha\ell_0  \\
    \implies U(A) = \frac {\ell - \ell_0} A
\end{align*}

Material density

\begin{align*}
    R &= \frac {R_0}{1 + U_A} \\
    \implies R &= \frac{R_0} {1 + \frac {\ell - \ell_0} {\ell_0}} = \frac {R_0} {\frac
    \ell {\ell_0}} = \frac {\ell_0}{\ell} R_0  
\end{align*}

\lecture{}

\subsubsection{Bungee rope at steady state}

\begin{itemize}
    \item []
$\ell$ is the length at rest while hanging: extended undner its own weight. 
\item[]
$\ell_0$ is the length of the reference configuration: bungee rope under no
force, at rest
\end{itemize}

$R_0 U_{t t} = R_0 F + EU_{ AA }$

$R_0 U_{t t} = 0$

$R_0 F$ is body force

The body force is obviously gravity so $F = g = 9.81 m/{s^2}$

So the momentum equation is $U = R_0 g + EU_{A A}$

We dont nkow the value of $U$ at $A = \ell_0$ but we know the boundary
condition: $T(\ell_0) = 0 \implies E U_A(\ell_0) = 0$, and $U(0) = 0$

Solve

$$\begin{cases}0 = R_0 g + EU_{A A}\quad\quad(*)\\ U(0) = 0 \\ U_A(A = \ell_0) = 0 \end{cases}$$

THe easiest way is to integrate $(*)$ on $(A, \ell_0)$ to get $0 = R_0 g(\ell_0
- A) + E(U_A(\ell_0)  - U_A(A))$

$\implies U_A = \frac {R_0}E g(\ell_0 - A) \quad \quad (* * )$

Ie, the tension on the bungee cord decreases linearly along the bungee cord
since $U \propto T$. (This makes sence since as you go down the cord there is
less weight beneath each cross section.)

        \begin{center}
            \begin{tikzpicture}{testp}
                \begin{axis}[ytick={0}, xtick={0.1},xlabel=$A$,
                    xticklabels={$\ell_0$},yticklabels={$\frac {R_0}E
                    g\ell_0$},]
                \addplot[domain=0:0.1, samples=90,color=blue,]{ -x / 4};
            \end{axis}
        \end{tikzpicture}
        \end{center}

        Integrate $(* * ) $ on $(0, A)$ to get

        $U = \frac {R_0} E g (\ell_0 A - \frac {A^2}{2})$

        $U(\ell_0) = \frac{R_0 g}{E} \frac {{\ell_0}^2}  2$

        \begin{center}
            \begin{tikzpicture}{hy}
                \begin{axis}[ytick={1}, yticklabels={$\frac{R_0 g}{E} \frac
                {{\ell_0}^2} 2$}, xtick={0,2},xlabel=$A$,
                    xticklabels={0, $2\ell_0$},
                    ]
                    \addplot[domain=0:2, samples=90,color=blue,]{ -(x -1)^2 + 1};
            \end{axis}
        \end{tikzpicture}
        \end{center}

        So clearly $\ell$ grows like ${\ell_0}^2$ because of the cord's own weight.



\section{Hyperelastic Materials}

\begin{itemize}
    \item Nonlinear stress-strain relations. 
    \item Energy balance
\end{itemize}


For example for rubber the stress strain raltion looks something like



\begin{figure}[htpb]
    \centering
    \begin{tikzpicture}{testppop}
        \begin{axis}[ytick={},xtick={},xlabel=$\epsilon$, ylabel=$T$,]
            \addplot[domain=-1:2, samples=90,color=blue,]{ (x -0)^3 + 0};
    \end{axis}
\end{tikzpicture}
    \caption{Stress-Strain for rubber}
\end{figure}

\begin{align*}
    T(\epsilon) &= \psi'(\epsilon) \text{ where }  \\
    \psi(\epsilon) &:= \text{ Strain energy density} \\
    \text{For example} &T(\epsilon) = K\epsilon^3 \\
    T(\epsilon) &= K\epsilon^3 \implies \psi(\epsilon) = K \frac {\epsilon^4}4 \\
    T(\epsilon) &= E\arctan(\epsilon) \\
    T(\epsilon) &= E \epsilon + K \epsilon^3 \\
                &\implies \psi(\epsilon) = E \frac {\epsilon^2} 2 + K \frac
                {\epsilon^4} 4
\end{align*}

Implies the momenntum equation 

\begin{align*}
    R_0 U{z z} &= R_0 F + \p_A \psi'(U_A) \\
    &= R_0 F + \psi'' (U_A) U_{A A} 
\end{align*}

Typically need that $\psi'' > 0$
so that the problem is well posed: $T(\epsilon) = \psi'(\epsilon) $ is monotone
in $\epsilon$. 

\subsubsection*{Energy balance} 

For a $A_\ell $ and $A_r$.  (note $\sigma$ is the cross sectional area).

Multiply momentum equatino by $\sigma U_t$.

\begin{align*}
    \sigma R_0 U_t U_{tt} &= \sigma R_0 F V + \sigma \p_A \psi'(U_A) U_z \\
\text{integrate over} A_\ell \text{ to } A_r \\
    \sigma R_0 \int_{A_\ell}{A_r} U_t U_{T t} \: dA &= \sigma R_0
    \int_{A_\ell}^{A_R} FV \: dA + \sigma \int_{A_\ell} ^{A_r} \p_A \psi'(U_A)
    U_z \: dA \\ 
                                                    &\text{Now integrate by parts}
    \\
    \sigma R_0 \int_{A_\ell}^{A_r} \frac d {dt} \frac {U_t^2} 2 \: da &= \sigma
    R_0 \int_{A_\ell} ^{A_r} FV \:dA + r[\psi'(U_A) U_t]_{A_\ell}^{A_R} - \sigma
    \int_{A_l}^{A_r} \psi'(U_A) U_{At} \: da  \\
    ... \\
    \implies
\end{align*}
$$\frac d{dt} \int_{A_\ell}^{A_r} \left[\sigma R_0 \frac {V^2}{2} + \sigma
\psi(U_A)\right] \: dA = \sigma \left(R_0 \int_{A_\ell}^{A_r} FV \: da +
TV|_{A_r} - TV |_{A_\ell}\right) $$

Where we can see

\begin{itemize}
    \item 
$\sigma R_0 \frac {V^2}{2}$ :  kinetic energy 
 \item 
      $ \sigma \psi(U_A)$ : strain energy
  \item 
$R_0 \int_{A_\ell}^{A_r} FV \: da $ : Power delivered by body forces
\item
$TV|_{A_r} - TV |_{A_\ell} $ : Power delivered by force acting on the tips

\end{itemize}

\section{Time-dependent elastic models}

For a linearly elastic material:

$$\begin{cases} R(A,t) = \frac {R_0} {1 + U_A} \quad(i)\\ 
    R_0 U_{t t} = R_0 F + EU_{AA}\quad(ii)\end{cases}$$

These are not coupled. You should solve $(ii)$ then $(i)$

From math2100: The wave equation on $\R$ 

$$(*) \begin{cases}
U_{t t} - c^2 U_{A A} = 0\\
U(t = 0, A) = f(A)\\
\p_t U(t = 0, A) = \rho(A) 
\end{cases}$$

You can get the solution using the d'Alembert solution to the wave equation.
This superimposes two waves, one that moves right and one that moves left at
the same speed. 

$$U(t,a) = \frac 12 (f(A - ct) + f(A + ct)) + \frac 1 {2c} \int_{A - ct}^{A +
ct} g(\hat A) \: d\hat A$$

Remark: This formula also gives a solution to $(*)$ on an interval with derichlet zero
boundary conditions ($U(t, A = 0) = 0 = U(t, A = \ell)$: when there is no

displacement at the endpoints of $(0,\ell)$). 

To get this initial condition you can extend $f,g$ as odd $2\ell$ periodic
functions. 

\begin{figure}[htpb]
    \centering
\incfig{exten}
    \caption{Extended function that is periodic on $2\ell$}
\end{figure}

\begin{ex}

    With a slinky: 

    $A = 0$ and $A = 10 = \ell_0$

    Intially move the loop of $A = 4$ to the left by $\Delta A = 2$, then
    release it. 

    \incfigw{exspr}{0.5}

    How do you write $U(t = 0, A) = U_0(A)$?
    
\end{ex}

\week{}
\lecture{}

We know 
\begin{align*}
    U_0(A = 0) &= 0 \\
    U_0(A = 10) &= 0  \\
    U_0(A = 4) &= -2 \\
\end{align*}

At $t = 0$ the slinky satisfies $EU_{A A} = 0 \implies U_A$ is a constant.   

\begin{align*}
    U_A &= \frac {-2 - 0}{4 - 0} = -\frac 1 2 \implies U = -\frac 12 A \\
    U_A &= \frac {0 - (-2)}{10 - 4} = \frac 13 \implies U = \frac 13 (A - 10) \\
    U_0 &=\begin{cases} -\frac 12 A & 0 \le A \le 4 \\ \frac 13 (A - 10) & 4 < A
    \le 10\end{cases}  \\
\end{align*}
The time-dependent problem is given by

\begin{align*}
    U_{t t} &= \frac E{R_0} U_{A A} \\
    U(0) &= 0 = U(10) \\
    U(t = 0, A) &=\begin{cases} -\frac 12 A & 0 \le A \le 4 \\ \frac 13 (A - 10)
    & 4 < A \le 10 \end{cases} \\
        \p_t U(t = 0, A) &= 0
\end{align*}

$\implies$ the solution is given by D'Alembert's formula where $\rho=0$ and $f$
 is the $2\ell$ periodic odd extension of $U_0$.  


 \begin{ex}
     Finite-length bar under external forcing (resonance)

     $U_{t t} = \frac E {R_0} U_{A A} + F(A,t) $ on the interval $(0,\ell)$.

     $U_A(t = 0, A = 0) = 0 = U_A(t = 0, A = \ell) $

     $F(A,t) = a(t) \cos(\kappa_n A) $ where $a(t) = \sin(\omega t)$

     $\omega$ is angular frequency of the forcing. 

     $\kappa_n = \frac {n \pi}\ell$ where $n$ is the wave number of the specific
     mode of the forcing term.

        \begin{center}
            \begin{tikzpicture}{notes-figure9}
                \begin{axis}[xtick={3.141},xticklabels={$\ell$}, ytick={},]
                    \addplot[domain=0:3.141 , samples=90, color=green]{ cos(deg(x))};
            \addlegendentry{$n = 1$}
                    \addplot[domain=0:3.141 , samples=90, color=black]{ cos(2 * deg(x))};
            \addlegendentry{$n = 2$}
                    \addplot[domain=0:3.141 , samples=90, color=blue]{ cos(3 * deg(x))};
            \addlegendentry{$n = 3$}
            \end{axis}
        \end{tikzpicture}
        \end{center}


    We can solve with Laplace transform to turn the pde into an ODE, transform
    the entire system of equations to get rid of the time derivatives


$$
F(s)=L(f(t))=\int_{0}^{\infty} e^{-s t} f(t) d t
$$

and the opposite 

$$
L(t f(t))=-\frac{d F(s)}{d s}
$$

\begin{align*}
    (*) \begin{cases} 
    s^2 \hat U - c^2 \hat U_{A A} &= \hat F (A,s) = \hat a(s) \cos(\kappa_n A) \\
                                  &\text{Boundary conditions:} \\
    \hat U_A(s, A=0) &= 0 = \hat U_A(s, A = \ell) 
\end{cases}
\end{align*}

To solve $(*)$ first consider the homogeneous equation: $\hat U_{hom} - \frac
{c^2}{s^2} \hat U_{A A}  = 0$

$$\hat U_{hom} = \alpha \exp(\frac cs A) + \beta \exp(-\frac cs A)$$

Find a particular solution : 

$$\hat U _{part} = \gamma \sin(\kappa_n A) + \delta \cos(\kappa_n A)$$


\begin{align*}
    \implies \p_{A A} U_{part} &= -\kappa_n^2 (\gamma \sin(\kappa_n A) + \delta
\cos(\kappa_n A)) \\
    s^2 \hat U_{part} - c^2 \p_{A A} \hat U_{part} &= (s^2 + c^2
    \kappa_n^2)(\sin(\kappa_n a) + \delta \cos(\kappa_n A)) \\
    &= \hat a(s) \cos(\kappa_n A) \\
    \implies \gamma = 0 \text{ and} \delta = \frac {\hat a(s)}{ s^2 + c^2
    \kappa_n^2}
\end{align*}

$\implies $ the general solution of $(*)$ is 

$$
\hat U = \frac {\hat a(s)}{ s^2 + c^2
    \kappa_n^2} \cos(\kappa_n A) + \alpha \exp(\frac cs A) + \beta \exp(-\frac cs A)
$$

Boundary conditions:

\begin{align*}
    \hat U_A &= -\kappa_n \frac {\hat a}{s^2 + \kappa_n^2 c^2} \sin(\kappa_n A )
    + \frac cs (\alpha \exp(\frac cs A) - \beta \exp(-\frac cs A)) \\
    U = \hat U_A(A = 0) &= \frac cs (\alpha - \beta) \implies \alpha = \beta  \\
    U = \hat U_A(A = \ell) &= \frac cs (\alpha \exp(\frac cs \ell) - \beta
    \exp(-\frac cs \ell))\alpha  \\
                           & \text{since }\ell \ne 0, \; \alpha \exp(\frac cs \ell) - \beta
    \exp(-\frac cs \ell)  
 \implies   \alpha = 0 \\
                           &\implies \hat U = \frac {\hat a} {s^2 \kappa_n^2
                           c^2} \cos (\kappa_n A)
\end{align*}

 \lecture{}

Now we compute the inverse laplace transform using the convolution theorem.

\begin{align*}
      \hat U &= \frac {\hat a} {s^2 \kappa_n^2 c^2} \cos (\kappa_n A) \\
             &\text{We have to use the convolution theorem} \\
      L^{-1}(\hat U) &= \cos(\kappa_n A) L^{-1}(\frac 1 {s^2 + \kappa_n^2
      c^2}) \times u(t)  
      \quad\quad  \text{(note $\cos(\kappa_n A)$ does not depend on $s$)} \\
     L^{-1}(\frac 1 {s^2 + \kappa_n^2
      c^2}) \times u(t) &= \frac 1 {k_n c} \int _0^t \sin(\omega(t - \hat t))
      \sin(\kappa_n c \hat t) \: :d\hat t \\
                        &= \frac 1 {\kappa_n c} \begin{cases} \frac 1 {\omega^2
                        \kappa_n^2 c^2}(\omega \sin(\kappa_n ct) - \kappa_n c
                        \sin(\omega t)) & \omega \ne \kappa_n c\\
                    -\frac 12 
                    \mathcolorbox{yellow}{t} 
                    \cos(\kappa_n ct) + \frac 1 {2\kappa_n c}
                    \sin(\kappa_n ct)   & \omega = \kappa_n c
                \end{cases} \\
\end{align*}

$\implies \omega = \kappa_n c$ is the resonance cases; the amplitude of $U$
grows linearly in time (due to the highlighted \mathcolorbox{yellow}{$t$}). The
reason is $\kappa_n c$ are the eigenmodes of the structure. 


\marginnote{Remember $R_0$ is density, $E$ is elastic modulus.}

Since the eigenfrequenceies $\kappa_n c$ are proportional to $c = \sqrt{\frac E
{R_0}}$ one can probe the material propertiy $\sqrt{\frac E {R_0}}$ through an
acoustic signal (frequancies of resonances in the material will be $\kappa_n
\sqrt{\frac E {R_0}}$). 

 \end{ex}

 \subsection{Viscoelastic Materials}

 \begin{itemize}
     \item Introducing damping to the stress-strain relation. 
 \end{itemize}

 \begin{ex}

 For example, considering the slinky again. It is unrealistic for it to
 oscillate forever, instead some sort of damping, eg air-external friction,
 internal friction (resistance to deformation) will slow it down.

 If you have a spring hanging freely with a weight and it is in equilibrium, you
 have the spring force $F_s = -K_s u = mg$, and a damping force $F_d-c\dot u = 0$. 
 When it moves the damping force creates a resistance. Resitive force is proportional 
 to the velocity: $F_d = -c \dot u$. 

 Newtond's Law: $m \ddot u = F_d + F_s$  gives the damped oscillation equation

 $$m \ddot u + c\dot u +  K u = 0$$

 These systems provide a framework to model viscoelasticity.
 \end{ex}

 \subsubsection{Kelvin-Voigt element}

 You have a (horizontally) a spring and dashpot in parallel. 

 $$F = - (F_s + F_d)$$

 So displacements are equal: $u_s = u_d = u$.   

 Spring
 constant $K$.

 $$F_s = -K u$$
 $$F_d = -c \dot u$$

 To create a constitutive law for stress, identify $F$ by stress $T$ and $u$ b
 strain $\epsilon$. Also replace $\dot F$ by $\partial_t T$ (stress rate), and
 $\dot u$ by $\p_t \epsilon$ (strain rate).

 $$T = E(\epsilon + \tau_1 \p_t \epsilon ) $$

 $E$ elastic modulus

 $\tau_1$ dissipation time of the strain.

 \subsubsection{Maxwell element}

 Spring and dashpot in series. 

Displacements sum to $u$:  $u = u_s + u_d$

Forces are equal: $F = - F_s = -F_d$

Constitutive law: 

$F = -F_s= ku_s \implies u_s = \frac FK \implies \dot u_s =
\frac {\dot F }K $

$F = -F_d = c \dot u_d \implies \dot u_d = \frac FC$ 

So 
$\dot u = \dot u _s + \dot u_d = \frac {\dot F}k + \frac F c$

$(* * ) \quad\quad c \dot u = \frac ck \dot F + F$

Which gives the constituive law

$$T + \tau_0 \p_t T = E\tau_1 \p_t \epsilon$$


\lecture{}

\subsubsection{Standard Linear Model}

\incfigw{dsprlin}{0.7}

To model the spring and despot in series use the resulf for the Maxwell element
$(* * ) $.

$$c_1 u' = \frac {c_1}{k_1} F_1' + F_1 = \frac {c_1}{k_1} (F' - k_0 u') + F -
k_0 u$$

$$\implies k_0 u + u(c_1 + \frac {k_0 }{k_1} c_1) = F + \frac {c_1}{k_1} F'$$

Rewrite this as

$$F + a_1 F' = a_2 u + a_3 u'$$

Where $a_1 = \frac {c_1} {k_1}$,  $a_2 = k_0$ and $a_3 = c_1 + \frac
{k_0}{k_1}c_1$ so that $a_3 = c_1 + \frac {k_0}{k_1} c_1 > \frac {k_0}{k_1} c_1
= a_1 a_2$

This gives the constitutive law

$$T + \tau_0 T_t = E(\epsilon + \tau_1 \p_t \epsilon)$$

\begin{align*}
    \tau_0 &= a_1 = \frac {c_1}{k_1} \\
    E &= a_2 = k_0 \\
    E T_1 &= a_3 
\end{align*}

And as we showed 

$ET_1 = a_3 > a_1 a_2 = ET_0$

$E = a_2 = k_0$  

Which $\implies \tau_1 > \tau_0$, where $\tau_0$ and $\tau_1$ are dissipation
times, $E$ the elastic modulus.


$$\begin{cases}
    R = \frac {R_0}{1 + U_A} \\
    R_0 U_{t t} = R_0 F + T_A \\
    T + \tau_0 T_t = E(U_A + \tau_1 U_{At})
\end{cases}$$

\subsubsection{Summary}

\emph{Standard linear model}
$$T + \tau_0 \p_t T = E(U_A + \tau_1 \p_t \epsilon)$$

\emph{Kelvin-Voigt model}
$$T = E(\epsilon + \tau_1 \p_t \epsilon ) $$

Note the above two models involve the strain (in the right hand side), so they will eventually go back to
the reference configuration. These are viscoelastic solids. 

\emph{Maxwell Model}
$$T + \tau_0 \p_t T = E\tau_1 \p_t \epsilon$$

This is a visco-elastic fluid since it only includes the strain rate, it will
not return to the reference configuration. It is elastic on fast time scales,
but on a long timescale you have a creeping motion that deforms the fluid. 

\subsubsection{Relaxation Function}

This unifies the above three models, by adding constant $C$ to the standard linear model. 

$$T + \tau_0 \p_t T = E(C\epsilon \tau_1 \p_t \epsilon)$$

\begin{itemize}
    \item[] $\tau_0 =0, C = 1$ Kelvin-Voigt model
    \item[] $C = 0$ Maxwell model
    \item[] $C = 1$ Standard linear model
\end{itemize}

Solving 
$$\begin{cases} T + \tau_0 \p_t T = h \\
T(t=0) = 0
\end{cases}$$ 

where $h = E(C\epsilon + \tau_1 \p_t \epsilon)$, (and assuming $\epsilon(0) = 0$)


$$T = \int_0^t E(C + (\frac {\tau_1}{\tau_0} - c)e^{\frac {-(t -
\hat t)}{\tau_0}})\p_t \epsilon{\hat t}\: d\hat t$$

Where the relaxation function 

$$G(t - \hat t) = E(C + (\frac {\tau_1}{\tau_0} - c)e^{\frac {-(t -
\hat t)}{\tau_0}})$$

\begin{itemize}
    \item Standard linear model
        $$G(t) = E(1 + (\frac {\tau_1}{\tau_0} - 1) e^{\frac {-t}{\tau_0}})$$
    \item Maxwell element

        $$G(t) = E \frac {\tau_1}{\tau_0} e^{\frac {-t}{\tau_0}}$$

    \item Kelvin-Voigt

        Note cannot simpoly sat $\tau_0$ to zero but it is attainable using a
        different method.

        $$T = E(\epsilon + \tau_1 \p_t \epsilon) = \int_0^1 E(\p_t \epsilon(\hat
        t) + \tau_1 \rho_0 (t - \hat t) \p _t \epsilon(\hat t) ) \: d \hat t
        \implies$$

        $$G(t) = E(1 + \tau_1 \delta_0(t))$$
\end{itemize}

A general formulation for constitutive relation for stress ($T$) is

$$T = \int_0^t G(t - \hat t) \epsilon_t (\hat t) \: d\hat t = G * \p_t \epsilon $$

Here $G \ge 0$ can be choses quite arbitrarily (want it to be integrable, decay
to zero).

This is the ``visco-elastic law of relaxation type''.

\week{}
\lecture{}

Instead you may also write the stress as an integral against the strain and
integrate by parts in $T_2$ using $\epsilon (0) = 0$:

$$T_2 = \int_0^t E \frac {T_1}{T_0} E^{-\frac {t - \hat t}{\tau_0}} \p_t
\epsilon(\hat t) \: dt = 
 E\frac {T_1}{T_0} \epsilon(t) - \int_0^t E\frac {T_1}{T_2} e^{-t - \hat t}
{\tau_0} \epsilon(\hat t) \: d \hat t$$

For example for the standard linear model ($C = 1$) this implies 

$$T = T_1 + T_2 = \underbrace{E \frac {T_1}{T_0}
\epsilon(t)}_{\begin{minipage}[adjusting]{5em}instantaneous elastic stress \end{minipage}}
- \underbrace{\frac E {\tau_0} \int_0^t \left( \frac {\tau_1}{\tau_0} - 1\right ) e^{\frac
{t  - \hat t}{\tau_0}}  \epsilon(\hat t) \: d \hat
t}_{\begin{minipage}{5em}damping\end{minipage}}
$$

{The computations clearly get complicated so often you can't solve
viscoelastic models explicitly very easiliy}

\begin{ex}
    Maxwell-viscoelastic bar occupying $0 \le A < \infty$ with body force $F =
    0$.

    Momentum equation: 

    $$(*) \qquad\quad R_{0} U_{t t} - T_A = \p_A (G \star \epsilon_t) = \int_0^t G(t - \hat t)
    U_{AA_t}$$

    Where $G(t) = E \frac {\tau_1}{\tau_0} e^{\frac {-t} {\tau_0}}$.

    We have the ``far field condition'' on the right, the displacement vanishes at
    infinity: 

    $$\lim_{A \to \infty} U(A,t) = 0$$

    $$U(t=0, A) = 0 = \p_t U(t = 0, A)$$

    How to enforce $U_A(t , A=0) = F(t) $ experimentally? 

    Apply a force to the left endpoint: $T(A = 0, t) = \int_0^t G(t - \hat t) F_t(\hat t) \: d\hat t $

    Use laplace transform on $(*)$, $\hat U = L(U)$
    
    $$\begin{cases}
    R_0 s^2 \hat U = L(G) s \hat U_{A A} $ where $L(G) = \frac {\hat E}{s +
    \frac 1 {\tau_0}}  \\
    \hat U_A (A = 0, s) = \hat F(s) \\
    \lim_{A \to \infty} \hat U(A, s) = 0
    \end{cases}$$


    \begin{align*}
        R_0 s \hat U = \frac {\hat E}{s + \frac 1 {\tau_0} } \hat U_{A A}
        &\implies \frac {R_0 s}{\hat E}(s + \frac 1 {\tau_0}) \hat U = \hat U
        _{A A}  \\
        &\implies \hat U = \alpha (s) e ^{\omega (s) A} + \beta(s) e^{-\omega(s) A}
        \qquad {\text{ $($ODE for }\hat U)}\\
        \text{where } \omega(s) &=\sqrt{\frac {R_0} s {\hat E} (s + \frac 1
        {\tau_0})}  \\
    \end{align*}

    \marginnote{laplace transform only works when the real part of $s$ is large
    enough compared to the imaginary part.}

    Note that $Re(\omega) > 0$ if $Re(s) $ is large enough. 
    Therefore $\alpha(s) = $ to guarantee that $\lim_{A \to \infty} \hat U(A) =
    0$.

    Find the other constant: 

    \begin{align*}
        \implies & \hat U(A,s) = \beta(s) e^{-\omega(s) A} \\
        \implies & \hat U_A(A,s) = -\omega(s) \beta(s) e^{-\omega(s) A} \\
        {\text{2nd boundary cond. }} \hat F(s) &= \hat U_A(A = 0, s)  \\
                                               &=-\omega9s) \beta(s)   \\
        \implies  \beta(s) &= \frac {-\hat F(s)}{\omega(s)} \\
        \implies  \hat U(A, s) &=  -\frac 1 {-\omega (s)} \hat F(s) e^{-\omega(s) A}  \\
    \end{align*}

    Now using the convolution theorem to inverse laplace transform:

    \begin{align*}
        U(A,t) &= F(t) \star \underbrace{L^{-1}\left(-\frac 1 {-\omega (s)} e^{-\omega(s) A}
        \right)}_{(* *)} \\
                (* *) &=  \sqrt{\frac {\hat E}{R_0}}
        e^{\frac {-t}{\tau_0}}I_0 \left(\frac {-1}{2\tau_0} \sqrt{t^2 - \frac {R_0}{\hat E}
        A^2}\right) \times H\left(t - \sqrt{\frac {R_0}{\hat E}} A\right)\\
                I_0 & {\text{ is the Bessel function of the 1st kind }} \\
                H & \text{ is the heaviside function } \\
                U &= \int_0^t Q(A,\hat t) H(\hat t - \sqrt{\frac {R_0}{\hat E}} A) F (t - \hat t)
                \hat H \\
                  &=\int_{\sqrt{\frac {R_0}{\hat E}}A}^t Q(A, \hat t) F(t - \hat t) \: d\hat t \times H\left(t - \sqrt
                  {\frac{R_0}{\hat E}}A\right)  \\
                \implies & U \text{ is zero wherever } t < \sqrt{\frac {R_0}A} , \\ 
                \implies & A > t\sqrt {\frac {\hat E} {R_0}} = t \sqrt{\frac {E \tau_1}{\tau_0}
                \frac 1 {R_0}} \\
                    \text{For } F(t) &= \sin(t) \\
    \end{align*}

    The speed of propagation in the linearly elastic material $c = \sqrt{\frac E {R_0}}$ is smaller
    than the propagation in the viscoelastic model $\sqrt{\frac E{R_0} \frac {\tau_1}{\tau_0} }$ 
    $(\tau_1 > \tau_0)$. And for $\tau_0 \to 0$ (no memory, purely viscous) the speed of propagation
    become infinite just like in the diffusion equation. 
\end{ex}

\lecture{}

\subsubsection{Dynamic (complex) modulus}

\begin{itemize}
    \item Dynamic mechanical analysis
    \item ie ``The {Vibration test}''
\end{itemize}

A bar of infinite length connected by a rod to a wheel that forces the bar periodically. We are
aiming to determine the properties of the bar, so do not assume anything about its elastic
properties. 


$U(A = 0, t) = a \sim(\omega t)$  apply a sinusoidal displacement at the tip. 

$\lim_{A \to \infty} U(A,t) = 0$.

$T + \tau_0 \p_t T = E(\epsilon + \tau_1 n\p _t \epsilon) = E(U_a + \tau_1 U_{At}) \qquad (*)$

$R_0 U_{t t} = T_a \qquad(* *)$

it is easier to write everything in complex coordinates: 

$U(A = 0, t) = e^{i \omega t}$, $(\sin(\omega t))$ is the complex part 

$U(A,t) = \bar U(A) e^{i \omega t}$

$T(A,t) = \bar T (A) e^{i \omega t}$

Substitute in $(*)$

\begin{align*}
    \bar T(A) e^{i \omega t} (1 + i\omega \tau_0)) &= E \bar U_A (A) e^{i \omega t} (1 + \tau_1 i
    \omega) \\
    \bar T(A) &= E \frac {1 + i \omega \tau_1}{1 + i \omega \tau_0} \bar U_A(A) = E^* \bar U_A(A) 
\end{align*}

Which is the complex stress-strain relation under oscillatory conditions (because it looks very
similar to the linear law for it but in terms of complex numbers). 

\begin{align*}
    E^* &= E \frac {1 + i\omega \tau_1}{1 + i \omega \tau_0} = E \frac {(1 + i \omega \tau_1)(1 - i
    \omega \tau_0)}{(1 + i \omega \tau_0)(1 - i \omega \tau_0)} \\
        &= E \frac {1 + \omega^2 \tau_0 \tau_1}{1 + \omega^2 {\tau_0} ^2} + i E \frac {\omega (\tau_1
    - \tau_0)}{1 + \omega^2 {\tau_0}^2} \\
    &= E' + i E'' \\
        E': & \text{ storage modulus (`the elastic part') } \\
        E'': & \text{ loss modulus (`the viscous part') } \\
             & \text{from the momentum equation } (* * ) : R_0 U_{t t} = T_A \\
        R_0 \bar U(A) (-\omega^2) e^{i \omega t} &= \bar T_A(A) e^{i \omega t} 
                = E* \bar U_{A A} (A) e^{i \omega t} \\
        \frac {-\omega^2 R_0}{E^*} \bar U(A) &= \bar U_{A A} (A) \\
        \text{ make } \kappa &= \sqrt{\frac {\omega^2 R_0} {E^*}} \\
                             &\text{The general solution is: } 
        \bar U(A) = \alpha e^{i \kappa A} + \beta e^{i \kappa A} \\
        \kappa^2 &= \frac {\omega^2 R_0} E \left(\frac {1 + \omega^2 \tau_0 \tau_1}{1 + \omega^2
        \tau_1^2} + i \frac {\omega(\tau_0 - \tau_1)} {1 + \omega^2 \tau_1^2}\right)  \\
                Im(\kappa) < 0 &\implies Re (i\kappa) > 0\\
                               & \implies \alpha = 0 \text{ to satisfy }  \lim_{A \to \infty} \bar U(A) = 0 \\
                U(A = 0, t) &= \alpha e^{i \omega t} \\
                            & \implies  \bar U(0) e^{i \omega t} = \beta e^{-i \kappa 0} e^{i \omega t} 
                                =  \beta e^{i \omega t} \\
                            & \implies  a = \beta \\
                            &\implies \bar U(A) = a e^{i \kappa A} 
\end{align*}

Now want the stress at the endpoint $A = 0$

\begin{align*}
    T(A = 0, t) &= \bar T(A = 0) e^{i \omega t} = E^* \bar U_A(A = 0) e^{i \omega t} \\
                &=  E^*(-i) \kappa a e^{i \omega t}\\
\end{align*}

\begin{itemize}
    \item You increase the freqency forcing and for every frequency you measure the phase
        difference. 
\end{itemize}


\lecture{}

What is the displacement $U$ and stress $T$ at the left end point.

$U(A = 0, t) = \bar U(A = 0) e^{i \omega t} a e^{i \kappa 0}e^{i \omega t} = a e^{i \omega _t}$

$T(A = 0, t) = \bar T(A = 0) e^{i \omega t} = E^* (-i) \kappa a e^{i \omega t}$

$T$ In polar coordinates: 
\begin{align*}
    -i &= e^{-\frac \pi 2}, \\
    \kappa a e^{i \omega t} &= \kappa a e^{i \omega t} \quad\text{and, } \\
    E^* &= re^{i \delta} = E^* \sqrt{R_0 \frac {\omega^2}{E^*}} = \sqrt{R_0 \omega^2
E^*} \\ 
    r^2 e^{2 i \delta}  &= R_0 \omega^2 E^* = R_0 \omega^2(E' + iE'') \\
                        & \text{take the argument } 2\delta = \arctan \frac {E''}{E'} \\
    \implies \delta &= \frac 12 \arctan {\frac {\omega(\tau_1 - \tau_0)}{1 +
    \omega^2 \tau_0 \tau_1}} \\
            \implies T(A = 0, t) &= a r \exp(i (\omega t - \frac \pi 2 + \delta)) 
\end{align*}

Phase difference: $\delta - \frac \pi 2$

\subsubsection*{Interpretation}

\begin{enumerate}
    \item 
If you have an elastic material, then $\tau_1 = \tau_0 = 0$, so $\delta = 0$ and
the phase difference is $-\frac \pi 2$ independently of the frequency $\omega$
\item
if you have $\tau_0 = 0$ (the Kelvin-Voigt model of a viscoelastic solid) $\delta = \frac 12
\arctan(\omega \tau_1)$, which is a monotone function in $\omega$
\item 
    In general $\delta$ has a maximum of $\omega = \frac 1 {\sqrt{\tau_0
    \tau_1}}$
    \begin{align*}
        g &= \frac {\omega (\tau_1 - \tau_0)}{1 + \omega^2 \tau_0 \tau_1} \\
        g' &= \frac {\sqrt{\tau_1 - \tau_0}}{1 + \omega_2 \tau_0 \tau_1} - \frac
        {\omega(\tau_1 - \tau_0)}{(1 + \omega^2 \tau_0 \tau_1)^2 } 2 \omega
        \tau_0 \tau_1 = 0\\
        0 &= 1 + \omega^2 \tau_0 \tau_1 - 2\omega^2 \tau_0 \tau_1 \\
        \omega &= \frac 1 {\sqrt {\tau_0 \tau_1}}
    \end{align*}
\end{enumerate}


\pagebreak
\chapter{Continuum Mechanics in 3D}

\section{Navier-Stokes Derivation}

\begin{description}
    \item[Material Coordinates]  (lagrangian)

    \item $A \in \R^3$ the material coordinate/position at $t = 0$
    \item $X(A,t) \in \R^3$, $X(A,0) = A$
    \item $U(A,t) = X(A,t) - A \in \R^3$ displacement
    \item $V(A,t) = \p_t U(A,t) $ velocity
    \item[Spatial/Eulerian coordinates]  $x\in \R^3$
    \item $u(x,t) = U(a(x,t) , t) \Leftrightarrow u(X(A,t) , t) = U(A,t)$
    \item $v(x,t) = V(a(x,t) , t) \Leftrightarrow v(X(A,t) , t) = V(A,t)$
\end{description}

Usually we use material coordinates for solids and spatial coordinates for
fluids. 

\subsection{Deformation Gradient}

Consider $A_0 \in \R^3$ and $A = A_0 + \Delta A \in \R^3$

\marginnote{Use the taylor expansion only keeping the first derivative term}
$x = X(A_0 + \Delta A, t) = \underbrace{X(A_0, t)}_{X_0} + \underbrace{\nabla_A
X(A_0,t)}_{F(A_0, t)} \Delta A + o(\Delta A)$

\begin{definition}{}{} Deformation gradient
    $$\nabla _A X(A_0,t) = F(A_0, t) $$

$x - x_0 \sim F(A_0, t) (A - A_0)$

eg at $t = 0$ $X(A,t = 0) = A)$, and  $F = \nabla_A X = \nabla_A A = \mathbb I$

\end{definition}

\subsubsection{Impenetrability}

The assumption is you can always go from material coordinates to spatial
coordinates and back: so the $A$ must always have an inverse, so $\det F \ne 0$. 

$$J = \det F = \det \underbrace{\nabla_A X}_{\text{Jacobi matrix}} \ne 0$$

This is essentially the impenetrability assumption from 1D continuum mechanics in 
three dimensions.

At $t = 0$, $X(A, t= 0) = A \implies F = \mathbb I \implies J = 1$

Since $J(t = 0) = 1$, $J \ne 0$ always implies $J(t) > 0$ 

\week{}
\lecture{}


\begin{ex} Uniform dilatation

        Motion according to $A \mapsto X(A,t) = \alpha(t) A$ where
            $\alpha(t) \ge 0, \alpha(0) = 1$.

            Consider a sphere centred at $\hat 0 =  \begin{pmatrix} 0\\0\\0
            \end{pmatrix} $ with radius $r = 1$ at $t = 0$ and at $t = 1$,
            $\alpha(1) = 2$, (the circle is being stretched according to its
            radius linearly increasing).

            Material Coordinates

            $U(A,t) = X(A, t) - A = (\alpha(t) - 1)A $

            $V(A,t) = \p_t U = \alpha' (t) A$

            Spatial coordinates

            Invert $x$: $A = \frac 1 {\alpha (t)} x$

            $u(x,t) = U(\frac 1 {\alpha(t)}x, t) = \alpha(t) - 1 \frac 1\
            {\alpha(t)} x = (1 - \frac 1 {\alpha(t)}) x$

            $v(x,t) = v(\frac 1{\alpha(t)} x, t) = \alpha ' (t) \frac 1
            {\alpha(t)} x$

            $\implies \p_t u \ne v$ like in one dimension.

            $F = \nabla_A X = \nabla _A(\alpha(t) A) = \alpha(t) \nabla_A
            \begin{pmatrix} A_1\\A_2\\A_3 \end{pmatrix} = \alpha(t)
            \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix}$ 

            To guarantee $J > 0$, $J = \det F = \det (\alpha(t) \mathbb I) = \alpha(t)^3$

            $\implies$ we need $\alpha(t) > 0$.
\end{ex}

\begin{ex}
    Simple shear

    Motion according to:

    \marginnote{There is a shear along the $y$ axis}
    $$
    X(A,t) = \begin{pmatrix} A_1 + \alpha(t) A_2 \\ A_2 \\ A_3 \end{pmatrix}
    \quad\text{Where } \alpha(0) = 0 \text{ so } X(A,0) = A
    $$

    \incfig{shear}
    Material coordinate:

    $U(A,t) = \begin{pmatrix} \alpha(t) A_2 \\ 0 \\ 0 \end{pmatrix} 
    \qquad \qquad
    V(A,t) = \begin{pmatrix} \alpha'(t) A_2 \\ 0 \\0 \end{pmatrix}$

    Spatial coordinate:

    Invert $X = \begin{pmatrix} A_1 + \alpha(t) A_2 \\ A_2 \\ A_3 \end{pmatrix}
    = \begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix} $.

    $A_1 = x_1 - \alpha(t) A_2 = x_1 - \alpha(t) x_2$

     $A_2 = x_2$ and $A_3 = x_3$ So

     $$A = \begin{pmatrix} x_1 - \alpha(t) x_2 \\ x_2 \\ x_3 \end{pmatrix} $$

     So

     $$u(x,t) = \begin{pmatrix} \alpha(t) x_2 \\ 0 \\0 \end{pmatrix} $$

 $$v(x,t) = V\left( \begin{pmatrix} x_1 - \alpha(t) x_2 \\ x_2 \\ x_3  \end{pmatrix}
 , t\right) = \begin{pmatrix} \alpha ' (t)x_2 \\ 0\\ 0  \end{pmatrix} $$

 Deformation gradient

 $F = \nabla_A  X = \nabla _a \begin{pmatrix} A_1 + \alpha(t) A_2 \\ A-2 \\ A_3
     \end{pmatrix} = \begin{pmatrix} 1&\alpha(t)&0 \\ 0&1&0 \\ 0&0&1
 \end{pmatrix} $

 $\det F = 1 > 0$ which does not imply any additional constraints. 


\end{ex}


\subsection{Material Derivative}

Let $F = (F_1,F_2,F_3)$ respecitvely to $f$ be a quantity in material respective to spatial
coordinates then we ahve

\marginnote{$X = \begin{pmatrix} X_1\\X_2\\X_3 \end{pmatrix} $}

$$F(A,t) = f(X(A,t), t) $$

\begin{align*}
    \frac {\p F }{\p t} &= \frac d {dt} f (X(A,t) , t) \\ 
                        &= \p_t f + \p _{X_1} f
\frac {\p{X_1}}{\p t} + \p _{X_2} f \frac {\p X_2} {\p t} + \p _{X_3} f \frac
    {\p X_3}{\p t} \\
                        &= \p _t f + \nabla f \cdot \p _ t X \\
&= \p_t f + \nabla f \p _t U   \\
&= p_t f + \nabla f \cdot v \\
&= \p_t f + (v \cdot \nabla) f = (\p_t + v \cdot \nabla) f \\
\end{align*}

{$(v \cdot \nabla)$ is notation for a differential operator and does
not make sense on its own, rather it is notation for the previous statement}

\begin{align*}
    (v \cdot \nabla) f = \nabla f \cdot v &= 
    \begin{pmatrix} 
                f'_{X_1} & f'_{X_2} & f'_{X_3} \\
                f''_{X_1} & f''_{X_2} & f''_{X_3} \\
                f'''_{X_1} & f'''_{X_2} & f'''_{X_3} \\
            \end{pmatrix} \cdot \begin{pmatrix} v_1 \\ v_2 \\ v_3 \end{pmatrix}
            \\
            &=      \begin{pmatrix} 
                f'_{X_1} v_1 + f'_{X_2} v_2 + f'_{X_3} v_3 \\
                f''_{X_1}v_1 + f''_{X_2}v_2 + f''_{X_3}v_3 \\
                f'''_{X_1}v_1 + f'''_{X_2}v_2 + f'''_{X_3}v_3 \\
            \end{pmatrix} \\
                           &= \begin{pmatrix} \nabla f^1 \cdot v \\ \nabla f^2 \cdot v \\ \nabla
            f^3 \cdot v\end{pmatrix} 
\end{align*}

$\implies$ the material derivative: 

$$D_t f = \p_t f + \nabla f \cdot v$$

The rate of change change from the perspective of fixed world coordinates, 
so how both the thing moves as a result of its own movement and it being 
transported by the material.

\begin{ex}
    Uniform dilatation

    $u(x,t) = 1 - \frac 1 {\alpha(t)}x$

    $v(x,t) = \frac {\alpha'} {\alpha} x$


    \begin{align*}
        D_t u &= \p _t u + \nabla u \cdot v = \frac {\alpha^2} {\alpha^2} x  + (1 -
    \frac {\alpha (t)}) \mathbb I \cdot \frac {\alpha ' }{\alpha } x  \\
              &= \frac {\alpha^2} {\alpha^2}x + (1 - \frac 1 {\alpha(t)}) \cdot
              \frac {\alpha '}{\alpha} x \\
              &= \frac {\alpha'}{\alpha} x \\
    \end{align*}

    Which is just the velocity field.

\end{ex}

\lecture{}

\subsection{Reynold's Transport theorem in 3D}
\marginnote{October 8 Lec 29}

The major technical tool to derive it is the fact that for $M(t) \in \R^{3\times
3}$ (smooth, invertible) it holds that. 

$$\frac d {dt} \det (M(t)) tr(M^{-1} \frac d{dt} M)$$

For $M = F = \nabla_A X$ and $j = \det F$.

$$\frac d {dt} J = J tr (F^{-1} \frac d {dt} F) = \left|\frac d {dt} \nabla_A (X - A)
= \nabla_A V \right|$$

{In the trace of the matrix, the multiplication commutes $|tr(AB) =
tr(BA) | \implies J tr(F^{-1} \nabla_A V) = J tr(\nabla_A V F^{-1})$}


Chain rule: $\nabla_A V = \nabla_Av(X(A,t),t) = \nabla_{\vec x} \vec v \cdot \nabla_A X = \nabla_{\vec x} \vec v
\cdot F$


    $$J tr(\nabla_A V F^{-1}) = J tr(\nabla_{\vec x } \vec v) = \begin{pmatrix} 
        v_{x_1}^1 & v_{x_2}^1 &v_{x_3}^1 \\ 
        v_{x_1}^2 & v_{x_2}^2 &v_{x_3}^2 \\ 
        v_{x_1}^3 & v_{x_2}^3 &v_{x_3}^3 \\ 
\end{pmatrix} = J \nabla \cdot \vec v $$

Now consider $R(t = 0) \suybset \R^3$, a piece of volume with a piecewise smooth
surface and $R(t) = X(R(t = 0), t)$.

\marginnote{Want to create a framework to use the balance laws in 3D.}
Consider a quantity $f(x,t)$ (smooth).

\marginnote {Change of variables to use material coordinates}
$$\frac d{dt} \iiint _{R(t)} f(x,t) \; d\vec x = \frac d {dt} \iiint_{R(t=0)} f(X(A,t) ,
t) \det (F)\; d A_1\: d A_2 \: d A_3 $$

$$= \iiint_{R(0)} \underbrace{\frac d {dt} f(X(A,t),t)}_{D_t f} J + f(X(A,t),t) \p_t J\; d\vec A $$

$$= \iiint_{R(0)} (D_t f + f(\nabla \cdot \vec v)) J \; d\vec A = \iiint_{R(t)}
(D_t f + f (\nabla \cdot \vec v)) \; d \vec x$$

$$= \iiint_{R(t)} (\p_t f + \underbrace{\nabla \cdot (f \vec v))}_{=\nabla f
\vec v + f \nabla \cdot \vec v} \; d\vec x = \iiint_{R(t)} \p_t f \; d\vec x +
\iint_{\p R(t) f \vec v \cdot \vec nn \: dS(\vec x)}$$

Which uses the divergence theorem:

$$\iiint_R \nabla \cdot \vec g \: d \vec x = \iint_{\p R} \vec g \cdot \vec n dS(\vec
x)$$


We get two ways of writing \textbf{Reynold's Transport Theorem} in 3D:

\begin{align*}
    \frac {d} {dt} \iiint_{R(t)} f (\vec x ,t )\; d\vec x &= \iiint _{R(t)} \p_t f
d \vec x + \iint_{\p R(t)} f \vec v \cdot \vec n \; dS(\vec x) \\
                                                          &= \iiint_{R(t)} (D_t
                                                          f + f(\nabla \cdot
                                                          \vec v)) \; d\vec x
\end{align*}

\subsection{Balance Laws in 3D}

Note; $\vec J :=$ flux, $J :=$ jacobi determinant.

$$\frac d {dt} \iint_{R(t)} f(\vec x, t ) \; d\vec x = - \iint_{\p R(t)} \vec J
\cdot \vec n \; d S(\vec x) + \iiint_{R(t)} Q(x,t) \; dV$$

$Q(\vec x, t) :=$ creation or depletion. 

Using reynold's transport theorem and the divergence theorem: 

$$\iiint_{R(t)} (D_t f + f(\nabla \cdot \vec v) ) \; d \vec x = \iiint_{R(t)}
(-\nabla \cdot \vec J + Q(\vec x, t)) \; d\vec x$$

Since $R(t) = X(R(0), t) $ is arbitrary, we use the dubois Reymond lemma to
conclude that the integrand vanishes: 

$$D_t f + f(\nabla \cdot \vec v) = - \nabla \vec J + Q(\vec x, t)$$

$$\p_t f + \nabla \cdot (f \vec v) = - \nabla \cdot \vec J + Q(\vec x, t)$$

Which are the two differential versions of the balance law.

\subsubsection{Continuity Equation}

Consider the mass density $\varrho = \varrho(\vec x, t)$, (mass per volume). 

Mass conservation:  $\vec J = 0, Q = 0$.

$$\frac d {dt} \iiint_{R(t)} \varrho(\vec x ,t) \; d\vec x = 0$$ 

The corresponding differential formulation is:

$$D_t \varrho + \varrho (\nabla \cdot \vec v) = 0$$

$$\p_t \varrho \nabla \cdot (\varrho \vec v) = 0$$

If the material is incompressible $(\varrho = $ a constant $)$, it holds that
$\p_t \varrho = \p_{x_i} \varrho = 0$ and the continuity equation
$\p_t \varrho + \nabla_x \varrho \vec v + \varrho(\nabla \cdot \vec v) = 0$
reduces to

$$\nabla \cdot \vec v = 0$$

\begin{ex} $ $

    \begin{itemize}
        \item 
    Shear flow $X(A,t) = \begin{pmatrix} A_1 + \alpha(t) A_2 \\ A_2\\A_3
        \end{pmatrix}$. In this case we had $\vec v = \begin{pmatrix} \alpha'(t)
    X_2 \\ 0 \\ 0 \end{pmatrix}   \implies  \nabla  \cdot \vec v = \p_{x_1}
    (\alpha'(t) X_2) = 0 \implies$ simple sheer is possbile for an
    incompressible material. 
        \item 
            Uniform dilatation: $X(A,t) = \alpha(t) A$
            We found $\vec v (\vec x, t) = \frac {\alpha '}{\alpha} \vec x$. 

            $$\nabla \cdot v = \nabla \cdot (\frac {\alpha'}{\alpha} \vec x) =
            \frac {\alpha'}{\alpha}(\p_{x_1}x_1) + \p{x_2} x_2 + \p_{x_3} x_3 = 3
            \frac {\alpha ' } {\alpha}$$

            So uniform dilatation is not possible with an incompressible
            material.
    \end{itemize}
\end{ex}

\week{}
\lecture{}

\subsection{Conservation of linear (translational momentum)}

Balance law for momentum 

$$\iiint _{R(t) } \varrho (\vec x, t) \vec v (\vec x, t) \: d\vec V = \iint \vec
t \: ds + \iint \varrho(\vec x, t) \vec f (\vec x , t) \: d\vec V $$

$\vec f \in \R ^3$ is body forces

$\vec t \in \R^3 $ is stress $force / area$ through the boundary at $\vec x$.
(through an infinitesmilly small segment of the surface: the orientatino is the
important part.) $\vec t = \vec t(t, \vec x, \vec n)$. ($t$ means time $\vec t$
stress)

Consider the normal stresses, 

\begin{itemize}
    \item 
$\vec t_x$ stress through out the plane $x=0$
\item $\vec t_y $ stress throughout the plane $y=0$ with normal $\vec n_y = \begin{pmatrix}
    0\\1\\0 \end{pmatrix} $
\item $\vec t_z$ stress through the plane with the normal $\vec n_z = \begin{pmatrix} 0 \\
    0 \\ 1\end{pmatrix} $
\end{itemize}

To relate the $\vec t$ to $\vec t_x$, $\vec t_y$ and $\vec t_z$ consider the
cauchy tetrahedron. 

\subsubsection{Cauchy Tetrahedron}

A tetrahedron with side lengths $A,B,C$, where each side is one of these
fundamental planes, and the edges are $\vec A = (A,0,0),\; \vec B = (0,B,0),\;
\vec C = (0,0,C)$.  

Notation: $\Delta A$ is the area of the triangle $\vec A \vec B \vec C$ and
$\Delta A_x$ is the area of the triangle $\vec B \vec C \vec 0$.

Consider the forces acting on the (infinitessimally small) tetrahedron: 

$$-\vec t_x \Delta A_x - \vec t_y \Delta A_y - \vec t_z \Delta A_z  + \vec t
\Delta A = 0$$

Note: the reason for the negative signs is that the outer unit normals of
$A_x,A_y,A_z$ are $(-1,0,0), (0,-1,0) $ and $(0,0,-1)$.

\begin{align*}
    \Delta A_x = n_1 \Delta A \\
    \Delta A_y = n_2 \Delta A \\
    \Delta A_z = n_3 \Delta A \\
\end{align*}

\lecture{}

$t_x = \begin{pmatrix} T_{11} \\ T_{12} \\ T_{13}\end{pmatrix} $
$t_y = \begin{pmatrix} T_{21} \\ T_{22} \\ T_{23}\end{pmatrix} $
$t_y = \begin{pmatrix} T_{31} \\ T_{32} \\ T_{33}\end{pmatrix} $

Gives $\vec t = T^T \vec n$  where $T$ is the \emph{Cauchy stress tensor}

    $$ T = \begin{pmatrix} T_{11} & T_{12} & T_{13} \\ T_{21} & T_{22} &
    T_{23} \\ T_{31} & T_{32} & T_{33} \end{pmatrix} $$

    Implies $\vec t$ is give by a linear map .

    $$\vec t(t,\vec x, \vec n) = T^T(t,\vec x) \vec n$$


    And substituting this, the balance law for momentum now reads

    $$\vec d {dt} \iiint_{R(t)} \varrho\vec v d V(\vec x) = \iint_{\p R(t) } T^T
    \vec n \: dS + \iiint_{R(t)} \varrho\vec f \: dV$$

    In differential form 
    
    $$D_t (\varrho \vec v) + \varrho \vec v  (\Delta \cdot \vec v) = \nabla \cdot
    t + \varrho \vec f$$

    Where $\nabla \cdot$ acts column-wise.

    $$\nabla \cdot T = 
     \begin{pmatrix}\p_x T_{11} &\p_y T_{12} &\p_z T_{13} \\ 
                    \p_x T_{21} &\p_y T_{22} &\p_z  T_{23} \\ 
                    \p_x T_{31} &\p_y T_{32} &\p_z T_{33} \end{pmatrix} $$

    


The left hand side of $( * *  ) $ can be written as 

$$\varrho D_t \vec v + \vec v \cdot D_t \varrho + \vec v \varrho(\nabla \cdot
\vec v) = \varrho D_t \vec v + \vec v (D_t \verrho + \varrho (\nabla \cdot \vec
v))$$

Notice the rightmost part is the continuity equation from a couple lectures ago,
which is equal to zero. 

So the momentum equation in material coordinates is given by

$$\varrho D_t \vec v \nabla \cdot T + \varrho \vec f$$

\subsection{Conservation of angular momentum}

The angular momentum of a piece of mass at a position $\vec x$ with velocity
$\vec v$ is given by 
the cross product of $\vec v$ and momenum $m\vec v$.

$$\vec l = \vec x \times (\vec v m )$$

\marginnote{$\vec a \times \vec b = \begin{pmatrix} a_2 b_3 - a_3 b_2 \\ a_3 b_1 - a_1 b_3 \\ a_1 b_2 - a_2
b_1 \end{pmatrix}$}

Write the same balance law as for linear momentum but take the cross product of
everything with the position. 

$$\frac d {dt}\iiint_{R(t)} \vec x \times (\varrho \vec v) \: dV(\vec x) =
\iint_{\p R(t)} \vec x \times (T^T \vec n) \; dS(\vec x) + \iiint_{R(t)} \vec x
\times (\varrho f) \: dV(\vec x)$$

Use notation $T =  \begin{pmatrix} \vec t_x \\ \vec t_y \\ \vec t_z
\end{pmatrix}  = (\vec t_1 | \vec t_2 | \vec t_3)$
 and so $T^T \vec n = \begin{pmatrix} \vec t_1 \cdot \vec n \\  \vec t_2 \cdot \vec n \\  \vec t_3 \cdot \vec n \end{pmatrix} $

 Writing the frist component of the balance law only:

$$
\begin{gathered}
\frac{d}{d t} \iiint_{R(t)} \varrho\left(x_{2} v_{3}-x_{3} v_{2}\right) d
V-\iiint_{R(t)} \varrho\left(x_{2} f_{3}-x_{3} f_{2}\right) d V- \\
-\iint_{\partial R(t)}\left(x_{2} \vec{t}_{3} \cdot \vec{n}-x_{3} \vec{t}_{2} \cdot \vec{n}\right) d S=0
\end{gathered}
$$

Or in differential form

$x_1 \cdot $ (3ed component of $(*) ) - x_3 ($ 2nd component of $(*)) - T_{23} +
T_{32} \implies T_{23} = T_{32}$ and if you do the other components of the
equation
$T_{13} = T_{31}$ and $T_{12} = T_{21}$, so we find $T$ is symmetric.

\begin{itemize}
    \item As a consequence of the conservation of angular momentum, $T$ must be
        a symmetric matrix.
\end{itemize}

So the full system we have so far is

\begin{itemize}
    \item Continuity equation
    $$\p_t \rho + \nabla (\rho \vec v) = 0$$
\item $$\rho D_t \vec v = \nabla \cdot T + \rho \vec f,\quad T = T^T$$
\end{itemize}

What is missing is the constitutive relation for $T$

\lecture{}
\subsection{Material frame indifference}

Any constitutive law should not depend on translation or rotation: any rigid body motion. 


$Q(t) \vec c + b(t) $ where  $b(t) \in \R^3$  and $Q(t) \in \R^{3\times 3}$ are
smooth functions, and $Q$ is a rotation matrix. $Q Q^T = \mathbb I \forall t >
0$, and $\det Q = 1$.

\begin{enumerate}
    \item Objectivity: the stress under one coordinate frame according to the
        constituvie law is the same in the
        other coordinate frame. You should be able to compute the stress tensor
        in either and get the same result. 

    \item  Form invariance: any constitutive law $T = G(T) $ stress tensor in
        terms of another quantity $R$ should be the same in the transformed
        coordinate frame. 
\end{enumerate}

These conditions give:

$$G(R^*) = T^* = Q G(R) Q^T \qquad\text{For all } R \text{ and } T$$


\begin{ex}
    $X^* = QX + b(t) \quad |_{\p t}$

    $V^* - Q'(t) X + QV + b'(t) $ this is how the velocity field transforms.

    And a constitutive law such as $T = G(v)$ now in spatial coordinates. (We are
    using $R = \vec v$ and $R^* = \vec v^*$).

    MFI in this case is
    
    $$Q G(v) Q^T = G(Q'(t) x + Q v + b'(t))$$

    For all rotations $Q(t)$ and shifts $b(t)$, but even if $Q = \mathbb I$
    $C(v) = G(v + b'(t))$ is only true if $b'(t) = 0$ ie $b$ is constant.


\end{ex}

\begin{ex}
    $T = \mu (\nabla v + (\nabla v)^T)  $ ie $R = \nabla v + \nabla v ^T$ is
    materially frame indifferent. 

\end{ex}

\begin{thm}{}{} Rivilin Ericson Theorem

    Assuming that $T$ and $R$ are symmetric and objective and $T = G(R)$ is form
    invariant, then it can be written as 

    $$T = \alpha_0 \mathbb I + \alpha_1 R + \alpha_2 R^2$$

    Where $\alpha_0, \alpha_1,\alpha_2 \in \R$ are functions of 

    $I_R = tr(R) =
    \lambda_1 + \lambda_2 + \lambda_3$

    $II_R = \frac 12(tr(R)^2 - tr(R^2)) = \lambda_1 \lambda_2 + \lambda_1
    \lambda_3 + \lambda_2 \lambda_3$

    $III_R = \det(R) = \lambda_1 \lambda_2 \lambda_3$

    Which are the principle invariants of $R$, ie, they are conserved under
    rotations and $\lambda_1 , \lambda_2 \lambda_3$ are eigenvalues of $R$. 

    How to show this: 

    Either constrctive by picking a range of potential rotations, or assuming
    that $G$ has a (taylor) expansion $G(R) = \sum_{n = 0}^\infty \kappa_n R^n$
    with constants $\kappa_n$, and using the Cayley-Hamilton theorem: 

    $R \in \R^{3\times 3} $ satisfies its characteristic equation: $\det(R - X
    \mathbb I) = 0 \implies R^3 I_R R^2 - II_R + III_R \mathbb I$.

$... = \alpha_{2,n} R^2 + \alpha_{1,n} R + \alpha_{0,n} \mathbb I$.
\end{thm}

\week{}
\lecture{}

\subsection{Newtonian Fluids}

The stress depends linearly on the strain rate. (In 1D strain rate was $\p_t
\varepsilon$). In 3D strain rate is pressure. 

\begin{itemize}
    \item Pressure is modelled by $T = -p \mathbb I$. Acting through an
        arbitrary plane with unit normal $\vec n$, it is given by $-p \vec n$.
        Pressure is isotripic. 
    \item Viscous stress. In 1D this was $\varepsilon_t = \p_t (U_A) = V_A$
        to model relative velocity differences between neighbouring particles.
        In 3D spatial coordinates an analog would be $\p_x v, \p_y v, \p_z v$ ie
        $\nabla \vec v$.

        Decompose this into its symmetric and anti-symmetric component: $\nabla
        \vec v = D + W $ where:

        $D$ is the rate of deformation tensor defined as:
$$
\begin{aligned}
{D} & \equiv \frac{1}{2}\left(\nabla \vec{v}+\nabla \vec{v}^{T}\right) \\
&=\left(\begin{array}{ccc}
\frac{\partial v_{1}}{\partial x} & \frac{1}{2}\left(\frac{\partial v_{1}}{\partial y}+\frac{\partial v_{2}}{\partial x}\right) & \frac{1}{2}\left(\frac{\partial v_{1}}{\partial z}+\frac{\partial v_{3}}{\partial x}\right) \\
\frac{1}{2}\left(\frac{\partial v_{1}}{\partial y}+\frac{\partial v_{2}}{\partial x}\right) & \frac{\partial v_{2}}{\partial y} & \frac{1}{2}\left(\frac{\partial v_{2}}{\partial z}+\frac{\partial v_{3}}{\partial y}\right) \\
\frac{1}{2}\left(\frac{\partial v_{1}}{\partial z}+\frac{\partial v_{3}}{\partial x}\right) & \frac{1}{2}\left(\frac{\partial v_{2}}{\partial z}+\frac{\partial v_{3}}{\partial y}\right) & \frac{\partial v_{3}}{\partial z}
\end{array}\right)
\end{aligned}
$$

And $W$ is the vorticity or spin tensor:
$$
\begin{aligned}
{W} & \equiv \frac{1}{2}\left(\nabla \vec{v}-\nabla \vec{v}^{T}\right) \\
&=\left(\begin{array}{ccc}
0 & \frac{1}{2}\left(\frac{\partial v_{1}}{\partial y}-\frac{\partial v_{2}}{\partial x}\right) & \frac{1}{2}\left(\frac{\partial v_{1}}{\partial z}-\frac{\partial v_{3}}{\partial x}\right) \\
\frac{1}{2}\left(\frac{\partial v_{2}}{\partial x}-\frac{\partial v_{1}}{\partial y}\right) & 0 & \frac{1}{2}\left(\frac{\partial v_{2}}{\partial z}-\frac{\partial v_{3}}{\partial y}\right) \\
\frac{1}{2}\left(\frac{\partial v_{3}}{\partial x}-\frac{\partial v_{1}}{\partial z}\right) & \frac{1}{2}\left(\frac{\partial v_{3}}{\partial y}-\frac{\partial v_{2}}{\partial z}\right) & 0
\end{array}\right)
\end{aligned}
$$
\end{itemize}

This suggests a constitutive law of the form:

$$T = -p \mathbb I + G(D,w) \qquad(*)$$
    
However under rigid body motion $X^* = Q(t) X + b(t)$ and $D$ and $w$ transform
according to  $D^* = QDQ^T$ $w^* = QwQ^T +
\underbrace{Q'Q^T}_{\text{an antisimmetric
matrix}}$

$(*)$ is MFI: $QG(R)Q^T = G(R^*)$ if 

$$QG(D,w) Q^T = G(D^*, w^*) = G(Q D Q^T, QwQ^T, + Q'Q^T)\quad \forall Q(t)$$

But if $Q = \mathbb I$ and $Q' = M$ (anti-symmetric), 

$$G(D,w) = G(D, w + M)$$

Which is only true if $M = \hat 0 \in \R^{3 \times 3}$ is the zero matrix. 
$\implies T = -p \mathbb I + G(D)$. 

The rivilin-erickson theorem now tells us that 

$$G(D) = \alpha_0 \mathbb I + \alpha_q D + \alpha_2 D^2$$

And $\alpha_0 = 0$ if $D = \mathbf 0$, there is no viscous stress if $D = \mathbf 0$.

Finally, since linearity characterises newtonian fluids, we simplify the model
by assuming linear dependence of $T$ on $D$.

$\implies \alpha_2 = 0$ because $D^2$ is not linear. 

$\alpha_1 = 2\mu$, a constant, since $D$ is linear. 

$\alpha_0 = \lambda I_D = \lambda tr(D) = \lambda \nabla \cdot \vec v$ since
$I_d$ is the only principle invariant that is linear.     


\begin{itemize}
    \item  [$\mu$]  Dynamic shear velocity 
    \item [$p$] pressure
    \item [$\lambda$] bulk viscosity
\end{itemize}

$\implies \nabla \cdot T = - \nabla p + \lambda \nabla (\nabla \cdot \vec v )
+ \mu(\Delta \vec v + \nabla (\nabla \cdot \vec v))$

And substitute into the momentum equation to obtain the Navier-Stokes equation. 

$$\varrho D_t \vec v = -\nabla p (\lamba + \mu) \nabla (\nabla \cdot \vec v) +
\mu \Delta \vec v + \varrho \vec f$$

Typically couple to the continuity equation:

$$D_t \varrho + \varrho \nabla \cdot \vec v = 0$$

Ideal gas law (with known temperature $T$):

$$p = \varrho R T$$

Usually this system is solved numerically, for example using OpenFoam software
package. 

\lecture{}

A common simplification is to assume incompressibility:  $\p_t \varrho = 0 =
\p_{\vec x_i } \varrho \implies \p_t \varrho + \Nabla \cdot (\varrho \vec v) = 0$

$\implies \nabvla \cot \vec v = 0$

The stress in this case is $T = -p \I + 2\mu D$ and we obtain the incompressible
navier-stokes equation:

$$\varrho D_t \vec v = - \nabla p + \mu \Delta \vec v + \varrho \vec f$$

Coupled to : $\nabvla \cdot \vec v = 0$

And treat $p$ as lagrangian multiplier to enforce $\nabla\cdot \vec v = 0 $ 

Then for an incompressible newtonian fluid in cartesian coordinates, letting
$\vec v = (u,v,w)$ and $\vec f = (f,g,h)$:

\begin{gathered}
\rho\left(\frac{\partial u}{\partial t}+u \frac{\partial u}{\partial x}+v \frac{\partial u}{\partial y}+w \frac{\partial u}{\partial z}\right)=-\frac{\partial p}{\partial x}+\mu\left(\frac{\partial^{2} u}{\partial x^{2}}+\frac{\partial^{2} u}{\partial y^{2}}+\frac{\partial^{2} u}{\partial z^{2}}\right)+\rho f \\
\rho\left(\frac{\partial v}{\partial t}+u \frac{\partial v}{\partial x}+v \frac{\partial v}{\partial y}+w \frac{\partial v}{\partial z}\right)=-\frac{\partial p}{\partial y}+\mu\left(\frac{\partial^{2} v}{\partial x^{2}}+\frac{\partial^{2} v}{\partial y^{2}}+\frac{\partial^{2} v}{\partial z^{2}}\right)+\rho g \\
\rho\left(\frac{\partial w}{\partial t}+u \frac{\partial w}{\partial x}+v \frac{\partial w}{\partial y}+w \frac{\partial w}{\partial z}\right)=-\frac{\partial p}{\partial z}+\mu\left(\frac{\partial^{2} w}{\partial x^{2}}+\frac{\partial^{2} w}{\partial y^{2}}+\frac{\partial^{2} w}{\partial z^{2}}\right)+\rho h \\
\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}=0
\end{gathered}


Note: $D_t \vec v = \p_t \vec v + (\vec v \cdot \nabla) \vec v = \p_t
\vec v + (\nabla\vec v) \vec v $  

$\Delta \vec v = \begin{pmatrix} \Delta u\\ \Delta v \\ \Delta w \end{pmatrix} $

 
\subsubsection{Remarks}

\begin{itemize}
    \item 
Proving global existance in time for the incompressible navier stokes
equation is one of the millenium problems. 
\item
This model is the basis for
simulations of fluids on all spatial scales: galaxies, earths mantle, oceans,
hydraulics, in biology; blood flow, tissue, cytoplasm, etc. 
\item 
    $\mu$ is the shear viscosity: 

    \begin{center}
\begin{array}{l|l|l}
\hline \text { Fluid } & \text { Viscosity (Pa-s) } & \text { Density }\left(\mathrm{kg} / \mathrm{m}^{3}\right) \\
\hline \text { Air } & 1.8 \times 10^{-5} & 1.18 \\
\hline \text { Water } & 8.9 \times 10^{-4} & 0.997 \times 10^{3} \\
\hline \text { Mercury } & 1.5 \times 10^{-3} & 1.3 \times 10^{4} \\
\hline \text { Olive oil } & 0.08 & 0.92 \times 10^{3} \\
\hline \text { Honey } & 37 & 1.4 \times 10^{3} \\
\hline \text { Pitch } & 2.3\times 10^8&                   \\
\hline
\end{array}
\end{center}
\end{itemize}

\section{Solutions to steady flow problems}

Steady flow is characterised by $\p_t \vec v = 0$, and $\p_t p = 0$ and $\p_t
\varrho = 0$.

This implies that $D_t \vec v = \p_t \vec v + (\vec v \cdot \nabla) \vec v$ 

So we now the Newtonian Navier--Stokes becomes:

$$\begin{cases} \varrho(\vec v \cdot \nabla) \vec v = - \nabla p + \mu \Delta
\vec v + \varrho \vec f \\ \nabla \vec v = 0\end{cases}$$

Note that $\vec v, p \varrho $ do not change with time but single particles
still move, and their trajectories must satisfy in general:

$$\begin{cases}\frac {d X} {dt} = V(A,t) = v(X(A,t), t)\\ X(A,t=0) = A
\end{cases}$$

So in steady flow:

$$\begin{cases} \frac {dX} {dt} = V(A) = v(X(A,t)) \\ X(A,t=0) = A\end{cases}$$

The trajectory lines do not change with time. 


\subsection{Coette plane flow}

Steady flow of fluid between two parallel planes. 

\incfigw{coette}{0.6}

Lower plate is stationary, upper plate has speed in x direction of $u_0$. 
$dt y = h$.

Use the boundary condition: assume there is a no-slip condition between the fluid
and the boundary with the plate: the fluid travels at the speed of the plate. 

At the upper plate $\vec v (x, y = h, z) =  (u_0 , 0 , 0)$

At the lower plate $\vec v(x, y= 0, z) = (0,0,0)$

Assume laminar flow: $\vec v = (u,0,0)$: flow is only in the $x$ direction,
solution is homogeneous in $z$-direction: $\p_z u = 0, \p_z p  = 0$

Now our model becomes:

$$
(* * *)
\begin{cases}
   u_x = 0 \\
   \varrho u_x u = -p_x + \mu (u_{x x} + u_{y y } + u_{z z}) \\
   0 = -py + 0  \\
   0 = 0 + 0 \\
\end{cases}
$$

$\p_z u = 0 = \p_x u \implies u = u(y)$

$\p_z p = 0 = \p_y p \implies p = p(x)$

So $(* * * ) $ becomes $0 = -p_x + \mu u_{y y}$, i.e., $p'(x) = \mu u''(y)$.

Use method of separation: $p'(x) = const = \mu u''(y) \implies p' = const
\implies$ 

$$p(x) = p_0 + x p_1$$

Modelling intuition: $p_1 \ne 0$ would impl;y that the pressure blows up at $x =
\pm \infty \implies $ assume that $p_1 = 0$, with that $p(x) = p_0$.

\lecture{}

The second equation becomes $\mu u''(y) = 0 \implies u(y) = ay + b$. The
boundary conditions are $u(y = 0) = 0, u(h = 0) = u_0$, so $a = \frac {u_0} h
\implies u = \gamma y$ where $\gamma = \frac {u_0} h$ 

$\gamma = \frac {u_0} {h} $ is the \emph{shear rate} which has dimension $\frac
1t$.

So we have a solution: 

$\vec v = \begin{pmatrix} \gamma y \\ 0 \\0 \end{pmatrix}$  and $p = p_0$

Remark: The stress tensor is given by 

$$T = -p \mathbb I + 2\mu D = -p_0 \mathbb I + \mu(\nabla v + {\nabla v}^T) = -p_0
         \mathbb I + \mu \begin{pmatrix}0 & \gamma & 0 \\ \gamma & 0 & 0 \\ 0 & 0 & 0   \end{pmatrix} $$

         So the stress through a horizontal plane, for example $ \{ y = 0\} ,
         \vec n = (0,1,0)$, is given by


         $$
         T \vec n = \left(-p_0 \mathbb I + \mu \begin{pmatrix} 0 & \gamma & 0 \\
         \gamma & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}  \right) \begin{pmatrix} 0 \\ 1 \\
         0\end{pmatrix} = \underbrace{-p_0 \begin{pmatrix} 0\\1\\0
 \end{pmatrix}}_{\vec t_s}  +
     \underbrace{\mu
 \begin{pmatrix} \gamma \\ 0 \\ 0 \end{pmatrix} }_{\vec t_p}
         $$

         $\vec t_s$ --- the shear component.

         $\vec t_p$ --- the normal component.

\subsubsection{Shear Test}

A Shear test can be used to investgate whether a fluid is Newtonian (where
stress is linear with respect to $\gamma$) and what the shear viscosity is. 

Put the fluid between two plates and move the upper plate with differnet speeds, 
$\gamma = \frac {u_0}{h}$. Measure the force (stress $\times$ area), to see if
it is linear with respect to the velocity of the top plate.

Since shear stress $\vec t_s = \mu \gamma \begin{pmatrix} 1 \\ 0 \\ 0
\end{pmatrix} $, the gradient of your shear stress / shear rate graph is the
shear viscosity $\mu$.

\subsection{Poiseuille Flow}

Steady flow through a cylindrical pipe of length $L$ with radius $R$. 

Use the incompressible Navier-Stokes equation written in cylindrical
coordinates. 

Write the problem in polar coordinatnes about the centre of the pipe, with $z$
running along its length: 

$\vec v = (v_R, v_\theta, v_z) =  v_r \vec e_r + v_\theta \vec e_\theta
+ v_z + \vec e_z$, and 

$\vec f = f_r \vec e_r + f_\theta \vec e_\theta + f_z \vec e_z$, 

$x = r\cos \theta$ and $y = r \sin \theta$ and $z = z$. 

The incompressible Navier--Stokes system for newtonian fluid in
cyndrical coordinates is:

\begin{aligned}
&\rho\left(\frac{\partial v_{r}}{\partial t}+v_{r} \frac{\partial v_{r}}{\partial r}+\frac{v_{\theta}}{r} \frac{\partial v_{r}}{\partial \theta}-\frac{v_{\theta}^{2}}{r}+v_{z} \frac{\partial v_{r}}{\partial z}\right)\\
&\quad=-\frac{\partial p}{\partial r}+\mu\left[\frac{\partial}{\partial r}\left(\frac{1}{r} \frac{\partial}{\partial r}\left(r v_{r}\right)\right)+\frac{1}{r^{2}} \frac{\partial^{2} v_{r}}{\partial \theta^{2}}-\frac{2}{r^{2}} \frac{\partial v_{\theta}}{\partial \theta}+\frac{\partial^{2} v_{r}}{\partial z^{2}}\right]+\rho f_{r}\\\\
&\rho\left(\frac{\partial v_{\theta}}{\partial t}+v_{r} \frac{\partial v_{\theta}}{\partial r}+\frac{v_{\theta}}{r} \frac{\partial v_{\theta}}{\partial \theta}+\frac{v_{r} v_{\theta}}{r}+v_{z} \frac{\partial v_{\theta}}{\partial z}\right)\\
&\quad=-\frac{1}{r} \frac{\partial p}{\partial
\theta}+\mu\left[\frac{\partial}{\partial r}\left(\frac{1}{r}
\frac{\partial}{\partial r}\left(r v_{\theta}\right)\right)+\frac{1}{r^{2}}
\frac{\partial^{2} v_{\theta}}{\partial \theta^{2}}+\frac{2}{r^{2}}
\frac{\partial v_{r}}{\partial \theta}+\frac{\partial^{2} v_{\theta}}{\partial
z^{2}}\right]+\rho f_{\theta}\\ \\
&\rho\left(\frac{\partial v_{z}}{\partial t}+v_{r} \frac{\partial v_{z}}{\partial r}+\frac{v_{\theta}}{r} \frac{\partial v_{z}}{\partial \theta}+v_{z} \frac{\partial v_{z}}{\partial z}\right)\\
&\quad =-\frac{\partial p}{\partial z}+\mu\left[\frac{1}{r}
\frac{\partial}{\partial r}\left(r \frac{\partial v_{z}}{\partial
r}\right)+\frac{1}{r^{2}} \frac{\partial^{2} v_{z}}{\partial
\theta^{2}}+\frac{\partial^{2} v_{z}}{\partial z^{2}}\right]+\rho f_{z}\\  \\
&\frac{1}{r} \frac{\partial\left(r v_{r}\right)}{\partial r}+\frac{1}{r} \frac{\partial v_{\theta}}{\partial \theta}+\frac{\partial v_{z}}{\partial z}=0
\end{aligned}
  
Unlike before there is pressure difference across the length ($z$) of the pipe 
and the flow is coming from the pressure difference. 


The boundary conditions are $\vec v = \begin{pmatrix} 0\\0\\0 \end{pmatrix} $ at
$r = R$ (no-slip boundary condition). $p(z = 0) = p_0$ and $p(z = L) = p_1 <
p_0$.

$\vec v = \begin{pmatrix} v_r \\ v_\theta \\ v_z \end{pmatrix} $

Again assume laminar flow (trajectory not depending on $\theta$ or $r$).

The velocity only has a $z$ component: 

$v_r(z = 0) = 0$ and $v_\theta(z = 0) = 0$

$v_r(z = L) = 0$ and $v_\theta (z = L) = 0$

$v_\tehta = v_r = 0 \implies \vec v = \begin{pmatrix} 0 \\ 0 \\ v_z
\end{pmatrix} $

And assume the solution satisfies rotational symmetry 

$v_z = v_z(r,z)$ and $p = p(r,z)$.

And assume body forces are $0$.

Now we can get rid of most of the Navier Stokes terms because of these
assumptions which leaves the governing equations

$$
\begin{cases}
    0 = \p_r p \\
    \varrho v_z \frac {\p v_z} {\p z} =  -p_z + \mu (\frac 1r \p_r( r \p_r
    v_z ) + \p _{z z} v_z)\\
    0 = \p_z v_z 
\end{cases}
$$


$\implies p = p(z) $ and $v_z = (r)$  and $p'(z) = \mu(v_z'' (r) + \frac 1r v_z'
(r)$

Use method of separation: $p'(z) = constant = \mu (v_z '' + \frac 1r v_z')$

Use $p (z = 0) = p_0$ and $p(z = L) = p_1$ $\implies$

$p = p_0 + \alpha z$ and $l(L) = p_1 = p_0 + \alpha l \implies \alpha = \frac
{p_1 - p_0 } L$ 

$\implies p(z) = p_0 + \frac {p_1 - p_0}L z$

$\frac {p_1 - p_0} L = \mu (v_z'' + \frac 1r v_z ' ) $

$$\frac r \mu \frac {p_1 - p_0} L = r v_z '' + vz' = (rvz') ' $$

Integrate

$$\frac {r^ 2} 2 \frac {p_1 - p_0} {\mu L} = r v_z ' - C$$

Integrate

$$v_z ' = \frac r2 \frac {p_1 - p_0} {\mu L} + \frac 1r C$$

Integrate using $v_z(R) = 0$

$$v_z = (\frac {r^2}2 - \frac {R^2} 2) \frac {p_1 - p_0} {2\mu L} + C \log r$$

If $C \ne 0$ the velocity would blow up towards the centre of the pipe, so $C =
0$. 

 $$\implies v_z = 
 \begin{pmatrix} 0 \\ 0 \\ R^2 - r^2 \end{pmatrix} 
 \frac {p_0 - p_1} {4 \mu L} $$

 So you get a parabolic velocity profile. 

 The stress tensor in cylindrical coordinates is given by: 

     $$T = \begin{pmatrix}  -p & 0 & \mu v_z' \\ 0 & -p & 0 \\
     \mu v_z ' & 0 & -p \end{pmatrix} $$

     Eg the stress at the surface of a cylinder with radius $r < R$, with normal
     to the cylinder of $\vec n = \begin{pmatrix} 1\\0\\0 \end{pmatrix}
     \implies$ $T \vec n = \begin{pmatrix} -p \\ 0 \\ 0 \end{pmatrix}  + \begin{pmatrix}
     0 \\ 0 \\ \mu v_z '\end{pmatrix} = \vec t_p + \vec t_s$, shear stress and
     normal stress.

\section{Closing remarks}

\begin{enumerate}
    \item Vorticity: the curl of the velocity at a specific point

        $$\vec \omega = \nabla \times \vec v = \begin{pmatrix} \frac {\p
        w} {\p y} - \frac {\p v} {\p t} \\ \frac {\p v} {\p t} - \frac {\p w }
    {\p x} \\ \frac {\p v} {\p x} - \frac {\p u} {\p y}  \end{pmatrix} $$

    It is like `how much the particle rotates'. 

    For example in plane Couette flow: $\vec v = \begin{pmatrix} \gamma y \\ 0
    \\ 0 \end{pmatrix} $, the vorticity is $\vec \omega = \nabla \times \vec v
    = \begin{pmatrix}  0 \\ 0 \\ - \gamma\end{pmatrix} $

\item 
    Helmholtz Representation Theorem

    Let $\vec q (\vex x) : \R^3 \to \R^3$ diffble, then 

    $\vec q (\vec x) = \nabla \phi + \nabla\times \vec g $ where 

    $\phi:\R^3 \to
    R$ and $\vec g: \R^3 \to \R^3, \nabla \cdot \vec g = 0$.

    It can be split into the scalar potential with divergence and no rotational
    component (curl) and
    the vector potential that is pure rotation. 


\end{enumerate}


\week{}
\lecture{}


\pagebreak
\appendix

\chapter{Laplace Tranform Review}\label{laplace}

\marginnote{Taken from MATH2100 Workbook}

$$
F(s)=L(f(t))=\int_{0}^{\infty} e^{-s t} f(t) d t
$$
Powers
$$
L\left(t^{n}\right)=\frac{n !}{s^{n+1}} \quad L\left(t^{a}\right)=\frac{\Gamma(a+1)}{s^{a+1}}
$$
Sine and Cosine
$$
L(\cos (\alpha t))=\frac{s}{s^{2}+\alpha^{2}} \quad \text { and } \quad L(\sin (\alpha t))=\frac{\alpha}{s^{2}+\alpha^{2}}
$$
Dirac Delta function
$$
L(\delta(t-k))=e^{-k s}
$$
Linearity
$$
L(a g(t)+b h(t))=a G(s)+b H(s)
$$
First Shifting Theorem
$$
L\left(e^{a t} f(t)\right)=F(s-a)
$$
Second Shifting Theorem
$$
L(f(t-k) u(t-k))=e^{-k s} F(s)
$$
Convolution Theorem
$$
L\left(\int_{0}^{t} f(\tau) g(t-\tau) d \tau\right)=F(s) G(s)
$$
Differentials
$$
\begin{gathered}
L\left(f^{(n)}(t)\right)=s^{n} F(s)-s^{n-1} f(0)-s^{n-2} \dot{f}(0)-\ldots-s f^{n-2}(0)-f^{n-1}(0) \\
L(t f(t))=-\frac{d F(s)}{d s}
\end{gathered}
$$
Periodic Functions $f(t+p)=f(t)$
$$
L(f(t))=\frac{1}{1-e^{-s p}} \int_{0}^{p} e^{-s t} f(t) d t
$$

\chapter{Misc}

\section{Cross Product}





\end{document}
