xmds 1.6.6-4 / usr / share / doc / xmds / examples / chain.xmds

<?xml version="1.0"?>
<simulation>
  
<!-- $Id: chain.xmds 1526 2007-08-21 17:30:14Z paultcochrane $ -->

<!--  Copyright (C) 2000-2007                                           -->
<!--                                                                    -->
<!--  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., 51 Franklin St, Fifth Floor, Boston,            -->
<!--  MA  02110-1301, USA.                                              -->

  <name> chain </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 free-radical chain polymerisation.

    Adapted for xmds from "Mathematica computer programs for physical
    chemistry", William H. Cropper, Springer Verlag (1998)

    Calculates [M], concentration of monomer, [P.], concentration of
    chain radicals, [I], concentration of initiator, [R.], 
    concentration of initiator radicals, and [P], 
    concentration of polymer, during the first 10 s
    of a free-radical chain polymerization process.  
    The various rates involved are

	Rate of initiation: Ri = 2 f kd [I]
	Rate of propagation: Rp = kp [M] [P.] + kp [M] [R.]
	Rate of termination: Rt = 2 kt [P.]^2,

    where kp, kd and kt are rate constants and f is the 
    initiation efficiency factor. Rate equations are:

	d[M]/dt = -kp [M] [P.] - kp [M] [R.]
	d[P.]/dt = 2 f kd [I] - 2 kt [P.]^2		
	d[I]/dt = -2 f kd [I]
	d[R.]/dt = 2 f kd [I] - kp [M] [R.]
	d[P]/dt = kt [P.]^2.

    The program also calculates the average degree of polmerization
    according to Eq. (10.47), as given in 
    "Mathematica computer programs for physical chemistry"

	xave = -(d[M]/dt)/(d[P]/dt).
  </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 (in 1/s for unimolecular reactions, 
    // and L/mol/s for bimolecular reactions)
    // Values of kp and kt are for methyl acrylate as the 
    // monomer at 25oC, and the value of kd is typical.
    const double kp = 1580.0;
    const double kd = 1e-5;
    const double kt = 2.75e7;
    
    // initiation efficiency factor
    const double f = 1.0;

    // initial concentrations (in mol/L)
    const double Monomer0 = 0.8;
    const double Chain0 = 0.0;
    const double Initiator0 = 0.001;
    const double Radical0 = 0.0;
    const double Polymer0 = 0.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> Monomer Chain Initiator Radical Polymer </components>       <!-- names of components -->
      <![CDATA[
        Monomer = Monomer0;
	Chain = Chain0;
	Initiator = Initiator0;
	Radical = Radical0;
	Polymer = Polymer0;
      ]]>
    </vector>
  </field>
  
  <!-- The sequence of integrations to perform -->
  <sequence>
    <integrate>
      <algorithm> RK4IP </algorithm> <!-- RK4EX, RK4IP, SIEX, SIIP -->
      <interval> 10 </interval>   <!-- how far in main dim? -->
      <lattice> 100000 </lattice>     <!-- no. points in main dim -->
      <samples> 1000 </samples> <!-- no. pts in output moment group -->
      <vectors>main</vectors>
      <![CDATA[
        dMonomer_dt = -kp*Monomer*Chain - kp*Monomer*Radical;
	dChain_dt = 2.0*kd*f*Initiator - 2*kt*Chain*Chain;
	dInitiator_dt = -2.0*kd*f*Initiator;
	dRadical_dt = 2*f*kd*Initiator - kp*Monomer*Radical;
	dPolymer_dt = kt*Chain*Chain;
      ]]>
    </integrate>
  </sequence>
  
  <!-- The output to generate -->
  <output format="ascii">
    <group>
      <sampling>
        <moments> 
	  MonomerOut ChainOut InitiatorOut RadicalOut PolymerOut xAveOut
	</moments>           <!-- names of moments -->
        <![CDATA[
	  // output stuff, scaled nicely for plotting
          MonomerOut = Monomer;
	  ChainOut = 1e7*Chain;
	  InitiatorOut = 1e3*Initiator;
	  RadicalOut = 1e10*Radical;
	  PolymerOut = 1e7*Polymer;
	  xAveOut = (kp*Monomer*Chain + kp*Monomer*Radical)/(kt*Chain*Chain)/5000.0;
        ]]>
      </sampling>
    </group>
  </output>
  
</simulation>