A Newtonian model for a simple ionic reaction

METRIC Project

This problem forms part of your coursework, and should be handed in to your tutor by Friday 27th March. Hand in your answers, fully annotated to your tutor. Discuss with him/her in what format it is required. From past experience you would be well advised to save multiple copies of your work, frequently.

1. Objectives

In this exercise, you are asked to set up and explore a simulation of the reaction between a free \(TraditionalForm\`\(Cl\^-\)\) ion and a bonded NaCl ion pair. You're going to be using a Newtonian model: that is to say, the reaction is modelled as an interaction between charges moving under forces, and their motion is predicted using Newton's second law, F = ma. For simplicity, the charges will move and lie along a straight line, which allows you to treat the problem as a one-dimensional one.

This reaction is a very simple one in chemical terms, but it is still too complex to be predicted using traditional pencil-and-paper mathematics, or indeed any of Mathematica's exact, symbolic functions. Instead, you'll be using a method of numerical approximation, which will of necessity be an imperfect one. One of your objectives is to find out how well this numerical approximation seems to model physical reality, by finding out whether energy appears to be conserved within the model.

Depending on the positions and velocities of the three ions when they begin interacting, sometimes a reaction will occur and sometimes it will not. A third outcome is also possible: under certain conditions, all three ions continue to interact with each other, with no one ion becoming dissociated from the other two. Your other main objective is to set the initial conditions which will produce these three outcomes and devise a way of deciding which has occurred.

2. Preliminaries

In this section you'll find the Mathematica code you'll need before you can set up your simulation. However, you should note that we have set up the code to make a movie of the reaction, but for your report we want you to alter the code to produce plots of position as a function of elapsed time. To save you a lot of typing we have placed this  Mathematica file, MathLabQue2.nb, on the hard disk of the Maths Lab computers in the Mathematica Folder.

Physical constants

The first priority is to inform Mathematica of the values of the various physical constants that the simulation requires. We used the Basic Input Palette to produce the suffices, powers and Greek letters you see here.

The permittivity of the vacuum:

(Subscript[ϵ, 0]) =  8.85419*^-12

Avogadro's number:

(Subscript[N, A]) =  6.02205*^23

The relative atomic masses of sodium and chlorine:

(Subscript[RAM, Na]) =  22.9898(Subscript[RAM, Cl]) = 34.9688

The charge, in Coulombs, on an electron:

e  =  1.60219*^-19

The equilibrium bond length of the NaCl dipole, in metres:

(Subscript[r, 0]) =  2.5*^-10

You can now use the relative atomic masses, together with Avogadro's number, to calculate the masses of the sodium and chloride ions in kg.

(Subscript[m, Na]) =  ((Subscript[RAM, Na])/((10^3) (Subscript[N, A])))(Subscript[m, Cl]) =  ((Subscript[RAM, Cl])/((10^3) (Subscript[N, A])))

Attraction and repulsion forces

Since, as you'll see, the two chloride ions never approach each other very closely, we can model the repulsive force between them using a simple Coulomb (electrostatic) repulsion. In newtons, for two negative charges separated by a displacement r metres:

(Subscript[F, mm])[r_]  = 
      -((e^2)/(4 π (Subscript[ϵ, 0 ]) (r^2)))

Between the sodium ion and either one of the sodium ions, on the other hand, there is a force that is attractive when the ions are far apart and repulsive when they are close together. This force is extremely well modelled by the following formula:

(Subscript[F, pm])[r_]  = 
      ((e^2)/(4 π (Subscript[ϵ, 0]))) 
        ((1/(r^2)) - (((((Subscript[r, 0])))^8)/(r^10)))

Time

You need to tell Mathematica the duration of the simulated reaction, in seconds. For different reactions under different sets of starting conditions, different durations may well be appropriate, but here is a suggestion.

T = 1.2*^-12

Initial values of the variables

Initial positions (in metres) and velocities (in metres per second) of the three ions. By changing these values you can observe different behaviours.

(Subscript[x1, 0]) = -5.5 (Subscript[r, 0])(Subscript[x2, 0]) =  0.0(Subscript[x3, 0]) =  1.1  (Subscript[r, 0])(Subscript[v1, 0]) =  5000.0(Subscript[v2, 0]) =  0.0(Subscript[v3, 0]) =  0.0

You should save your notebook at this point so that you can retrieve the preliminary code without retyping it.

3. Simulating the reaction

Solving the equations of motion

Each of the three ions moves under the influence of the other two. More precisely, the force on each, and therefore the acceleration of each, depends on its position relative to the other two. This gives three coupled second-order equations of motion, which we solve approximately using Mathematica's NDSolve function, subject to the initial conditions specified above.

\(\(xrules = 
      NDSolve[{
			
          \(\(x1'\)'\)[t] == 
            \(1\/m\_Cl\((F\_pm[x2[t] - x1[t]]\  + \ F\_mm[x3[t] - x1[t]])
              \)\), 
			
          \(\(x2'\)'\)[t] == 
            \(1\/m\_Na\) 
              \((\ \(-F\_pm[x2[t] - x1[t]]\)\  + \ F\_pm[x3[t] - x2[t]])\), 

          			\(\(x3'\)'\)[t] == 
            \(1\/m\_Cl\) 
              \((\(-F\_mm[x3[t] - x1[t]]\)\  - \ F\_pm[x3[t] - x2[t]])\), 

          			x1[0] == x1\_0, x2[0] == x2\_0, x3[0] == x3\_0, 
			
          \(x1'\)[0] == v1\_0, \(x2'\)[0] == v2\_0, \(x3'\)[0] == v3\_0}, 

        		{x1, x2, x3}, {t, 0, T}, MaxSteps -> 2000]; \)\)

The solutions are given in the form of three "hidden" functions of t, defined only for 0 ≤ t ≤ T.

Positions and velocities in physical space

We need a way of getting actual values out of the three "hidden" functions that came out of the call to NDSolve.

pos1[t_]  :=  (x1[t] /. xrules)[[1]]pos2[t_]  :=  (x2[t] /. xrules)[[1]]pos3[t_]  :=  (x3[t] /. xrules)[[1]]vel1[t_]  :=  (x1'[t] /. xrules)[[1]]vel2[t_]  :=  (x2'[t] /. xrules)[[1]]vel3[t_]  :=  (x3'[t] /. xrules)[[1]]

You should satisfy yourself that these definitions have worked, by reading values with commands like

pos1[0.5  T]vel3[0.3  T]

and so on.

Creating an animation

You are now in a position to create a 401-frame animation of the motion of the three ions:

\(threeIonFrame[t_]\  := 
      Show[Graphics[{AbsolutePointSize[15], 

            				{RGBColor[0.8, 0.8, 0], Point[{pos1[t], 0}]}, {
              RGBColor[0.4, 0.4, 0.4], Point[{pos2[t], 0}]}, {
              RGBColor[0.8, 0.8, 0], Point[{pos3[t], 0}]}}], 
		
        PlotRange -> {{\(-15\) r\_0, 15  r\_0}, Automatic}]\)

Table[threeIonFrame[t], {t, 0, T, T/100}];

To view this animation, click on the cell bracket that groups all the frames together and choose Animate Selected Graphics from the Cell menu.

4. The problem

(i) Recoding

Your first task is to re-code the example program so that you can plot the positions of each ion as a function of time. Make sure that your plots are clearly and unambiguously labelled.

(ii) The three possible outcomes

Now find initial conditions for which a reaction appears to occur; for which it does not appear to occur; and a third set for which neither chloride ion appears to dissociate, and the three-ion interaction appears to continue indefinitely.

You should feel free to change the values of (Subscript[x1, 0]), (Subscript[x2, 0]), (Subscript[x3, 0]) and (Subscript[v1, 0]), but you should probably leave (Subscript[v2, 0]) and (Subscript[v3, 0]) set at zero in the first instance.

(iii) Energy conservation

The total mechanical energy of the system is formed of the kinetic energies of each of the three particles (which you can calculate using the velocity functions you have defined) and the potential energy that arises from their relative positions in the force field (which you can calculate using the position functions you have defined, based on what you know about potentials).

The simulation is based on a numerical approximation. In a fundamental sense, then, it's not to be trusted, since no numerical method is exact. If we found, for example, that energy appeared not be be conserved, at least approximately, in our model, then we should reject it.

Investigate how well energy appears to be conserved within this model. Report on the performance of the simulation in this regard.

(iv) A test for dissociation

In part (ii) of this problem, you uncovered three different types of apparent behaviour of the system. However, there's a problem. How can we tell, if an ion appears to be dissociated, whether it really is?

Devise a test for deciding whether, when dissociation appears to have occurred, it really has. Your test need not be perfect but it should take you further than mere observation. Apply your test to your "reaction" and "no reaction" examples from part (ii), and comment on your findings