diff --git a/README.md b/README.md index 51d474b..493be28 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,29 @@ # Chemical Engineering - Analysis Notebooks Literate Programming Examples for Chemical Engineering Analysis -# Click this button to begin: [![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/mewilhel/AnalysisNotebooks/master) +## Authors + +## Copyright Information + +## Intended Use and Goals + +## Overview of Content + +This repository contains a series of supplemental notebooks used in XXX. + +The course made use of Dorfman and Daotidis's "Numerical Methods and Chemical Engineering Applications" as a primary textbook. A list of the chapters of Dorfman and Daotidis which correspond to each tutorial is provided below. We should however note that the enclosed material was design to be used in a stand alone fashion and may be readily integrated into other courses. + +- Tutorial \#0 (Introduction to Matlab/Julia) - N/A +- Tutorial \#1 (Linear Algebraic Systems) - Chapter 2 +- Tutorial \#2 (Nonlinear Algebraic Systems) - Chapter 3 +- Tutorial \#3 (ODEs, Initial Value Problems) - Chapter 4 +- Tutorial \#4 (ODEs, Boundary Value Problems) - Chapter 6 +- Tutorial \#5 (Partial Differential Equations) - Chapter 7 +- Tutorial \#6 (Optimization) - N/A + +The Ju + +## Extendability + +## References +- Dorfman, Kevin D., and Prodromos Daoutidis. Numerical Methods with Chemical Engineering Applications. Cambridge University Press, 2017. diff --git a/tutorial/julia/Tutorial 5 - Partial Differential Equations/Initial Value PDEs.ipynb b/tutorial/julia/Tutorial 5 - Partial Differential Equations/Initial Value PDEs.ipynb index 8e56613..bf9179b 100644 --- a/tutorial/julia/Tutorial 5 - Partial Differential Equations/Initial Value PDEs.ipynb +++ b/tutorial/julia/Tutorial 5 - Partial Differential Equations/Initial Value PDEs.ipynb @@ -6,8 +6,6 @@ "source": [ "# Initial Boundary Value PDEs \n", "\n", - "*Supplemental material (Julia) for Chapters 7 of \"Numerical Methods and Chemical Engineering Applications\" by Dorfman and Daoutidis*\n", - "\n", "
\n", "\n", "# Learning Objectives\n", @@ -27,59 +25,30 @@ "source": [ "# Simulating a graded catalyst bed\n", "\n", - "Process engineers are often engaged in tasks centered around increasing the economic productivity of process equipment while ensuring safe operation. For the commodity chemical industry, increased conversion of reactant to product remains a desired outcome. One of the most common used continuously operating reactor in this sector is that of the catalytic packed bed reactor. In this configuration, a tubular reactor is employed. This tube is packed with catalyst covered particles and the reactant flows over the catalyst bed. In many cases, a large number of tubes will be used in parallel cooled by the same jacket. We'll just consider the single tube variant here. A simplifed depiction is given below. \n", - "\n", - "\n", - "\n", - "\n", - "One important industrial chemical is that ortho-xylene which is used to produce phthalic anhydride a common plastizer. \n", - "The desired reaction is:\n", - "\n", - "\n", - "However, two underdesirable side reactions are also known to occur:\n", + "A suction pyrometer is a commonly used instrument for measuring gas phase temperatures in extemely hot environments (one example given [here](https://pubmed.ncbi.nlm.nih.gov/24248279/)). In these environments, it's quite typical to observe a large amount of radiative heat flux. This will typically dominate the heat transfer between the gas phase and an object. As a consequence, directly immersing a measurement device, such as a thermocouple or resistance temperature detector (RTD) into the process gas, will yield meaningless results. Moreover, an object left in this environment for a long period of time will increase in temperature until it too emits a significant amount of thermal radiation. So simply shielding a sensor isn't a viable solution either. This issue is circumvented by using a suction pyrometer. This device consistes of a small jacketing tube attached to a sensor (thermocouple or RTD) and a vacuum pump. By drawing gas through the tube, the jacket is cooled preventing emission, and allowing for a valid reading from the internal sensor. A picture is included below:\n", "\n", - "\n", + "\n", "\n", - "\n", + "The key design considerations behind the suction pyrometer are the superficial velocity and Biot number needed to ensure that the temperature rise of the process gas in the jacket is minimal. We'll construct a basic model of heat transfer in this tube and use it to determine adequate settings for reading our temperature range of interest. \n", "\n", - "The reactor is typically operated such at a hotspot occurs in the reactor, the reactor is stable (and in turn safe), and a high degree of selectivity and conversion is achieved. A number of process parameters including pressure, temperature, residence time, cooling rate, and feed composition may be manipulated to improve reactor productivity, selectivity, conversion, or overall profitability. We'll build a simulation (based on the model in this [paper](https://pubs.acs.org/doi/abs/10.1021/ie4005699).) to investigate the how using a graded catalyst may impact reactor performance. That is to say a reactor in which differing zones are loaded with catalyst of differing activity levels. This can allow the operator to raise the effective reactor temperature while preventing thermal runaway and maintaining adequate selectivity.\n", + "We'll consider a constant heat capacity gas which we're measuring using a 1ft (0.3048 m) by 0.25 inch (0.00635 m) diameter tube. Heat transfer at the wall surface is proportional to the temperature difference and characterized by the Biot number.\n", "\n", "## A basic descriptive model\n", "\n", - "A two-dimensional steady-state pseudohomogeneous packed bed reactor equations consist of a mass balance and an energy balance:\n", - "\n", - "\\\\[u_s\\frac{\\partial \\mathbf{c}}{\\partial z} = \\frac{D}{R}\\frac{\\partial}{\\partial R}\\left(R\\frac{\\partial \\mathbf{c}}{\\partial R}\\right) + \\sum_{i}{\\nu_{ij}r_i(\\mathbf{c},T)} \\\\]\n", - "\n", - "\\\\[u_s \\rho_{f} \\left(\\sum_{i}{C_{p,i}c_i}\\right) \\frac{\\partial T}{\\partial z} = \\frac{\\Lambda}{R}\\frac{\\partial}{\\partial R}\\left(R\\frac{\\partial T}{\\partial R}\\right) + \\sum_{i}{(-\\Delta H_{i})r_i(\\mathbf{c},T)} \\\\]\n", - "\n", - "were OX is denotes ortho-xylene, PA denotes phthalic anhydride and \n", - "\n", - "\\\\[j \\in \\{c_{OX,1}, c_{PA,1}, c_{H_2 O,1}, c_{O_2,1}, c_{CO_2,1}, c_{N_2}\\}\\\\].\n", - "\n", - "For a full review of different packed bed reactor models, the reader is directed to [this reference](https://ris.utwente.nl/ws/portalfiles/portal/6073612/t0000040.pdf).\n", + "A two-dimensional steady-state energy balance is given by:\n", "\n", - "
\n", - "Form of the above equations: The above equations are coupled convection-diffusion equations. One of the most commonly encountered forms in fluid mechanics.\n", - "
\n", - "\n", - "
\n", - "Limiting cases: Consider briefly under what conditions we might assume that these equations reduce to the Poisson or Laplace equations below.\n", - "
\n", + "\\\\[u_s C_{p} \\frac{\\partial T}{\\partial z} = \\frac{\\Lambda}{R}\\frac{\\partial}{\\partial R}\\left(R\\frac{\\partial T}{\\partial R}\\right), \\\\]\n", "\n", - "1. **Poisson Equation**: \\\\[\\nabla^2 x = f(x)\\\\]\n", - "2. **Laplace Equation**: \\\\[\\nabla^2 x = 0\\\\]\n", + "**Limiting cases:** Consider briefly under what conditions we might assume that this equation reduces to the Poisson or Laplace equations below.\n", "\n", - "Next, we assume a sufficient cooling flow exists to hold the reactor wall temperature constant at a fixed value of $T^c$. Then we can write a symmetry and cooling boundary condition:\n", + "1. **Poisson Equation:** \\\\[\\nabla^2 x = f(x) \\\\]\n", + "2. **Laplace Equation:** \\\\[\\nabla^2 x = 0 \\\\] \n", "\n", - "\\\\[\\left.\\frac{dT}{dR}\\right\\vert_{R=0} = 0 \\qquad \\qquad \\left.\\frac{dT}{dR}\\right\\vert_{R=R_t} = -Bi(T - T^c) \\\\]\n", - "\n", - "
\n", - "Note: The wall boundary condition ($R = R_t$) is of the Robin type and the symmetry condition ($R = 0$) is of the Neumann type.\n", - "
\n", + "Next, we assume a sufficient cooling flow exists to hold the reactor wall temperature constant at a fixed value of . Then, we can write a symmetry and cooling boundary condition:\n", "\n", - "We further assume that the reactor wall is impermeable to all species (and note that similar symmetry conditions should also hold):\n", + "\\\\[\\left.\\frac{\\partial T}{\\partial R}\\right\\vert_{R=0} = 0 \\qquad \\qquad \\left.\\frac{\\partial T}{\\partial R}\\right\\vert_{R=R_t} = -Bi(T - T^c) \\\\]\n", "\n", - "\\\\[\\left.\\frac{d\\mathbf{c}}{dR}\\right\\vert_{R=0} = 0 \\qquad \\qquad \\left.\\frac{d\\mathbf{c}}{dR}\\right\\vert_{R=R_t} = 0 \\\\]" + "Note that the symmetry condition (at $R = 0$) is of the Neumann type and the wall boundary condition (at $R = R_t$) is of the Robin type. " ] }, { @@ -88,43 +57,15 @@ "source": [ "## Further model development\n", "\n", - "- Assume the inlet stream is well-mixed.\n", - "\\\\[\\begin{align} \\mathbf{c}(z=0) &= \\mathbf{c}_{in} \\\\\n", - " T(z=0) &= T_{in} \\end{align}\\\\]\n", - " \n", - "
\n", - "Note: The inlet conditions are examples of Dirichlet boundary conditions.\n", - "
\n", - "\n", - "
\n", - "Activity! If the type of the boundary condition isn't immediate obvious take a moment and see if you can rearrange these equations into one of the below forms. \n", - "
\n", - "\n", - "1. **Dirichlet Boundary Condition**: \\\\[x = a\\\\] \n", - "2. **Neumann Boundary Condition**: \\\\[\\hat{n}\\cdot\\nabla x = f(x)\\\\]\n", - "3. **Robin Boundary Condition**: \\\\[\\hat{n}\\cdot\\nabla x + \\alpha x = f(x)\\\\]\n", + "Assume the inlet stream is well-mixed and at known temperature:\n", "\n", - "- All reactions are treated as pseudo-first-order. Each reaction rate may be written with relation to it's partial pressure which in turn is related to the concentration assuming an ideal gas equaton of state.\n", + "\\\\[\\mathbf T(z=0)=\\mathbf T_{in}. \\\\]\n", "\n", - "\\\\[r_1 = (\\sigma P_{O_2}\\rho_s) k_1 P_{OX} = k_1 (\\sigma \\rho_s)/(R_g T)^2 c_{O_2} c_{OX} \\\\]\n", - "\\\\[r_2 = (\\sigma P_{O_2}\\rho_s) k_2 P_{PA} = k_2 (\\sigma \\rho_s)/(R_g T)^2 c_{O_2} c_{PA} \\\\]\n", - "\\\\[r_3 = (\\sigma P_{O_2}\\rho_s) k_3 P_{OX} = k_3 (\\sigma \\rho_s)/(R_g T)^2 c_{O_2} c_{OX} \\\\]\n", - "\n", - "The rate constant can be calculated using the Arrhenius relationship as such\n", - "\n", - "\\\\[k_i = k_i^r \\exp{\\frac{E^r_i(T - T^r)}{T R_g T^r}}, \\qquad i = \\{1, 2, 3\\} \\\\]\n", - "\n", - "Individual component rates of change may then be written as:\n", - "\n", - "\\\\[\\begin{align} r_{OX} &= -r_1 + r_3 \\\\\n", - "r_{PA} &= r_1 - r_2 \\\\\n", - "r_{H_2 O} &= 3r_1 + 2r_2 + 5r_3 \\\\\n", - "r_{O_2} &= -3r_1 - 7.5r_2 + 10.5r_3 \\\\\n", - "r_{CO_2} &= 8r_2 + 8r_3 \\\\\n", - "r_{N_2} &= 0\n", - "\\end{align}\\\\]\n", - "\n", - "We omit a further discussion of model parameters in favor of including stating there values in the included code." + "This is because boundary conditions are known in the radial dimension and the system is fully specified at the reactor entrance for all variables (initial temperature).\n", + " \n", + "
\n", + "Note: This problem is treated as an IBVP and not a 2D BVP.\n", + "
" ] }, { @@ -146,17 +87,12 @@ "metadata": {}, "source": [ "### Method of Lines Derivation - Interior\n", - "Applying these rules to the mass balance for the interior nodes, we have \n", - "\n", - "\\\\[ \\left.\\frac{\\partial \\mathbf{c}}{\\partial z}\\right\\vert_{R=R_j} = u_s^{-1}\\left(\\frac{D}{R}\\frac{\\partial \\mathbf{c}}{\\partial R} + \\frac{\\partial^2 \\mathbf{c}}{\\partial R^2} + \\sum_{i}{\\nu_{ij}r_i(\\mathbf{c},T)}\\right) \\qquad j = 2, \\ldots, N_{r-1}\\\\]\n", - "\n", - "\\\\[ \\left.\\frac{d \\mathbf{c}_j}{d z}\\right\\vert_{R=R_j} = u_s^{-1}\\left(\\frac{D}{R}\\frac{\\mathbf{c}_{j+1} - \\mathbf{c}_{j-1}}{2\\Delta R} + D\\frac{\\mathbf{c}_{j+1} - 2\\mathbf{c}_{j} + \\mathbf{c}_{j-1}}{\\Delta R^2} + \\sum_{i}{\\nu_{ij}r_i(\\mathbf{c},T)}\\right) \\qquad j = 2, \\ldots, N_{r-1}\\\\]\n", + "Applying these rules to the energy balance for the interior nodes, where is the number of spatial nodes for discretizing the dimension, we have\n", "\n", - "Accordlingly, the temperature balance becomes\n", + "\\\\[\\frac{\\partial T}{\\partial z} = \\left(u_s C_P\\right)^{-1}\\left(\\frac{\\Lambda}{R}\\frac{\\partial}{\\partial R}\\left(R\\frac{\\partial T}{\\partial R}\\right)\\right) \\\\]\n", "\n", - "\\\\[\\left.\\frac{\\partial T}{\\partial z}\\right\\vert_{R=R_j} = \\left(u_s\\sum_{i}{f_{i}c_{pj}}\\right)^{-1}\\left(\\frac{\\Lambda}{R}\\frac{\\partial}{\\partial R}\\left(R\\frac{\\partial T}{\\partial R}\\right) + \\sum_{i}{(-\\Delta H_{i})r_i(\\mathbf{c},T)}\\right) \\qquad j = 2, \\ldots, N_{r-1} \\\\]\n", "\n", - "\\\\[ \\left.\\frac{d T}{d z}\\right\\vert_{R=R_j} = \\left(u_s\\sum_{i}{f_{i}c_{pj}}\\right)^{-1}\\left(\\frac{D}{R}\\frac{T_{j+1} - T_{j-1}}{2\\Delta R} + D\\frac{T_{j+1} - 2T_{j} + T_{j-1}}{\\Delta R^2} + \\sum_{i}{\\nu_{ij}r_i(\\mathbf{c},T)}\\right) \\qquad j = 2, \\ldots, N_{r-1}\\\\]" + "\\\\[\\left.\\frac{d T}{d z}\\right\\vert_{R=R_j} = \\left(u_s C_{p}\\right)^{-1}\\left(\\frac{D}{R}\\frac{T_{j+1} - T_{j-1}}{2\\Delta R} + D\\frac{T_{j+1} - 2T_{j} + T_{j-1}}{\\Delta R^2}\\right) \\qquad j = 2, \\ldots, N_{r-1} \\\\]" ] }, { @@ -164,20 +100,17 @@ "metadata": {}, "source": [ "### Method of Lines Derivation - Center\n", - "The discretization at the center node yeilds:\n", - "\\\\[ \\frac{d \\mathbf{c}_1}{d z} = u_s^{-1}\\left(\\frac{D}{R}\\frac{\\mathbf{c}_{2} - \\mathbf{c}_{0}}{2\\Delta R} + D\\frac{\\mathbf{c}_{2} - 2\\mathbf{c}_{1} + \\mathbf{c}_{0}}{\\Delta R^2} + \\sum_{i}{\\nu_{ij}r_i(c_1,T_1)}\\right)\\\\]\n", + "The discretization at the center node ($R=0, j=1$) yields: \n", "\n", - "from the boundary condition we have \\\\[\\frac{d \\mathbf{c}_1}{d R} = 0 \\rightarrow \\frac{\\mathbf{c}_{2} - \\mathbf{c}_{0}}{2\\Delta R} = 0 \\rightarrow \\mathbf{c}_{0} = \\mathbf{c}_{2}\\\\]. Substituting this in yeilds:\n", + "\\\\[\\frac{d T_1}{d z} = \\left(u_s C_P\\right)^{-1}\\left(\\frac{\\Lambda}{R}\\frac{T_{2} - T_{0}}{2\\Delta R} + \\Lambda \\frac{T_{2} - 2T_{1} + T_{0}}{\\Delta R^2}\\right).\\\\]\n", "\n", - "\\\\[ \\frac{d \\mathbf{c}_1}{d z} = u_s^{-1}\\left(\\frac{2D}{\\Delta R^2}(\\mathbf{c}_{2} - \\mathbf{c}_{1}) + \\sum_{i}{\\nu_{ij}r_i(\\mathbf{c}_1,T_1)}\\right)\\\\]\n", + "From the boundary condition we have \n", "\n", - "Accordingly for the energy balance we have:\n", + "\\\\[\\left.\\frac{\\partial T}{\\partial R}\\right|_{R=0} = 0 \\Rightarrow \\frac{T_{2} - T_{0}}{2\\Delta R} = 0 \\Rightarrow T_{0} = T_{2} \\\\]\n", "\n", - "\\\\[ \\frac{d T_1}{d z} = \\left(u_s\\sum_{i}{f_{i}c_{pj}}\\right)^{-1}\\left(\\frac{\\Lambda}{R}\\frac{T_{2} - T_{0}}{2\\Delta R} + \\Lambda \\frac{T_{2} - 2T_{1} + T_{0}}{\\Delta R^2} + \\sum_{i}{(-\\Delta H_{i})r_i(\\mathbf{c}_1,T_1)}\\right)\\\\]\n", + "Substituting this into the previous expression yields:\n", "\n", - "from the boundary condition we have \\\\[\\frac{d T_1}{d R} = 0 \\rightarrow \\frac{T_{2} - T_{0}}{2\\Delta R} = 0 \\rightarrow T_{0} = T_{2}\\\\]. Substituting this in yeilds:\n", - "\n", - "\\\\[ \\frac{d T_1}{d z} = \\left(u_s\\sum_{i}{f_{i}c_{pj}}\\right)^{-1}\\left(\\frac{2\\Lambda}{\\Delta R^2}(T_{2} - T_{1}) + \\sum_{i}{(-\\Delta H_{i})r_i(\\mathbf{c}_1,T_1)}\\right)\\\\]\n", + "\\\\[\\frac{d T_1}{d z} = \\left(u_s C_P\\right)^{-1}\\left(\\frac{2\\Lambda}{\\Delta R^2}(T_{2} - T_{1})\\right). \\\\]\n", "\n", "
\n", "Activity! This derivation proceeded by solving obtaining algebraic equations from the boundary condition, solving them analytically for the fictive node value, and then substituting this expression in. We could have instead solved a system of coupled differential and algebraic equations (referred to as a differential-algebraic system of equations or DAEs). In many cases, the algebraic equations formed may not have an analytic (closed form) solution and we must resort to solving DAEs. What type of boundary condition(s) may lead algebraic equations with no analytic solution?\n", @@ -190,639 +123,380 @@ "source": [ "### Method of Lines Derivation - Wall\n", "The discretization at the wall node yeilds:\n", - "\\\\[ \\frac{d \\mathbf{c}_R}{d z} = u_s^{-1}\\left(\\frac{D}{R}\\frac{\\mathbf{c}_{R+1} - \\mathbf{c}_{R-1}}{2\\Delta R} + D\\frac{\\mathbf{c}_{R+1} - 2\\mathbf{c}_{R} + \\mathbf{c}_{R-1}}{\\Delta R^2} + \\sum_{i}{\\nu_{ij}r_i(\\mathbf{c}_R,T_R)}\\right)\\\\]\n", - "\n", - "from the boundary condition we have \n", - "\n", - "\\\\[\\frac{d \\mathbf{c}_R}{d R} = 0 \\rightarrow \\frac{\\mathbf{c}_{R+1} - \\mathbf{c}_{R-1}}{2\\Delta R} = 0 \\rightarrow \\mathbf{c}_{R+1} = \\mathbf{c}_{R-1}.\\\\] \n", - "\n", - "Substituting this in yeilds:\n", - "\n", - "\\\\[ \\frac{d \\mathbf{c}_R}{d z} = u_s^{-1}\\left(\\frac{2D}{\\Delta R^2}(\\mathbf{c}_{R-1} - \\mathbf{c}_{R}) + \\sum_{i}{\\nu_{ij}r_i(\\mathbf{c}_R,T_R)}\\right)\\\\]\n", "\n", - "Accordingly for the energy balance we have:\n", + "\\\\[\\frac{d T_{N_r}}{d z} = \\left(u_s C_p\\right)^{-1}\\left(\\frac{\\Lambda}{R}\\frac{T_{N_{r+1}} - T_{N_{r-1}}}{2\\Delta R} + \\Lambda \\frac{T_{N_{r+1}} - 2T_{N_r} + T_{N_{r-1}}}{\\Delta R^2}\\right) \\\\]\n", "\n", - "\\\\[ \\frac{d T_R}{d z} = \\left(u_s\\sum_{i}{f_{i}\\mathbf{c}_{pj}}\\right)^{-1}\\left(\\frac{\\Lambda}{R}\\frac{T_{R+1} - T_{R-1}}{2\\Delta R} + \\Lambda \\frac{T_{R+1} - 2T_{R} + T_{R-1}}{\\Delta R^2} + \\sum_{i}{(-\\Delta H_{i})r_i(\\mathbf{c}_R,T_R)}\\right)\\\\]\n", + "From the boundary condition we have\n", "\n", - "from the boundary condition we have \n", + "\\\\[\\left.\\frac{\\partial T}{\\partial R}\\right|_{R=R_t} = 0 \\Rightarrow \\frac{T_{N_{r+1}} - T_{N_{r-1}}}{2\\Delta R} = -Bi(T_{N_r} - T^c) \\Rightarrow T_{N_{r+1}} = T_{N_{r-1}} -2Bi\\Delta R (T_{N_r} - T^c). \\\\]\n", "\n", - "\\\\[\\frac{d T_R}{d R} = 0 \\rightarrow \\frac{T_{R+1} - T_{R-1}}{2\\Delta R} = -Bi(T_R - T^c) \\rightarrow T_{R+1} = T_{R-1} -2Bi\\Delta R (T_R - T^c).\\\\]\n", - " \n", - "Substituting this in yeilds:\n", + "Substituting this in yields:\n", "\n", - "\\\\[ \\frac{d T_R}{d z} = \\left(u_s\\sum_{i}{f_{i}c_{pj}}\\right)^{-1}\\left(\\frac{Bi\\Lambda}{R}(T_{C} - T_{R}) + \\frac{2\\Lambda}{\\Delta R^2}(T_{R-1} - T_{R} + Bi\\Delta R (T_{C} - T_{R})) + \\sum_{i}{(-\\Delta H_{i})r_i(\\mathbf{c}_R,T_R)}\\right)\\\\]\n", + "\\\\[\\frac{d T_{N_r}}{d z} = \\left(u_s C_p\\right)^{-1}\\left(\\frac{Bi\\Lambda}{R}(T^c - T_{N_r}) + \\frac{2\\Lambda}{\\Delta R^2}(T_{N_{r-1}} - T_{N_r} + Bi\\Delta R (T^c - T_{N_r}))\\right). \\\\]\n", "\n", - "We now need to interlace the variables letting \\\\[ y = \\{T_1, c_{OX,1}, c_{PA,1}, c_{H_2 O,1}, c_{O_2,1}, c_{CO_2,1}, c_{N_2}, \\ldots, T_{N_R}, c_{OX,N_R}, c_{PA,N_R}, c_{H_2 O,N_R}, c_{O_2,N_R}, c_{CO_2,N_R}, c_{N_2,N_R}\\}\\\\]\n", - "\n", - "We'll now use a nonlinear equation solving package in Julia to write a simple a fixed-step size implicit Euler integration method (progressively, solving the below equation for $y_{z+1}$)\n", - "\n", - "\\\\[ y_{z+1} - y_z - \\Delta z f(z_{k+1}) = 0 \\\\]\n", + "We'll now use a ODE solver readily available in DifferentialEquations.jl to integrate the above system of differential equations.\n", "\n", "
\n", - "Activity! Take some time to look over the code for the residual. Where does the right-hand side of the ODE system appear? Note how the physical property calculations are implemented.\n", - "
\n" + "Activity! Manipulate the code below, to determine the superficial velocity at which a temperature rise of 2K is observed in the tube for gas temperatures ranging from 900 to 1000.\n", + "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "
" + "## Solution\n", + "\n", + "For simplicity's sake, we'll assume the following physical constants are valid over the range of interest (constant heat capacities, and thermal heat dispersion coefficients)." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "# Input model parameters\n", + "tube_length = 0.3048 # tube length [m]\n", + "tube_radius = 0.00635/2 # radius [m]\n", + "\n", + "heat_coeff = 7.3871/3.6 # heat dispersion coefficient [J/(m s K)]\n", + "\n", + "Tc = 1400 # Tube temperature [K]\n", + "T0 = 900 # Gas design temperature [K] MANIPULATE THIS \n", + "biot_number = 0.001 # Biot number, Bi \n", + "\n", + "us = 0.1 # superficial flow rate (m/h) MANIPULATE THIS \n", + "Cp = 0.74 # pseudo heat capacity\n", + "\n", + "nr = 20 # number of discretization points in radial direction\n", + "delR = tube_radius/(nr - 1); " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "For simplicity sake, we'll assume the follow physical constants are valid over the range of interest (constant heat capacities, enthalpy of reaction, diffusion, and thermal heat dispersion coefficients but we'll allow reaction rate to vary with temperature." + "We'll use the adaptive step-size implicit Euler method built into DifferentialEquations.jl to solve the problem. You can read more about the built-in solvers and functionality for DifferentialEquations.jl [here](https://diffeq.sciml.ai/v2.0/). We first define the right-hand side function for the differential equation below. " ] }, { "cell_type": "code", - "execution_count": 1, + "execution_count": 7, "metadata": {}, "outputs": [ { - "name": "stdout", - "output_type": "stream", - "text": [ - "residual function defined\n" - ] + "data": { + "text/plain": [ + "f (generic function with 1 method)" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" } ], "source": [ - "#=\n", - "Define rhs of implicit Euler step.\n", - "=#\n", - "function residual!(out, y1, y2::Vector{W}, a, params) where W\n", - "\n", - " # unpack parameters\n", - " nr, nz, delR, delZ, catalyst_density, reactor_diam, biot_number,\n", - " Cp, delH, Tc, diff_coeff, heat_coeff, Rg, k_ref, E_ref, T_ref, us = params\n", - "\n", - " # compute step size\n", - " dr = 1/(nr - 1)\n", + "function f(out, y, p, t)\n", + " # Set up center boundary condition rhs equations\n", + " out[1] = (Cp*us)^(-1)*((2*heat_coeff/delR^2)*(y[3] - y[1]))\n", "\n", - " # preallocate storage\n", - " r_rxn = zeros(W,3)\n", - " r_prod = zeros(W,6)\n", - " k = zeros(W,3)\n", - "\n", - " # pollulate the rhs of the ODEs for the interior nodes\n", + " # Define interior node right-hand-side equations\n", " for i = 2:(nr-1)\n", - "\n", - " # get node radial position and temperature\n", - " r = (i-1)*nr\n", - " T = y2[7*(i-1) + 1]\n", - "\n", - " # compute reaction rate constant\n", - " k[1] = k_ref[1]*exp(E_ref[1]*(T - T_ref)/(T*Rg*T_ref))\n", - " k[2] = k_ref[2]*exp(E_ref[2]*(T - T_ref)/(T*Rg*T_ref))\n", - " k[3] = k_ref[3]*exp(E_ref[3]*(T - T_ref)/(T*Rg*T_ref))\n", - "\n", - " # compute reaction rate\n", - " r_rxn[1] = a*catalyst_density*k[1]*y2[7*(i-1) + 2]*y2[7*(i-1) + 5]\n", - " r_rxn[2] = a*catalyst_density*k[2]*y2[7*(i-1) + 3]*y2[7*(i-1) + 5]\n", - " r_rxn[3] = a*catalyst_density*k[3]*y2[7*(i-1) + 2]*y2[7*(i-1) + 5]\n", - "\n", - " # compute production rate\n", - " r_prod[1] = -r_rxn[1] - r_rxn[3]\n", - " r_prod[2] = r_rxn[1] - r_rxn[2]\n", - " r_prod[3] = 3.0*r_rxn[1] + 2.0*r_rxn[2] + 5.0*r_rxn[3]\n", - " r_prod[4] = -3.0*r_rxn[1] - 7.5*r_rxn[2] - 10.5*r_rxn[3]\n", - " r_prod[5] = 8.0*r_rxn[2] + 8.0*r_rxn[3]\n", - " r_prod[6] = 0.0\n", - "\n", - " # energy balance discretized equation\n", - " prefT = delZ*Cp*us\n", - " sourceT = zero(T)\n", - " for j = 1:3\n", - " sourceT += (-delH[j])*r_rxn[j]\n", - " end\n", - " Tp1 = y2[7*(i-1) + 8]\n", - " Tm1 = y2[7*(i-1) - 6]\n", - " term1 = heat_coeff*(Tp1 - Tm1)/(2*dr*r)\n", - " term2 = heat_coeff*(Tp1 - 2*T + Tm1)/dr^2\n", - " out[7*(i-1) + 1] = y2[7*(i-1) + 1] - y1[7*(i-1) + 1] + prefT*(term1 + term2 + sourceT)\n", - "\n", - " # mass balance discretized equations\n", - " prefC = delZ*us\n", - " for j = 0:5\n", - " Cp1 = y2[7*(i-1) + 9 + j]\n", - " C1 = y2[7*(i-1) + 2 + j]\n", - " Cm1 = y2[7*(i-1) - 5 + j]\n", - " term1 = diff_coeff*(Tp1 - Tm1)/(2*dr*r)\n", - " term2 = diff_coeff*(Tp1 - 2*T + Tm1)/dr^2\n", - " out[7*(i-1) + 2 + j] = y2[7*(i-1) + 2 + j] - y1[7*(i-1) + 2 + j] + prefC*(term1 + term2 + r_prod[6 - j])\n", - " end\n", - " end\n", - "\n", - " # pollulate the rhs of the ODEs for the center nodes\n", - " T = y2[1]\n", - "\n", - " # compute reaction rate constant\n", - " k[1] = k_ref[1]*exp(E_ref[1]*(T - T_ref)/(T*Rg*T_ref))\n", - " k[2] = k_ref[2]*exp(E_ref[2]*(T - T_ref)/(T*Rg*T_ref))\n", - " k[3] = k_ref[3]*exp(E_ref[3]*(T - T_ref)/(T*Rg*T_ref))\n", - "\n", - " # compute reaction rate\n", - " r_rxn[1] = a*catalyst_density*k[1]*y2[2]*y2[5]\n", - " r_rxn[2] = a*catalyst_density*k[2]*y2[3]*y2[5]\n", - " r_rxn[3] = a*catalyst_density*k[3]*y2[2]*y2[5]\n", - "\n", - " # compute production rate\n", - " r_prod[1] = -r_rxn[1] - r_rxn[3]\n", - " r_prod[2] = r_rxn[1] - r_rxn[2]\n", - " r_prod[3] = 3.0*r_rxn[1] + 2.0*r_rxn[2] + 5.0*r_rxn[3]\n", - " r_prod[4] = -3.0*r_rxn[1] - 7.5*r_rxn[2] - 10.5*r_rxn[3]\n", - " r_prod[5] = 8.0*r_rxn[2] + 8.0*r_rxn[3]\n", - " r_prod[6] = 0.0\n", - "\n", - " # energy balance discretized equation\n", - " prefT = delZ*Cp*us\n", - " sourceT = zero(T)\n", - " for i = 1:3\n", - " sourceT += (-delH[i])*r_rxn[i]\n", - " end\n", - " out[1] = y2[1] - y1[1] + prefT*((2*heat_coeff/dr^2)*(y2[8] - T) + sourceT)\n", - "\n", - " # mass balance discretized equations\n", - " prefC = delZ*us\n", - " for i = 0:5\n", - " CR1 = y2[i + 9]\n", - " CR = y2[i + 2]\n", - " out[i + 2] = y2[i + 2] - y1[i + 2] + prefC*((2*diff_coeff/dr^2)*(CR1 - CR) + r_prod[6 - i])\n", + " r = i*delR; # get node radial position\n", + " # energy balance discretized equation \n", + " out[i] = (heat_coeff*(Cp*us)^-1)*((y[i+1] - y[i])/(2*delR*r) + (y[i+1] - 2*y[i] + y[i-1])/delR^2)\n", " end\n", "\n", - " # pollulate the rhs of the ODEs for the walls nodes\n", - " T = y2[7*nr - 6]\n", - "\n", - " # compute reaction rate constant\n", - " k[1] = k_ref[1]*exp(E_ref[1]*(T - T_ref)/(T*Rg*T_ref))\n", - " k[2] = k_ref[2]*exp(E_ref[2]*(T - T_ref)/(T*Rg*T_ref))\n", - " k[3] = k_ref[3]*exp(E_ref[3]*(T - T_ref)/(T*Rg*T_ref))\n", - "\n", - " # compute reaction rate\n", - " r_rxn[1] = a*catalyst_density*k[1]*y2[1]*y2[4]\n", - " r_rxn[2] = a*catalyst_density*k[2]*y2[2]*y2[4]\n", - " r_rxn[3] = a*catalyst_density*k[3]*y2[1]*y2[4]\n", - "\n", - " # compute production rate\n", - " r_prod[1] = -r_rxn[1] - r_rxn[3]\n", - " r_prod[2] = r_rxn[1] - r_rxn[2]\n", - " r_prod[3] = 3.0*r_rxn[1] + 2.0*r_rxn[2] + 5.0*r_rxn[3]\n", - " r_prod[4] = -3.0*r_rxn[1] - 7.5*r_rxn[2] - 10.5*r_rxn[3]\n", - " r_prod[5] = 8.0*r_rxn[2] + 8.0*r_rxn[3]\n", - " r_prod[6] = 0.0\n", - "\n", - " # energy balance discretized equation\n", - " prefT = delZ*Cp*us\n", - " sourceT = zero(T)\n", - " for i = 1:3\n", - " sourceT += (-delH[i])*r_rxn[i]\n", - " end\n", - " TR_m_1 = y2[7*nr - 13]\n", - " term1 = (biot_number*heat_coeff/reactor_diam)*(Tc - T)\n", - " term2 = (2*heat_coeff/dr^2)*(TR_m_1 - T + biot_number*dr*(Tc - T))\n", - " out[7*nr - 6] = y2[7*nr - 6] - y1[7*nr - 6] + prefT*(term1 + term2 + sourceT)\n", - "\n", - " # mass balance discretized equations\n", - " prefC = delZ*us\n", - " for i = 0:5\n", - " CR1 = y2[7*nr - i]\n", - " CR = y2[7*(nr-1) - i]\n", - " out[7*nr - i] = y2[7*nr - i] - y1[7*nr - i] + prefC*((2*diff_coeff/dr^2)*(CR1 - CR) + r_prod[6 - i])\n", - " end\n", - "\n", - " return nothing\n", - "end\n", - "\n", - "println(\"residual function defined\")" + " # Set up wall boundary condition rhs equations\n", + " r = (nr-1)*delR\n", + " out[nr] = (2*heat_coeff*(Cp*us)^(-1))*(y[nr - 1]/delR^2 - y[nr]*(biot_number*(1/r + 1/delR) + 1/delR^2) + \n", + " biot_number*(1/(2*r) + 1/delR)*Tc)\n", + " return \n", + "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "We'll use an off the shelf implementation of Newton's method through the package NLSolve.jl. " + "Now we build and solve the corresponding ODEs problem. \n", + "\n", + "
\n", + "Note: If this is your first time using DifferentialEquations.jl or Plots.jl some time may be spent precompiling. So taking a coffee break after clicking run on the below cell wouldn't be ill-advised.\n", + "
" ] }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ - "\u001b[32m\u001b[1m Updating\u001b[22m\u001b[39m registry at `C:\\Users\\wilhe\\.julia\\registries\\General`\n", "\u001b[32m\u001b[1m Resolving\u001b[22m\u001b[39m package versions...\n", "\u001b[32m\u001b[1mNo Changes\u001b[22m\u001b[39m to `C:\\Users\\wilhe\\Project.toml`\n", - "\u001b[32m\u001b[1mNo Changes\u001b[22m\u001b[39m to `C:\\Users\\wilhe\\Manifest.toml`\n" + "\u001b[32m\u001b[1mNo Changes\u001b[22m\u001b[39m to `C:\\Users\\wilhe\\Manifest.toml`\n", + "\u001b[32m\u001b[1m Resolving\u001b[22m\u001b[39m package versions..." ] }, { "name": "stdout", "output_type": "stream", "text": [ - "defined integration scheme\n" + "ODE Integrated" ] - } - ], - "source": [ - "import Pkg; Pkg.add(\"NLsolve\")\n", - "using NLsolve\n", - "\n", - "function activity_profile(v)\n", - " if v < 0.25\n", - " return 0.9 # keep between 0.9 and 1.0\n", - " elseif v < 1.0\n", - " return 0.7\n", - " elseif v < 1.5\n", - " return 1.0\n", - " end\n", - " return 1.0\n", - "end\n", - "\n", - "function nonlinear_step_solve!(out, indx, z, initial_y, params)\n", - " a = activity_profile(z)\n", - " residual_func! = (out, y2) -> residual!(out, initial_y, y2, a, params)\n", - " result = nlsolve(residual_func!, initial_y, autodiff = :forward, method = :newton)\n", - " out[indx,:] = result.zero\n", - " return nothing\n", - "end\n", - "\n", - "function integrate_ODEs!(out, initial_y, params)\n", - " nr, nz, delR, delZ, catalyst_density, reactor_diam, biot_number,\n", - " Cp, delH, Tc, diff_coeff, heat_coeff, Rg, k_ref, E_ref, T_ref, us = params\n", - " z = 0.0\n", - " out[1,:] = initial_y # store initial condition\n", - " for indx = 2:nz\n", - " z += delZ # advance axial dimension value\n", - " nonlinear_step_solve!(out, indx, z, out[indx-1,:], params) # perform a nonlinear solve for next step\n", - " end\n", - " return nothing\n", - "end\n", - "\n", - "println(\"defined integration scheme\")" - ] - }, - { - "cell_type": "code", - "execution_count": 3, - "metadata": {}, - "outputs": [ + }, { - "name": "stdout", + "name": "stderr", "output_type": "stream", "text": [ - "parameter defined\n", - "ODE integrated\n" + "\n", + "\u001b[32m\u001b[1mNo Changes\u001b[22m\u001b[39m to `C:\\Users\\wilhe\\Project.toml`\n", + "\u001b[32m\u001b[1mNo Changes\u001b[22m\u001b[39m to `C:\\Users\\wilhe\\Manifest.toml`\n" ] - } - ], - "source": [ - "# set up discretization points\n", - "nr = 20 # number of discretization points in radial direction\n", - "nz = 80 # number of discretization points in axial direction\n", - "\n", - "# set reactor design parameters\n", - "reactor_length = 3 # length of reactor \n", - "reactor_diam = 0.0254 # diameter of reactor \n", - "\n", - "delR = reactor_diam/(nr - 1)\n", - "delZ = reactor_length/(nz - 1)\n", - "\n", - "# Input model parameters\n", - "k_ref = [6.519E-2; 5.698E-3; 6.442E-3] # Reaction rate at refence temperature\n", - "E_ref = [113.57; 129.71; 119.68] # Activation energy at reference temperature\n", - "T_ref = 600 # Reference temperature\n", - "Rg = 8.3144*1000 # ideal gas constant\n", - "\n", - "catalyst_density = 1300 # Density of packed catalytic material\n", - "diff_coeff = 0.5\n", - "heat_coeff = 7.3871 # heat dispersion coefficient\n", - "\n", - "Tc = 323.15 # 50C water used for cooling\n", - "T0 = 625\n", - "biot_number = 1.6 # Biot number, Bi\n", - "\n", - "us = 0.8\n", - "Cp = 0.992 # approximate\n", - "delH = [-1.361; -14.329; -18.809]*1000 # heat of reaction\n", - "\n", - "# define initial condition\n", - "feed_rate = 0.174 # total molar feed rate\n", - "feed_conc_split = [0.011; 0.208; 0.781; 0.0; 0.0; 0.0] # mole fraction\n", - "\n", - "initial_y = zeros(7*nr) # create vector of zeros\n", - "for i = 1:nr\n", - " initial_y[7*(i-1) + 1] = T0\n", - " initial_y[7*(i-1) + 2] = feed_conc_split[1]*feed_rate\n", - " initial_y[7*(i-1) + 3] = feed_conc_split[2]*feed_rate\n", - " initial_y[7*(i-1) + 4] = feed_conc_split[3]*feed_rate # only assign nonzero values\n", - "end\n", - "\n", - "# specify parameters\n", - "params = nr, nz, delR, delZ, catalyst_density, reactor_diam, biot_number,\n", - "Cp, delH, Tc, diff_coeff, heat_coeff, Rg, k_ref, E_ref, T_ref, us\n", - "\n", - "println(\"parameter defined\")\n", - "\n", - "# preallocate storage\n", - "out = zeros(nz, 7*nr)\n", - "\n", - "# solve the system\n", - "integrate_ODEs!(out, initial_y, params)\n", - "\n", - "println(\"ODE integrated\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Extract coupled values and generate 2D plot" - ] - }, - { - "cell_type": "code", - "execution_count": 4, - "metadata": {}, - "outputs": [ + }, { "name": "stdout", "output_type": "stream", "text": [ - "properties extracted\n" + "\n" ] } ], "source": [ - "function extract_values(out, nr, nz)\n", - " T = zeros(nr,nz)\n", - " cOX = zeros(nr,nz)\n", - " cPA = zeros(nr,nz)\n", - " cH2O = zeros(nr,nz)\n", - " cO2 = zeros(nr,nz)\n", - " cCO2 = zeros(nr,nz)\n", - " cN2 = zeros(nr,nz)\n", - " for i = 1:nr\n", - " for j = 1:nz\n", - " T[i,j] = out[j, 7*(i-1) + 1]\n", - " cOX[i,j] = out[j, 7*(i-1) + 2]\n", - " cPA[i,j] = out[j, 7*(i-1) + 3]\n", - " cH2O[i,j] = out[j, 7*(i-1) + 4]\n", - " cO2[i,j] = out[j, 7*(i-1) + 5]\n", - " cCO2[i,j] = out[j, 7*(i-1) + 6]\n", - " cN2[i,j] = out[j, 7*(i-1) + 7]\n", - " end\n", - " end\n", - " return T, cOX, cPA, cH2O, cO2, cCO2, cN2\n", - "end\n", - "\n", - "T, cOX, cPA, cH2O, cO2, cCO2, cN2 = extract_values(out, nr, nz)\n", + "import Pkg; Pkg.add(\"DifferentialEquations\"); Pkg.add(\"Plots\")\n", + "using DifferentialEquations\n", "\n", - "println(\"properties extracted\")" + "y0 = [T0 for i = 1:nr] # set initial condition to inlet temperature\n", + "prob = ODEProblem(f, y0, (0, tube_length)) # create an ODE problem and solve it\n", + "sol = solve(prob)\n", + "println(\"ODE Integrated\") # Displays when complete" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "We now generate a plot of the centerline temperature versus position. If Plots.jl hasn't been installed yet this may take some time." + "We now generate a plot of the centerline temperature versus position. If Plots.jl hasn't been used before yet this may also take some time." ] }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 15, "metadata": {}, "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "\u001b[32m\u001b[1m Resolving\u001b[22m\u001b[39m package versions...\n", - "\u001b[32m\u001b[1mNo Changes\u001b[22m\u001b[39m to `C:\\Users\\wilhe\\Project.toml`\n", - "\u001b[32m\u001b[1mNo Changes\u001b[22m\u001b[39m to `C:\\Users\\wilhe\\Manifest.toml`\n" - ] - }, { "data": { "image/svg+xml": [ "\n", "\n", "\n", - " \n", + " \n", " \n", " \n", "\n", - "\n", "\n", - " \n", + " \n", " \n", " \n", "\n", - "\n", "\n", - " \n", - " \n", + " \n", + " \n", " \n", "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n" + "\n", + "\n", + "\n", + "\n" ] }, - "execution_count": 9, + "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "import Pkg; Pkg.add(\"Plots\")\n", "using Plots\n", - "\n", - "Tcenter = [T[1,i] for i in 1:nz]\n", - "zpnt = [i for i in range(0.0, reactor_length, length = nz)]\n", - "plot(zpnt, Tcenter, label=\"\", linealpha = 0.5, marker = (:hexagon, 2, 0.6, :green))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "
\n", - "Activity! In the below box, compute the average temperature, concentrations, and then conversion at each spatial position. Next try, to modifying the activity_profile and re-running the example. Can you pick out can clear trend in the behavior of the system with respect to temperature/conversion with respect to activity profile?\n", - "
" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "If you manipulated this example substantially, you may have discovered that some of the simulation results you reached were not physically meaningful. Large oscillations may be observed or even negative values. This is actually expected behavior for many PDE systems which may require specialized discretization methods to be dealt with (such as those encountered in [this reference]((https://ris.utwente.nl/ws/portalfiles/portal/6073612/t0000040.pdf))). Alternatively, one can formulate this simulation as an optimization problem (which can accept constraints such as nonnegative concentration). This later approach can be of particular interest when we are attempting to choose design/process parameters from a broad set of possibility in which many options are invalid." + "plot(sol, label=\"Centerline Temperature\", legend=:bottomright, vars = 1)" ] }, { diff --git a/tutorial/julia/Tutorial 5 - Partial Differential Equations/pyrometer_pde_pic.png b/tutorial/julia/Tutorial 5 - Partial Differential Equations/pyrometer_pde_pic.png new file mode 100644 index 0000000..3b8f136 Binary files /dev/null and b/tutorial/julia/Tutorial 5 - Partial Differential Equations/pyrometer_pde_pic.png differ diff --git a/tutorial/julia/Tutorial 5 - Partial Differential Equations/react1edited.png b/tutorial/julia/Tutorial 5 - Partial Differential Equations/react1edited.png deleted file mode 100644 index 04bd482..0000000 Binary files a/tutorial/julia/Tutorial 5 - Partial Differential Equations/react1edited.png and /dev/null differ diff --git a/tutorial/julia/Tutorial 5 - Partial Differential Equations/react2edited.png b/tutorial/julia/Tutorial 5 - Partial Differential Equations/react2edited.png deleted file mode 100644 index e489b5d..0000000 Binary files a/tutorial/julia/Tutorial 5 - Partial Differential Equations/react2edited.png and /dev/null differ diff --git a/tutorial/julia/Tutorial 5 - Partial Differential Equations/react3edited.png b/tutorial/julia/Tutorial 5 - Partial Differential Equations/react3edited.png deleted file mode 100644 index 8a264dc..0000000 Binary files a/tutorial/julia/Tutorial 5 - Partial Differential Equations/react3edited.png and /dev/null differ diff --git a/tutorial/matlab/Tutorial 5 - Partial Differential Equations/M5_PDEsystems.mlx b/tutorial/matlab/Tutorial 5 - Partial Differential Equations/M5_PDEsystems.mlx index 482556f..4bc1301 100644 Binary files a/tutorial/matlab/Tutorial 5 - Partial Differential Equations/M5_PDEsystems.mlx and b/tutorial/matlab/Tutorial 5 - Partial Differential Equations/M5_PDEsystems.mlx differ