thermkin.xmds

XMDS is a code generator that integrates equations. You write them down in human readable form in a

<?xml version="1.0"?>
<simulation>
  
<!-- $Id: thermkin.xmds,v 1.2 2004/08/03 05:29:49 paultcochrane Exp $ -->

<!--  Copyright (C) 2000-2004                                           -->
<!--                                                                    -->
<!--  Code contributed by Greg Collecutt, Joseph Hope and Paul Cochrane -->
<!--                                                                    -->
<!--  This file is part of xmds.                                        -->
<!--                                                                    -->
<!--  This program is free software; you can redistribute it and/or     -->
<!--  modify it under the terms of the GNU General Public License       -->
<!--  as published by the Free Software Foundation; either version 2    -->
<!--  of the License, or (at your option) any later version.            -->
<!--                                                                    -->
<!--  This program is distributed in the hope that it will be useful,   -->
<!--  but WITHOUT ANY WARRANTY; without even the implied warranty of    -->
<!--  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the     -->
<!--  GNU General Public License for more details.                      -->
<!--                                                                    -->
<!--  You should have received a copy of the GNU General Public License -->
<!--  along with this program; if not, write to the Free Software       -->
<!--  Foundation, Inc., 59 Temple Place - Suite 330, Boston,            -->
<!--  MA  02111-1307, USA.                                              -->

  <name> thermkin </name>      <!-- the name of the simulation -->
  
  <author> Paul Cochrane </author>  <!-- the author of the simulation -->
  <description>
    <!-- a description of what the simulation is supposed to do -->
    Example simulation of a thermokinetic oscillator
    The reaction scheme is
		A ----> X       (1)
		X ----> R       (2)
    Reaction 1 is neither exothermic nor endothermic and
    its activation energy is equal to zero.  Reaction 2
    is exothermic with a reaction enthalpy deltarH, and 
    it has a substantial activation energy Ea. Assume 
    that the concentration of component A is constant. 
    The rate equation for X is
	d[X]/dt = k1 [A] - k2(T) [X].
    Express the Arrhenius equation for Reaction 2 in terms
    of the rate constant at T0, the temperature of the
    surroundings and deltaT = T - T0,
	k2(T) = k2(T0) Exp[alpha deltaT/(1 + deltaT/T0)],
    where
	alpha = Ea/(R T0^2).
    Then the rate equation for X is
	d[X]/dt = k1 [A] - k2(T0)[X] Exp[alpha deltaT/(1 + deltaT/T0)].
    The rate equation for deltaT at time t is
	d[deltaT]/dt = gamma [X] Exp[alpha deltaT/(1 + deltaT/T0)] - beta deltaT,
    in which
	beta = A h/(V rho cp)
	gamma = |deltarH| k2(T0)/(rho cp),
    with rho and cp the density and specific heat of the 
    reaction mixture,  A the area for heat transfer from 
    the system to the surroundings, h the heat transfer 
    coefficient, and V the volume of the system.  For more 
    on this system see P. Gray and S.K. Scott, "Chemical 
    Oscillations and Instabilities," Oxford Univ. Press, 
    1990, Chap. 4.

    Adapted for xmds from "Mathematica computer programs for physical
    chemistry", William H. Cropper, Springer Verlag (1998)
  </description>
  
  <!-- Global system parameters and functionality -->
  <prop_dim> t </prop_dim>    <!-- name of main propagation dim -->
  
  <error_check> yes </error_check>   <!-- defaults to yes -->
  <use_wisdom> yes </use_wisdom>     <!-- defaults to no -->
  <benchmark> yes </benchmark>       <!-- defaults to no -->
  <use_prefs> yes </use_prefs>       <!-- defaults to yes -->
  
  <!-- Global variables for the simulation -->
  <globals>
  <![CDATA[
    // rate constants for the two
    // reactions in s^-1, with k2 given a value k20 for
    // T = T0
    const double k1 = 0.1;
    const double k20 = 0.5;

    // a, the constant concentration
    // of component A in mol L^-1.  Try values of a in
    // the range .3 to 1.5 mol L^-1.e
    const double a = 0.6;

    // initial values for [X] in mol L^-1 and deltaT in K
    const double CX0 = 0.0;
    const double deltaT0 = 100.0;

    // temperature of the surroundings T0
    const double T0 = 300.0;

    // the parameter alpha = Ea/(R T0^2)
    // in K^-1, assuming that Ea = 166 kJ mol^-1 and T0 = 
    // 400 K
    const double Alpha = 0.1248;

    // the parameter beta = A h/(V rho cp)
    // in s^-1, assuming that A = .05 m^2, h = 30 J m^-2 K^-1
    // s^-1, V = .001 m^-3 and (rho cp) = 150 J m^-3 K^-1
    const double Beta = 10.0;

    // the parameter gamma = |deltarH|*
    // k20/(rho cp) in m^3 mol^-1 K s^-1, assuming that 
    // |deltarH| = 400 kJ mol^-1, k20 = .5 s^-1 and (rho cp) = 
    // 150 J m^-3 K^-1
    const double Gamma = 1333.0;
  ]]>
  </globals>
  
  <!-- Field to be integrated over -->
  <field>
    <name> main </name>
    <samples> 1 </samples>       <!-- sample 1st point of dim? -->
    
    <vector>
      <name> main </name>
      <type> double </type>           <!-- data type of vector -->
      <components> CX deltaT </components>       <!-- names of components -->
      <![CDATA[
        CX = CX0;
	deltaT = deltaT0;
      ]]>
    </vector>
  </field>
  
  <!-- The sequence of integrations to perform -->
  <sequence>
    <integrate>
      <algorithm> RK4IP </algorithm> <!-- RK4EX, RK4IP, SIEX, SIIP -->
      <interval> 20 </interval>   <!-- how far in main dim? -->
      <lattice> 1000000 </lattice>     <!-- no. points in main dim -->
      <samples> 1000 </samples> <!-- no. pts in output moment group -->
      
      <![CDATA[
        dCX_dt = k1*a - k20*CX*exp(Alpha*deltaT/(1.0 + deltaT/T0));
	ddeltaT_dt = Gamma*CX*exp(Alpha*deltaT/(1.0 + deltaT/T0)) - Beta*deltaT;
      ]]>
    </integrate>
  </sequence>
  
  <!-- The output to generate -->
  <output format="ascii">
    <group>
      <sampling>
        <lattice> 1000 </lattice>           <!-- no. points to sample -->
        <moments> X delT </moments>           <!-- names of moments -->
        <![CDATA[
          X = CX;
	  delT = deltaT;
        ]]>
      </sampling>
    </group>
  </output>
  
</simulation>