diff --git a/tutorial/julia/Tutorial 4 - ODEs, Boundary Value Problems/.ipynb_checkpoints/J4_ODEBVPsystems_solved-checkpoint.ipynb b/tutorial/julia/Tutorial 4 - ODEs, Boundary Value Problems/.ipynb_checkpoints/J4_ODEBVPsystems_solved-checkpoint.ipynb
new file mode 100644
index 0000000..1ffefc4
--- /dev/null
+++ b/tutorial/julia/Tutorial 4 - ODEs, Boundary Value Problems/.ipynb_checkpoints/J4_ODEBVPsystems_solved-checkpoint.ipynb
@@ -0,0 +1,433 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Ordinary Differential Equations, Boundary Value Problems \n",
+ "\n",
+ " \n",
+ "\n",
+ "# Learning Objectives \n",
+ "\n",
+ "### ODE boundary value problems (BVPs)\n",
+ "- Be able to distinguish bewteen ODE IVPs and BVPs\n",
+ "- Be able to formulate three common types of boundary conditions: Dirichlet, Neumann, Robin\n",
+ "\n",
+ "### Finite difference method (FDM)\n",
+ "- Be able to derive forward, backward, and centered finite difference approximations for the first derivatives\n",
+ "- Be able to derive centered finite difference approximations for the second derivatives\n",
+ "- Be able to state the truncation errors for different finite difference approximations\n",
+ "\n",
+ "### Solving ODE BVPs\n",
+ "- Demonstrate the ability to discretize ODE BVPs into a set of algebraic equations\n",
+ "- Be able to deal with Dirichlet boundary conditions\n",
+ "- Be able to deal with Neumann and Robin boundary conditions by making fictious nodes\n",
+ "- Solve algebraic equations using appropriate numerical methods (Gauss-Siedel for linear systems or Newton-Raphson for nonlinear systems)\n",
+ "\n",
+ "### Solving coupled BVPs*\n",
+ "- Be able to solve coupled BVPs by interlacing the unknown dependent variables\n",
+ "- Be able to identify the bandwidth of the Jacobin matrix of the reformulated equations\n",
+ "\n",
+ "* Not covered in tutorial\n",
+ "\n",
+ " "
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "
\n",
+ "INTERACTIVE! Before starting this example, import the Plots.jl library by running the cell below.\n",
+ "
"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "\u001b[32m\u001b[1m Updating\u001b[22m\u001b[39m registry at `C:\\Users\\stub0\\.julia\\registries\\General.toml`\n",
+ "\u001b[32m\u001b[1m Resolving\u001b[22m\u001b[39m package versions...\n",
+ "\u001b[32m\u001b[1m No Changes\u001b[22m\u001b[39m to `C:\\Users\\stub0\\.julia\\environments\\v1.7\\Project.toml`\n",
+ "\u001b[32m\u001b[1m No Changes\u001b[22m\u001b[39m to `C:\\Users\\stub0\\.julia\\environments\\v1.7\\Manifest.toml`\n",
+ "\u001b[32m\u001b[1mPrecompiling\u001b[22m\u001b[39m project...\n",
+ "\u001b[32m ✓ \u001b[39m\u001b[90mStrideArraysCore\u001b[39m\n",
+ "\u001b[32m ✓ \u001b[39m\u001b[90mIntervalContractors\u001b[39m\n",
+ "\u001b[32m ✓ \u001b[39m\u001b[90mPolyester\u001b[39m\n",
+ "\u001b[32m ✓ \u001b[39m\u001b[90mGroebner\u001b[39m\n",
+ "\u001b[32m ✓ \u001b[39m\u001b[90mFastBroadcast\u001b[39m\n",
+ "\u001b[32m ✓ \u001b[39m\u001b[90mDistributions\u001b[39m\n",
+ "\u001b[32m ✓ \u001b[39m\u001b[90mSciMLBase\u001b[39m\n",
+ "\u001b[32m ✓ \u001b[39m\u001b[90mTriangularSolve\u001b[39m\n",
+ "\u001b[32m ✓ \u001b[39m\u001b[90mDEDataArrays\u001b[39m\n",
+ "\u001b[32m ✓ \u001b[39m\u001b[90mRecursiveFactorization\u001b[39m\n",
+ "\u001b[32m ✓ \u001b[39m\u001b[90mNonlinearSolve\u001b[39m\n",
+ "\u001b[32m ✓ \u001b[39m\u001b[90mLinearSolve\u001b[39m\n",
+ "\u001b[32m ✓ \u001b[39m\u001b[90mMathOptInterface\u001b[39m\n",
+ "\u001b[32m ✓ \u001b[39m\u001b[90mDiffEqBase\u001b[39m\n",
+ "\u001b[32m ✓ \u001b[39mIpopt\n",
+ "\u001b[32m ✓ \u001b[39m\u001b[90mCbc\u001b[39m\n",
+ "\u001b[32m ✓ \u001b[39mGurobi\n",
+ "\u001b[32m ✓ \u001b[39m\u001b[90mSymbolics\u001b[39m\n",
+ "\u001b[32m ✓ \u001b[39mJuMP\n",
+ "\u001b[32m ✓ \u001b[39m\u001b[90mDiffEqCallbacks\u001b[39m\n",
+ "\u001b[32m ✓ \u001b[39m\u001b[90mDiffEqJump\u001b[39m\n",
+ "\u001b[32m ✓ \u001b[39m\u001b[90mBoundaryValueDiffEq\u001b[39m\n",
+ "\u001b[32m ✓ \u001b[39m\u001b[90mSundials\u001b[39m\n",
+ "\u001b[32m ✓ \u001b[39m\u001b[90mDiffEqNoiseProcess\u001b[39m\n",
+ "\u001b[32m ✓ \u001b[39mDiffEqGPU\n",
+ "\u001b[32m ✓ \u001b[39m\u001b[90mSteadyStateDiffEq\u001b[39m\n",
+ "\u001b[32m ✓ \u001b[39m\u001b[90mMINLPTests\u001b[39m\n",
+ "\u001b[32m ✓ \u001b[39mModelingToolkit\n",
+ "\u001b[32m ✓ \u001b[39mEAGO\n",
+ "\u001b[32m ✓ \u001b[39mOrdinaryDiffEq\n",
+ "\u001b[32m ✓ \u001b[39m\u001b[90mDelayDiffEq\u001b[39m\n",
+ "\u001b[32m ✓ \u001b[39m\u001b[90mStochasticDiffEq\u001b[39m\n",
+ "\u001b[32m ✓ \u001b[39mDifferentialEquations\n",
+ " 33 dependencies successfully precompiled in 78 seconds (324 already precompiled)\n"
+ ]
+ }
+ ],
+ "source": [
+ "import Pkg; Pkg.add(\"Plots\"); using Plots"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "In the last tutorial, we introduced the general ODE form,\n",
+ "\n",
+ "\\\\[F\\left(x,y,\\frac{dy}{dx},\\frac{d^2y}{dx^2},\\ldots,\\frac{d^ny}{dx^n}\\right)=0,\\\\]\n",
+ "\n",
+ "which when combined with initial conditions, formed the ODE-IVP problem, that we then proceeded to solve.\n",
+ "\n",
+ "In this unit, we consider ODE problems in which $n$ auxiliary conditions are specified at different values of the independent variable $x$ rather than just the specification of the initial condition. Such conditions are termed boundary conditions and the resulting problems are called boundary value problems (BVPs). In most of the problems of interest in chemical engineering the ODE will be of the form,\n",
+ "\n",
+ "\\\\[\\frac{d^2y}{dx^2}=f\\left(x,y,\\frac{dy}{dx}\\right),\\\\]\n",
+ "\n",
+ "where $x$ belongs to an interval $[x_l,x_r]$. The positions $x=x_l$ and $x=x_r$ define the boundaries of the spatial domain of interest. The three common types of boundary conditions are as follows:\n",
+ "- **Dirichlet** (boundary condition of the first kind): For this type, we specify the value of $y$ on the boundary, e.g. $y(x_l)=\\alpha$.\n",
+ "- **Neumann** (boundary condition of the second kind): For this type, we specify the derivative of $y$ on the boundary, e.g. $\\left. \\frac{dy}{dx} \\right |_{x_l} = \\alpha$.\n",
+ "- **Robin** (boundary condition of the third kind): This type involves a combination of $y$ and the derivative of $y$ on the boundary, e.g. $\\alpha_1 \\frac{dy}{dx} + \\alpha_2 y(x_r) = \\alpha_3$.\n",
+ "\n",
+ "## Finite Difference Method\n",
+ "The finite difference method is a common approach for solving ODE-BVPs. The first step to apply the finite difference method, is to discretize the spatial domain. This begins by defining $n$ (often equidistant) nodes $x_1,\\ldots,x_n$. The discretization step is then given by\n",
+ "\n",
+ "\\\\[\\Delta x = \\frac{x_r-x_l}{n-1}.\\\\]\n",
+ "\n",
+ "First-order derivatives and second-order derivatives appearing in $F$ will be replaced with finite difference approximations. Common finite difference approximations for the first-order derivatives are list below.\n",
+ "- Forward difference approximation: $\\left.\\frac{dy}{dx}\\right |_{x_i}=\\frac{y_{i+1}-y_i}{\\Delta x}$ \n",
+ "- Backward difference approximation: $\\left.\\frac{dy}{dx}\\right |_{x_i}=\\frac{y_{i}-y_{i-1}}{\\Delta x}$ \n",
+ "- Centered finite difference approximation: $\\left.\\frac{dy}{dx}\\right |_{x_i}=\\frac{y_{i+1}-y_{i-1}}{2\\Delta x}$\n",
+ "The centered finite difference approximation is the most common finite difference approximation used for second-order derivatives and is defined as\n",
+ "\n",
+ "\\\\[\\left.\\frac{d^2y}{dx^2}\\right |_{x_i}=\\frac{y_{i+1} +2y_i -y_{i-1}}{\\Delta x^2}.\\\\]\n",
+ "\n",
+ "## Solving ODE-BVPs\n",
+ "After the spatial discretization and the derivation of the finite difference approximations, we can proceed with the transformation of the ODE-BVP into a set of algebraic equations. At each interior node, we use the centered finite difference approximations to convert the ODE $\\frac{d^2y}{dx^2}=f(x,y,\\frac{dy}{dx})$ into the discretized form\n",
+ "\n",
+ "\\\\[\\frac{y_{i+1} - 2y_i +y_{i-1}}{\\Delta x^2}=f(x_i,y_i,\\frac{y_{i+1} - y_{i-1}}{2 \\Delta x}).\\\\]\n",
+ "\n",
+ "The result is a set of algebraic equations in $n-2$ unknowns. The remaining two equations will come from the boundary conditions. The resulting set of (in general, nonlinear) algebraic equations will have to be solved simultaneously to obtain the values of the unknowns. An outline for general approach to solving BVPs is:\n",
+ "1. Discretize the domain into n nodes.\n",
+ "2. For the interior nodes, discretize the ODE using finite difference approximations of the derivatives.\n",
+ "3. For the boundary nodes, if there is a Dirichlet boundary, set the boundary node to the specified value. Otherwise, make a fictitious node and discretize both the boundary condition and the ODE.\n",
+ "4. If the problem is linear, solve using any method the resulting system $\\mathbf A \\mathbf x = \\mathbf b$. If the problem is nonlinear, form the Jacobian and solve by Newton–Raphson.\n",
+ "5. Check for mesh refinement.\n",
+ "\n",
+ "When solving ODE-BVPs, it is very important to apply the appropriate procedure for dealing with each different of type of **boundary condition**.\n",
+ "When a **Dirichlet** boundary condition is present, the value of the boundary node is a specified, e.g. if we have $y(x_l)=\\alpha$, then we set $y_1=\\alpha$.\n",
+ "When a **Neumann** and **Robin** boundary conditions with a derivative term is present, extra consideration is merited. When using the centered finite difference method for approximations, a fictious node should be introduced in order to to keep the same level of accuracy at the boundary. \n",
+ "Consider a Neumann boundary condition at the left boundary $\\left. \\frac{dy}{dx} \\right |_{x_l} = \\alpha$, we need to introduce a fictious node $x_0$ on the left of $x_l$, with a corresponding value of $y(x=x_0)=y_0$. Subsequently, we implement the centered finite difference approximation of the boundary condition:\n",
+ "\n",
+ "\\\\[\\frac{y_2-y_1}{2\\Delta x} = \\alpha,\\\\]\n",
+ "\n",
+ "which we can solve for the value of the fictitious node,\n",
+ "\n",
+ "\\\\[y_0 = y_2 - 2 \\alpha \\Delta x.\\\\]\n",
+ "\n",
+ "We now substitute this result into the differential equation at $x_1$.\n",
+ "\n",
+ "## Case Study: Laminar Flow\n",
+ "Suppose we have laminar flow of water, at 20 °C, through a horizontal cylindrical tube.\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "Assume that we only consider the kinematics in the z-direction, and the flow is fully developed and at steady state. Then we can have the equation:\n",
+ "\n",
+ "\\\\[0=-\\frac{\\partial P}{\\partial z} + \\mu \\left [ \\frac{1}{r} \\frac{\\partial}{\\partial r} \\left(r \\frac{\\partial u_z}{\\partial r} \\right) \\right ]. \\\\]\n",
+ "\n",
+ "We can rearrange and simplify the above equation as\n",
+ "\n",
+ "\\\\[\\frac{1}{r} \\frac{\\partial u_z}{\\partial r} + \\frac{\\partial^2 u_z}{\\partial r^2} = \\frac{1}{\\mu} \\frac{\\partial P}{\\partial z}.\\\\]\n",
+ "\n",
+ "The viscosity of water at $20 ^\\circ \\text C$ is $\\mu = 8.9\\times 10^{-4} \\; \\text{Pa} \\cdot \\text s$, and the pressure drop through the pipe is $\\frac{\\partial P}{\\partial z} = -1000 \\; \\text{Pa}/\\text{m}$. We plug in this information and the equation simplifed to\n",
+ "\n",
+ "\\\\[\\frac{1}{r} \\frac{\\partial u_z}{\\partial r} + \\frac{\\partial^2 u_z}{\\partial r^2} = \\frac{-1000}{8.9\\times 10^{-4}}.\\\\]\n",
+ "\n",
+ "At the inner surface of the pipe, a no-slip boundary condition holds, which is a Dirichlet condition:\n",
+ "\n",
+ "\\\\[u_z|_{r=R} = 0.\\\\]\n",
+ "\n",
+ "At the center of the pipe, a no-flux boundary condition holds due to symmetry, which is a Neumann condition:\n",
+ "\n",
+ "\\\\[\\left. \\frac{du_z}{dr} \\right |_{r=0} = 0.\\\\]\n",
+ "\n",
+ "We may now solve for the $z$-direction velocity $u_z$ profile.\n",
+ "\n",
+ "\n",
+ "## Solution:\n",
+ "We first discretize the spatial domain into $n$ nodes, $r_1,\\ldots,r_n$, with step size $\\Delta r =R/(n-1)$.\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "Then, we use centered finite difference approximations for the derivative terms at each interior node:\n",
+ "\n",
+ "\\\\[\\left. \\frac{d u_z}{dr} \\right|_{r_i} \\approx \\frac{u_{z,i+1} - u_{z,i-1}}{\\Delta r},\\\\]\n",
+ "\n",
+ "\\\\[\\left. \\frac{d^2 u_z}{dr^2} \\right|_{r_i} \\approx \\frac{u_{z,i+1} - 2 u_{z,i} + u_{z,i-1}}{\\Delta r^2}.\\\\]\n",
+ "\n",
+ "Thus, the equations for interior nodes become\n",
+ "\n",
+ "\\\\[\\frac{1}{r_i} \\frac{u_{z,i+1} - u_{z,i-1}}{\\Delta r}+\\frac{u_{z,i+1} - 2 u_{z,i} + u_{z,i-1}}{\\Delta r^2}=\\frac{-1000}{8.9\\times 10^{-4}}.\\\\]\n",
+ "\n",
+ "For the boundary nodes, at $r=R$ where the Dirchlet boundary condition is applied, we have\n",
+ "\n",
+ "\\\\[u_{z,n}=0.\\\\]\n",
+ "\n",
+ "At $r=0$ where the Neumann boundary condition is applied, a fictious node $u_{z,0}$ is introduced, thus we have\n",
+ "\n",
+ "\\\\[\\frac{u_{z,2}-{u_{z,0}}}{\\Delta r}=0.\\\\]\n",
+ "\n",
+ "Solving for the fictious node, we get $u_{z,0} = u_{z,2}$. The equation for the first node has to be modified because the radius is equal to zero at that node. We apply L’Hôpital’s rule:\n",
+ "\n",
+ "\\\\[\\lim_{r \\to 0} \\frac{1}{r} \\frac{d u_z}{dr} = \\lim_{r \\to 0} \\frac{\\frac{d}{dr}\\left(\\frac{du_z}{dr}\\right)}{\\frac{d}{dr}(r)} = \\lim_{r\\to 0}\\frac{d^2 u_z}{dr^2}.\\\\]\n",
+ "\n",
+ "So the equation at $r=0$ becomes $2\\frac{d^2 u_z}{dr^2} = \\frac{-1000}{8.9\\times 10^{-4}}$.\n",
+ "We apply centered finite differencing and plug the fictious node into the first equation to obtain\n",
+ "\n",
+ "\\\\[2\\cdot\\frac{2u_{z,2}-2u_{z,1}}{\\Delta r^2}=\\frac{-1000}{8.9\\times 10^{-4}}. \\\\]\n",
+ "\n",
+ "Now, we have a system of linear equations $\\mathbf A \\mathbf u_z =\\mathbf b$ and we can solve for the velocity profile."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "done\n"
+ ]
+ }
+ ],
+ "source": [
+ "# process specifications and constants\n",
+ "mu = 8.9e-4 # Pa*s\n",
+ "dPdz = -1000 # Pa/m\n",
+ "delta_r = 0.001 # m\n",
+ "r = range(0.0, 0.05, step = delta_r) # m\n",
+ "println(\"done\")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "
\n",
+ "INTERACTIVE! Now, we enter the discretized system of linear algebraic equations in matrix-vector form.\n",
+ "You should write the equation for the first node (center of the pipe) by yourself..\n",
+ "
"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n"
+ ]
+ },
+ "execution_count": 3,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "# creating a matrix and vector for the equations\n",
+ "n = length(r)\n",
+ "A = zeros(n,n)\n",
+ "b = zeros(n)\n",
+ "\n",
+ "### FILL IN THE ENTRIES OF A[1,i] and b[1] for i = 1,...,n by yourself here ###\n",
+ "A[1,1] = -4/delta_r^2\n",
+ "A[1,2] = 4/delta_r^2\n",
+ "b[1] = -1000/8.9e-4\n",
+ "\n",
+ "# interior nodes\n",
+ "for i = 2:n-1\n",
+ " A[i,i-1] = -1/(r[i]*2*delta_r) + 1/delta_r^2\n",
+ " A[i,i] = -2/delta_r^2\n",
+ " A[i,i+1] = 1/(r[i]*2*delta_r) + 1/delta_r^2\n",
+ " b[i] = dPdz/mu\n",
+ "end\n",
+ "\n",
+ "# last node\n",
+ "A[n,n] = 1\n",
+ "b[n] = 0\n",
+ "\n",
+ "# solving the matrix for the velocity profile\n",
+ "u = A\\b\n",
+ "\n",
+ "# plot of velocity profile\n",
+ "plt = plot(u,r,label=\"\")\n",
+ "plot!(plt,u,-r,label=\"\", title= \"Velocity Profile\")\n",
+ "xlabel!(\"Velocity (m/s)\")\n",
+ "ylabel!(\"Radius (m)\")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ " \n",
+ "\n",
+ "# Questions for reflection \n",
+ "\n",
+ "- How can you tell if the spatial discretization chosen for this example was small enough?\n",
+ "- In general, a number of higher-order spatial discretizations may be used when discretizing a problem (see [slide 9 here](https://www.mathematik.uni-dortmund.de/~kuzmin/cfdintro/lecture4.pdf)). What potential trade-offs may occur when using a higher-order accuracy scheme to discretize the problem solved in the=is tutorial? "
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Julia 1.7.2",
+ "language": "julia",
+ "name": "julia-1.7"
+ },
+ "language_info": {
+ "file_extension": ".jl",
+ "mimetype": "application/julia",
+ "name": "julia",
+ "version": "1.7.2"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
![](\"pipe1.png\")
\n",
+ "\n",
+ "\n",
+ "Assume that we only consider the kinematics in the z-direction, and the flow is fully developed and at steady state. Then we can have the equation:\n",
+ "\n",
+ "\\\\[0=-\\frac{\\partial P}{\\partial z} + \\mu \\left [ \\frac{1}{r} \\frac{\\partial}{\\partial r} \\left(r \\frac{\\partial u_z}{\\partial r} \\right) \\right ]. \\\\]\n",
+ "\n",
+ "We can rearrange and simplify the above equation as\n",
+ "\n",
+ "\\\\[\\frac{1}{r} \\frac{\\partial u_z}{\\partial r} + \\frac{\\partial^2 u_z}{\\partial r^2} = \\frac{1}{\\mu} \\frac{\\partial P}{\\partial z}.\\\\]\n",
+ "\n",
+ "The viscosity of water at $20 ^\\circ \\text C$ is $\\mu = 8.9\\times 10^{-4} \\; \\text{Pa} \\cdot \\text s$, and the pressure drop through the pipe is $\\frac{\\partial P}{\\partial z} = -1000 \\; \\text{Pa}/\\text{m}$. We plug in this information and the equation simplifed to\n",
+ "\n",
+ "\\\\[\\frac{1}{r} \\frac{\\partial u_z}{\\partial r} + \\frac{\\partial^2 u_z}{\\partial r^2} = \\frac{-1000}{8.9\\times 10^{-4}}.\\\\]\n",
+ "\n",
+ "At the inner surface of the pipe, a no-slip boundary condition holds, which is a Dirichlet condition:\n",
+ "\n",
+ "\\\\[u_z|_{r=R} = 0.\\\\]\n",
+ "\n",
+ "At the center of the pipe, a no-flux boundary condition holds due to symmetry, which is a Neumann condition:\n",
+ "\n",
+ "\\\\[\\left. \\frac{du_z}{dr} \\right |_{r=0} = 0.\\\\]\n",
+ "\n",
+ "We may now solve for the $z$-direction velocity $u_z$ profile.\n",
+ "\n",
+ "\n",
+ "## Solution:\n",
+ "We first discretize the spatial domain into $n$ nodes, $r_1,\\ldots,r_n$, with step size $\\Delta r =R/(n-1)$.\n",
+ "\n",
+ "\n",
+ "![](\"pipe2.png\")
\n",
+ "\n",
+ "\n",
+ "Then, we use centered finite difference approximations for the derivative terms at each interior node:\n",
+ "\n",
+ "\\\\[\\left. \\frac{d u_z}{dr} \\right|_{r_i} \\approx \\frac{u_{z,i+1} - u_{z,i-1}}{\\Delta r},\\\\]\n",
+ "\n",
+ "\\\\[\\left. \\frac{d^2 u_z}{dr^2} \\right|_{r_i} \\approx \\frac{u_{z,i+1} - 2 u_{z,i} + u_{z,i-1}}{\\Delta r^2}.\\\\]\n",
+ "\n",
+ "Thus, the equations for interior nodes become\n",
+ "\n",
+ "\\\\[\\frac{1}{r_i} \\frac{u_{z,i+1} - u_{z,i-1}}{\\Delta r}+\\frac{u_{z,i+1} - 2 u_{z,i} + u_{z,i-1}}{\\Delta r^2}=\\frac{-1000}{8.9\\times 10^{-4}}.\\\\]\n",
+ "\n",
+ "For the boundary nodes, at $r=R$ where the Dirchlet boundary condition is applied, we have\n",
+ "\n",
+ "\\\\[u_{z,n}=0.\\\\]\n",
+ "\n",
+ "At $r=0$ where the Neumann boundary condition is applied, a fictious node $u_{z,0}$ is introduced, thus we have\n",
+ "\n",
+ "\\\\[\\frac{u_{z,2}-{u_{z,0}}}{\\Delta r}=0.\\\\]\n",
+ "\n",
+ "Solving for the fictious node, we get $u_{z,0} = u_{z,2}$. The equation for the first node has to be modified because the radius is equal to zero at that node. We apply L’Hôpital’s rule:\n",
+ "\n",
+ "\\\\[\\lim_{r \\to 0} \\frac{1}{r} \\frac{d u_z}{dr} = \\lim_{r \\to 0} \\frac{\\frac{d}{dr}\\left(\\frac{du_z}{dr}\\right)}{\\frac{d}{dr}(r)} = \\lim_{r\\to 0}\\frac{d^2 u_z}{dr^2}.\\\\]\n",
+ "\n",
+ "So the equation at $r=0$ becomes $2\\frac{d^2 u_z}{dr^2} = \\frac{-1000}{8.9\\times 10^{-4}}$.\n",
+ "We apply centered finite differencing and plug the fictious node into the first equation to obtain\n",
+ "\n",
+ "\\\\[2\\cdot\\frac{2u_{z,2}-2u_{z,1}}{\\Delta r^2}=\\frac{-1000}{8.9\\times 10^{-4}}. \\\\]\n",
+ "\n",
+ "Now, we have a system of linear equations $\\mathbf A \\mathbf u_z =\\mathbf b$ and we can solve for the velocity profile."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "done\n"
+ ]
+ }
+ ],
+ "source": [
+ "# process specifications and constants\n",
+ "mu = 8.9e-4 # Pa*s\n",
+ "dPdz = -1000 # Pa/m\n",
+ "delta_r = 0.001 # m\n",
+ "r = range(0.0, 0.05, step = delta_r) # m\n",
+ "println(\"done\")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "