Modeling & Simulation of Chemical Reaction | Continuous System Simulation | C++ Implementation

Chemical reactions exhibit dynamic equilibrium, which means that a combination reaction is also accomplished by the reverse process of decomposition reaction. At the steady state the rates of the forward and the backward reaction is same.

In addition to the forward reaction where Ch1 and Ch2 react to produce Ch3, there may also be backward reaction, where by, Ch3 decomposes back into Ch1 and Ch2. Let the rate of formation of Ch3 be proportional to the product of the amounts Ch1 and Ch2 present in the mixture and let the rate of decomposition of Ch3 be proportional to its amount in the mixture. 

Let us consider C1, C2  and C3 Amount of Ch1, Ch2 and Ch3 at any instant of time t, the rate of increase of C1, C2& C3 are given by the following differential equations:

dC1/dt = K1C3– K1C1C2

dC2/dt = K2C3– K1C1C2

dC3/dt=2K1C1C2-2K2C3

Where K1and K2 are constants.

As soon as the chemicals Ch1 and Ch2are mixed, the reaction starts and the amount of C1, C2, C3 in the mixture goes on changing as time progresses.The simulation of reaction will determine the state of the system. i.e. value of quantities C1, C2 and C3at different points in time. Starting at zero time, a very small increment of time is taken in each step. It is assumed to be so small, that all changes in the mixture can be taken to occur at theend of each increment. If C1(t), C2(t) and C3(t) are the quantities to there chemicals at time t, then at time t +  Δt, the quantities are ;

C1(t +  Δt) = C1(t) + dC1(t)/dt

C2(t +  Δt) = C2(t) + dC2(t)/dt

C3(t + Δt) = C3(t) + dC3(t)/d

After forming the the mathematical model, we can simulate the system using following program written in C++:

// Simulate Chemical Reaction,  +Jivan Nepali  

 #include<iostream>  
 #include<fstream>  
 using namespace std;  
 float k1 = 0.008,k2 = 0.002,del_t = 0.1;  
 float calc_c1(float c1, float c2, float c3)  
 {  
   return (c1+(k2*c3-k1*c1*c2)*del_t);  
 }  
 float calc_c2(float c1, float c2, float c3)  
 {  
   return (c2+(k2*c3-k1*c1*c2)*del_t);  
 }  
 float calc_c3(float c1, float c2, float c3)  
 {  
   return (c3+(2*k1*c1*c2-2*k2*c3)*del_t);  
 }  
 int main()  
 {  
   float c1 = 25,c2 = 80, c3 =0;  
   int i;  
   float temp_c1,temp_c2,temp_c3;  
   ofstream outfile("simulation_lab.txt");  
   outfile<<c1<<","<<c2<<","<<c3<<"\n";  
   for(i=1;i<400;i++)  
   {  
     temp_c1 = calc_c1(c1,c2,c3);  
     temp_c2 = calc_c2(c1,c2,c3);  
     temp_c3 = calc_c3(c1,c2,c3);  
     c1 = temp_c1;  
     c2 = temp_c2;  
     c3 = temp_c3;  
     outfile<<c1<<","<<c2<<","<<c3<<"\n";  
   }  
   outfile.close();  
   cout<<"\nAll data printed on the file \"Simulation_Lab.txt\" on your project folder.\n";  
   cout<<"\nHave a look!!!\n\n";  
   return 0;  
 }  

To simulate the result obtained through the above model, we can use following MATLAB code:

%Simulate Chemical Reaction, +Jivan Nepali  

data = load('simulation_lab.txt');  
 plot(data(:,1),'b-');  
 hold on;  
 plot(data(:,2),'r-');  
 hold on;  
 plot(data(:,3),'k-');  
 legend('reactant c1','reactant c2','product c3');  

The graph obtained using the MATLAB will be as shown:

Enjoy Simulating....!!!