You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
3710 lines
127 KiB
3710 lines
127 KiB
|
|
\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}
|
|
|