MediaWiki API result

{
    "batchcomplete": "",
    "continue": {
        "gapcontinue": "Suggested_readings",
        "continue": "gapcontinue||"
    },
    "warnings": {
        "main": {
            "*": "Subscribe to the mediawiki-api-announce mailing list at <https://lists.wikimedia.org/postorius/lists/mediawiki-api-announce.lists.wikimedia.org/> for notice of API deprecations and breaking changes."
        },
        "revisions": {
            "*": "Because \"rvslots\" was not specified, a legacy format has been used for the output. This format is deprecated, and in the future the new format will always be used."
        }
    },
    "query": {
        "pages": {
            "64": {
                "pageid": 64,
                "ns": 0,
                "title": "Software for integrating ODEs",
                "revisions": [
                    {
                        "contentformat": "text/x-wiki",
                        "contentmodel": "wikitext",
                        "*": "== Background ==\n\nNumerically integrating ODE's is something you should have mastered in your pre-requisite 3E4 course. If you are not familiar with this topic, it is your responsibility to quickly catch up, because we will use this intensively for the rest of the course.\n\nHere is a tutorial a wrote a few years ago when I [http://modelling3e4.connectmv.com/wiki/Software_tutorial/Integration_of_ODEs taught 3E4 in 2010]. The tutorial below is similar, but uses a reactor design example.\n\n== Example ==\n\nConsider the pair of differential equations (covered in class on [[Isothermal_reactor_design_-_2013|11 February]]):\n<rst>\n<rst-options: 'toc' = False/>\n<rst-options: 'reset-figures' = False/>\n.. math::\n\n\t\\dfrac{dX}{dW} &= \\dfrac{-r_A'}{F_{A0}} \\\\ \t\t\\dfrac{dy}{dW} &= -\\dfrac{\\alpha}{2y}\\left(1 + \\varepsilon X\\right)\n\t\nSome terminology (recap from your pre-requisite math courses)\n\n\t*\tThe independent variable is \\(W\\)\n\t*\tThe two dependent variables (variables being integrated with respect to the independent variable) are \\(X\\) and \\(y\\)\n\t\nSince there are two dependent variables, we require initial conditions for each variable. In this case:\n\n\t*\t:math:`X(0) = 0.0`\n\t*\t:math:`y(0) = 1.0`\n\t\nWe also need to specify initial **and final conditions** for the independent variable:\n\n\t*\t:math:`W` initially is 0.0 at the reactor entrance\n\t* \t:math:`W` finally is 28.0 kg at the reactor exit\n\t\nBut there are a few other unknowns we must first specify:\n\n\t*\t:math:`r_A' = -k' \\dfrac{(1-X)y}{(1+\\varepsilon X)}`\n\t*\t:math:`q = \\dfrac{(1 + \\varepsilon X)}{y}`\n\t*\t:math:`F_{A0} = 0.1362\\,\\text{mol.s}^{-1}`\n\t*\t:math:`k' = 0.0074\\,\\text{mol.s}^{-1}\\text{.(kg catalyst)}^{-1}`\n\t*\t:math:`\\alpha = 0.0367\\,\\text{kg}^{-1}`\n\t*\t:math:`\\varepsilon = -0.15`\n</rst>\n\n== Polymath==\n\nDownload and install Polymath from the CD/DVD included with your course textbook (Windows only; sorry Mac and Linux versions are not available). Unfortunately this is only a 15 day version, which is pretty useless for this course, since you will require the code for at least the next 30 to 45 days. So rather use one of the other options below, if possible.\n\nEnter the following Polymath code (as demonstrated in class on [[Isothermal_reactor_design_-_2013|11 February]]. Check the video out for a detailed explanation of the code.\n<syntaxhighlight lang=\"text\">\n# Differential equations\nd(X) / d(W) = -rAdash / FA0 \nX(0) = 0\nd(y) / d(W) = -alpha/(2*y) * (1+eps*X) \ny(0) = 1.0\n\n# Constants\nFA0   = 0.1362  # [mol/s]\nkdash = 0.0074  # [mol/(kg catalyst . s)]\nalpha = 0.0367  # [1/kg]\neps   = -0.15   # [-]\n\n# Algebraic equations\nrAdash = -kdash * (1-X)/(1+eps*X) * y\nflow_ratio = (1 + eps*X)/y\n\n# Initial and final values for independent variable:\nW(0) = 0\nW(f) = 28\n</syntaxhighlight>\nEnsure that you obtain a graphical output as shown here. This represents the profile of conversion, \\(X\\) and the pressure drop ratio \\(y\\) throughout the reactor.\n[[Image:Polymath-output.jpg]]\n\n== Python ==\n=== Installing Python ===\nInstallation differs from Windows to Linux to Mac. One Python distribution that works well for many people is the [https://www.enthought.com/ Enthought Python Distribution] (EPD) and the other is [https://www.continuum.io/downloads Anaconda]. Another version, especially for Windows users is [https://python-xy.github.io/ Python(x,y)].\n\n=== Libraries required ===\nYou are ready to go if you enter the following without getting error messages:\n<syntaxhighlight lang=\"python\">\nimport numpy as np\nfrom scipy import integrate\nfrom matplotlib.pylab import (plot, grid, xlabel, ylabel)\n</syntaxhighlight>\n\nIf you experience problems importing these libraries, then re-install your Python distribution using one of the two options listed above.\n\n=== Source code ===\n<syntaxhighlight lang=\"python\">\nimport numpy as np\nfrom scipy import integrate\nfrom matplotlib.pylab import (plot, grid, xlabel, ylabel, show, legend)\n \ndef PBR_model(indep, depnt):\n    \"\"\"\n    Dynamic balance for the packed bed reactor (PBR); demo problem class 05C\n \n    indep: the independent ODE variable, such as time or length\n    depnt: a vector of dependent variables\n \n    X = depnt[0] = the conversion\n    y = depnt[1] = the pressure ratio = P/P_0 = y\n \n    Returns d(depnt)/d(indep) = a vector of ODEs\n    \"\"\"\n \n    # Assign some variables for convenience of notation\n    X = depnt[0]\n    y = depnt[1]\n \n    # Constants\n    FA0   = 0.1362  # [mol/s]\n    kdash = 0.0074  # [mol/(kg catalyst . s)]\n    alpha = 0.0367  # [1/kg]\n    eps   = -0.15   # [-]\n \n    #Algebraic equations\n    rAdash = -kdash * (1-X)/(1+eps*X) * y\n \n    # Output from this ODE function must be a COLUMN vector, with n rows\n    n = len(depnt) \n    d_depnt__d_indep = np.zeros((n,1))\n    d_depnt__d_indep[0] = -rAdash / FA0 \n    d_depnt__d_indep[1] = -alpha/(2*y) * (1+eps*X) \n    return d_depnt__d_indep\n \n# The \"driver\" that will integrate the ODE(s):\n# -----------\n# Start by specifying the integrator:\n# use ``vode`` with \"backward differentiation formula\"\nr = integrate.ode(PBR_model).set_integrator('vode', method='bdf')\n \n# Set the independent variable's range\nindep_start = 0.0   # kg\nindep_final = 28.0  # kg\ndelta = 1.0         # the results will be shown only at these ``delta`` points \nnum_steps = np.floor((indep_final - indep_start)/delta) + 1 # Number of steps: 1 extra for initial condition\n \n# Set initial condition(s): for integrating variables (dependent variables)\nX_depnt_zero = 0.0   # i.e. X(W=0) = 0.0\ny_depnt_zero = 1.0   # i.e. y(W=0) = 1.0\nr.set_initial_value([X_depnt_zero, y_depnt_zero], indep_start)\n \n# Create vectors to store trajectories\nW = np.zeros((num_steps, 1))\nX = np.zeros((num_steps, 1))\ny = np.zeros((num_steps, 1))\nW[0] = indep_start\nX[0] = X_depnt_zero\ny[0] = y_depnt_zero\n \n# Integrate the ODE(s) across each delta\nk = 1\nwhile r.successful() and k < num_steps:\n    r.integrate(r.t + delta)\n \n    # Store the results to plot later\n    W[k] = r.t\n    X[k] = r.y[0]\n    y[k] = r.y[1]\n    k += 1\n \n# All done!  Plot the trajectories:\nplot(W, X)\ngrid('on')\nxlabel('Catalyst weight, W [kg]')\nplot(W, y)\nylabel('X and y [both dimensionless]')\nlegend(['X', 'y'])\nshow()\n</syntaxhighlight>\n\nEnsure that you obtain a graphical output as shown here. This represents the profile of conversion, \\(X\\) and the pressure drop ratio \\(y\\) throughout the reactor.\n[[Image:Python_output.jpg]]\n\n<!-- === Mac computers ===\n* Start the built-in Mac program called <tt>Terminal</tt>\n* Type the following commands:\n<syntaxhighlight lang=\"bash\">\n# first install \"easy install\" \nsudo easy_install ipython\nsudo easy_install readline\nsudo easy_install numpy\n# Download and install \"gfortran-4.2.3.dmg\"\nsudo easy_install scipy\nipython\n</syntaxhighlight> -->\n\n== MATLAB ==\n\nWith MATLAB we create two files: a file containing the model equation, and another that drives the ODE integrator. Once you have both files on your computer, call <tt>ODEdriver</tt> from the MATLAB command line.\n\n=== Model file: <tt>pbr.m</tt>===\n\n<syntaxhighlight lang=\"MATLAB\">\nfunction d_depnt__d_indep = pbr(indep, depnt)\n \n% Dynamic balance for the packed bed reactor (PBR); demo problem class 05C\n% \n%    indep: the independent ODE variable, such as time or length\n%    depnt: a vector of dependent variables\n% \n%    X = depnt(1) = the conversion\n%    y = depnt(2) = the pressure ratio = P/P_0 = y\n% \n%    Returns d(depnt)/d(indep) = a vector of ODEs\n \n% Assign some variables for convenience of notation\nX = depnt(1);\ny = depnt(2);\n\n% Constant(s)\nFA0   = 0.1362;  % [mol/s]\nkdash = 0.0074;  % [mol/(kg catalyst . s)]\nalpha = 0.0367;  % [1/kg]\neps   = -0.15;   % [-]\n\n% Algebraic equation(s)\nrAdash = -kdash * (1-X)/(1+eps*X) * y;\n\n% Output from this ODE function must be a COLUMN vector, with n rows\nn = numel(depnt);\nd_depnt__d_indep = zeros(n,1);\nd_depnt__d_indep(1) = -rAdash / FA0;\nd_depnt__d_indep(2) = -alpha/(2*y) * (1+eps*X);\n</syntaxhighlight>\n\n=== ODE integrator: <tt>ODEdriver.m</tt> ===\n<syntaxhighlight lang=\"MATLAB\">\n% Integrate the ODE\n% -----------------\n \n% The independent variable always requires an initial and final value:\nindep_start = 0.0;   % kg\nindep_final = 28.0;  % kg\n \n% Set initial condition(s): for integrating variables (dependent variables)\nX_depnt_zero = 0.0;   % i.e. X(W=0) = 0.0\ny_depnt_zero = 1.0;   % i.e. y(W=0) = 1.0\n \n% Integrate the ODE(s):\n[indep, depnt] = ode45(@pbr, [indep_start, indep_final], [X_depnt_zero, y_depnt_zero]);\n \n% Plot the results:\nclf;\nplot(indep, depnt(:,1), 'b')\ngrid('on')\nhold('on')\nplot(indep, depnt(:,2), 'g')\nxlabel('Catalyst weight, W [kg]')\nylabel('X and y')\nlegend('X', 'y')\n</syntaxhighlight>"
                    }
                ]
            },
            "86": {
                "pageid": 86,
                "ns": 0,
                "title": "Steady-state nonisothermal reactor design - 2013",
                "revisions": [
                    {
                        "contentformat": "text/x-wiki",
                        "contentmodel": "wikitext",
                        "*": "{{ClassSidebarYouTube\n| date = 18 March\n| dates_alt_text = \n| vimeoID1 = FnsyFP1QoP0\n| vimeoID2 = 8YVhVL5dVIU\n| vimeoID3 = GJxYdIZ-nCI\n| vimeoID4 = -kMr6kk54Pk\n| vimeoID5 = yLXZSyJK05E\n| vimeoID6 = bAPUtqIWAsI\n| vimeoID7 = 6JIZAhz7Y6o\n| vimeoID8 = xy6N_-Ob4GA\n| vimeoID9 = gaEvAyo29Yw\n| course_notes_PDF = \n| course_notes_alt = Course notes\n| overheads_PDF = \n| assignment_instructions = \n| assignment_solutions = \n| video_download_link_MP4 = http://learnche.mcmaster.ca/media/3K4-2013-Class-10A.mp4\n| video_download_link_MP4_size = 399 M\n| video_notes1 =\n| video_download_link2_MP4 = http://learnche.mcmaster.ca/media/3K4-2013-Class-10B.mp4\n| video_download_link2_MP4_size = 266 M\n| video_notes2 =\n| video_download_link3_MP4 = http://learnche.mcmaster.ca/media/3K4-2013-Class-10C.mp4\n| video_download_link3_MP4_size = 392 M\n| video_notes3 =\n| video_download_link4_MP4 = http://learnche.mcmaster.ca/media/3K4-2013-Class-11A.mp4\n| video_download_link4_MP4_size = 342 M\n| video_notes4 =\n| video_download_link5_MP4 = http://learnche.mcmaster.ca/media/3K4-2013-Class-11B.mp4\n| video_download_link5_MP4_size = 355 M\n| video_notes5 =\n| video_download_link6_MP4 = http://learnche.mcmaster.ca/media/3K4-2013-Class-11C.mp4\n| video_download_link6_MP4_size = 392 M\n| video_notes6 =\n| video_download_link7_MP4 = http://learnche.mcmaster.ca/media/3K4-2013-Class-12A.mp4\n| video_download_link7_MP4_size = 331 M\n| video_notes7 =\n| video_download_link8_MP4 = http://learnche.mcmaster.ca/media/3K4-2013-Class-12B.mp4\n| video_download_link8_MP4_size = 556 M\n| video_notes8 =\n| video_download_link9_MP4 = http://learnche.mcmaster.ca/media/3K4-2013-Class-12C-part1.mp4\n| video_download_link9_MP4_size = 180 M\n| video_notes9 =}}\n__NOTOC__\n\n== Textbook references ==\n\n* F2011: Chapter 11 and 12\n* F2006: Chapter 8\n\n== Suggested problems ==\n\n{| class=\"wikitable\"\n|-\n! F2011\n! F2006\n|-\n| Problem 12-6\n| Problem 8-5\n|-\n| Problem 12-15 (a)\n| Problem 8-16 (a)\n|-\n| Problem 12-16 (a)\n| Problem 8-18 (a)\n|-\n| Problem 12-24 (set up equations) \n| Problem 8-26 (set up equations)\n|}\n\n==Class materials ==\n\n* '''18 March 2013 (10A)''': [http://learnche.mcmaster.ca/media/3K4-2013-Class-10A.mp3 Audio] and [http://learnche.mcmaster.ca/media/3K4-2013-Class-10A.mp4 video] recording of the class.\n* '''20 March 2013 (10B)''': [http://learnche.mcmaster.ca/media/3K4-2013-Class-10B.mp3 Audio] and [http://learnche.mcmaster.ca/media/3K4-2013-Class-10B.mp4 video] recording of the class; and the [[Media:3K4-2013-Class-10B.pdf|handout used in class]].\n* '''21 March 2013 (10C)''': [http://learnche.mcmaster.ca/media/3K4-2013-Class-10C.mp3 Audio] and [http://learnche.mcmaster.ca/media/3K4-2013-Class-10C.mp4 video] recording of the class; and the [[Media:3K4-2013-Class-10C.pdf|revised handout]]\n* '''25 March 2013 (11A)''': [http://learnche.mcmaster.ca/media/3K4-2013-Class-11A.mp3 Audio] and [http://learnche.mcmaster.ca/media/3K4-2013-Class-11A.mp4 video] recording of the class.\n* '''27 March 2013 (11B)''': [http://learnche.mcmaster.ca/media/3K4-2013-Class-11B.mp3 Audio] and [http://learnche.mcmaster.ca/media/3K4-2013-Class-11B.mp4 video] recording of the class.\n* '''28 March 2013 (11C)''': [http://learnche.mcmaster.ca/media/3K4-2013-Class-11C.mp3 Audio] and [http://learnche.mcmaster.ca/media/3K4-2013-Class-11C.mp4 video] recording of the class.\n* '''01 April 2013 (12A)''': [http://learnche.mcmaster.ca/media/3K4-2013-Class-12A.mp3 Audio] and [http://learnche.mcmaster.ca/media/3K4-2013-Class-12A.mp4 video] recording of the class.\n* '''03 April 2013 (12B)''': [http://learnche.mcmaster.ca/media/3K4-2013-Class-12B.mp3 Audio] and [http://learnche.mcmaster.ca/media/3K4-2013-Class-12B.mp4 video] recording of the class.\n* '''04 April 2013 (12C-part1)''': [http://learnche.mcmaster.ca/media/3K4-2013-Class-12C-part1.mp3 Audio] and [http://learnche.mcmaster.ca/media/3K4-2013-Class-12C-part1.mp4 video] recording of the class.\n\n==Source code ==\n\n=== Example on 21 March (class 10C) ===\n\n<tt>'''pfr_example.m'''</tt>\n<syntaxhighlight lang=\"matlab\">\nfunction [d_depnt__d_indep, CA] = pfr_example(indep, depnt, param)\n \n% Dynamic balance for the reactor\n% \n%    indep: the independent ODE variable, such as time or length\n%    depnt: a vector of dependent variables\n%    Returns d(depnt)/d(indep) = a vector of ODEs\n \n% Assign some variables for convenience of notation\nX = depnt(1);\n\n% Constants. Make sure to use SI for consistency\nFT0 = param.FT0;  % mol/s.  Note how we use the \"struct\" variable to access the total flow\nFA0 = 0.9 * FT0;  % mol/s\nT1 = 360;         % K\nT2 = 333;         % K\nE = 65700;        % J/mol\nR = 8.314;        % J/(mol.K)\nHR = -6900;       % J/(mol of n-butane)\nCA0 = 9300;       % mol/m^3\nk_1 = 31.1/3600;  % 1/s (was 1/hr originally)\nK_Cbase = 3.03;   % [-]\n\n% Equations\nT = 43.3*X + param.T_0;  % derived in class, from the heat balance\nk1 = k_1 * exp(E/R*(1/T1 - 1/T)); % temperature dependent rate constant\nKC = K_Cbase * exp(HR/R*(1/T2 - 1/T)); % temperature dependent equilibrium constant\nk1R = k1 / KC; % reverse reaction rate constant\nCA = CA0 * (1 - X); % from the stoichiometry (differs from Fogler, but we get same result)\nCB = CA0 * (0 + X);\nr1A = -k1 * CA; % rate expressions derived in class \nr1B = -r1A;\nr2B = -k1R * CB;\nr2A = -r2B;\nrA = r1A + r2A; % total reaction rate for species A\n\nn = numel(depnt);\nd_depnt__d_indep = zeros(n,1);\nd_depnt__d_indep(1) = -rA / FA0;\n</syntaxhighlight>\n\n'''<tt>driver_ode.m</tt>'''\n<syntaxhighlight lang=\"matlab\">\n% Integrate the ODE\n% -----------------\n \n% The independent variable always requires an initial and final value:\nindep_start = 0.0;  % m^3\nindep_final = 5.0; % m^3\n \n% Set initial condition(s): for integrating variables (dependent variables)\nX_depnt_zero = 0.0;   % i.e. X(V=0) = 0.0\n\n% Other parameters. The \"param\" is just a variable in MATLAB.\n% It is called a structured variable, or just a \"struct\"\n% It can have an arbitrary number of sub-variables attached to it.\n% In this example we have two of them.\nparam.T_0 = 330;         % feed temperature [K]\nparam.FT0 = 163000/3600; % mol/s (was kmol/hour originally)\n\n% Integrate the ODE(s):\n[indep, depnt] = ode45(@pfr_example, [indep_start, indep_final], [X_depnt_zero], odeset(), param);\n\n% Deal with the integrated output to show interesting plots\nX = depnt(:,1);\nT = 43.3.*X + param.T_0;          % what was the temperature profile?\nrA_over_FA0 = zeros(numel(X), 1); % what was the rate profile?\nC_A = zeros(numel(X), 1);         % what was the concentration profile?\nfor i = 1:numel(X)\n    [rA_over_FA0(i), C_A(i)] = pfr_example([], X(i), param);\nend % the above can be done more efficiently in a single line, but the code would be too confusing\n\n% Plot the results\nf=figure;\nset(f, 'Color', [1,1,1])\nsubplot(2, 2, 1)\nplot(indep, X); grid\nxlabel('Volume, V [kg]', 'FontWeight', 'bold')\nylabel('Conversion, X [-]', 'FontWeight', 'bold')\n\nsubplot(2, 2, 2)\nplot(indep, T); grid\nxlabel('Volume, V [kg]', 'FontWeight', 'bold')\nylabel('Temperature profile [K]', 'FontWeight', 'bold')\n\nsubplot(2, 2, 3)\nplot(indep, C_A); grid\nxlabel('Volume, V [kg]', 'FontWeight', 'bold')\nylabel('Concentration C_A profile [K]', 'FontWeight', 'bold')\n\nsubplot(2, 2, 4)\nplot(indep, rA_over_FA0); grid\nxlabel('Volume, V [kg]', 'FontWeight', 'bold')\nylabel('(Reaction rate/FA0) profile [1/m^3]', 'FontWeight', 'bold')\n\n% Now plot one of the most important figures we saw earlier in the course: \n% F_A0 / (-rA) on the y-axis, against conversion X on the x-axis. This plot\n% is used to size various reactors.\n\n% The material leaves the reactor at equilibrium; let's not plot\n% that far out, because it distorts the scale. So plot to 95% of \n% equilibrium\nf = figure; set(f, 'Color', [1,1,1])\nindex = find(X>0.95 * max(X), 1);\nplot(X(1:index), 1./rA_over_FA0(1:index)); grid\nxlabel('Conversion, X [-]', 'FontWeight', 'bold')\nylabel('FA0/(-r_A) profile [m^3]', 'FontWeight', 'bold')\n\n% Updated code to find the optimum inlet temperature:\ntemperatures = 300:3:400;\nconv_at_exit = zeros(size(temperatures));\ni = 1;\nfor T = \n   param.T_0 = T;         % feed temperature [K]\n   [indep, depnt] = ode45(@pfr_example, [indep_start, indep_final], ...\n                       [X_depnt_zero], odeset(), param);\n   conv_at_exit(i) = depnt(end, 1);\n   i = i + 1;\nend\nplot(temperatures, conv_at_exit); grid;\nxlabel('Inlet temperature')\nylabel('Conversion at reactor exit')\n</syntaxhighlight>"
                    }
                ]
            }
        }
    }
}