Chemical Reactor Simulation in C++ – Simulation Example

Chemical Reactor basically simulate the data flow of chemical reagent in a chemical reaction. The process of forming product from some reactant is known as chemical reaction. In a chemical reaction the amount of reactant is decreased with time. On the other hand the amount of product increases with time.

Chemical Reactor Simulation Problem.

There are two reactants A and B. When they form a chemical reaction the amount of A and B decreases. After having a reaction between A and B a product C is formed. So, the amount of C is increased by time.

If the reaction would flow in one direction then it was easy to complete the reaction. But the reaction flows in bi-direction. That means, C is decomposed to A and B also.

The composition rate of C from A and B is = k₁

The decomposition rate of C to A and B is = k₂

The rate of the reaction is dependent to the amount of present substance.

It is given that,

1 mole A + 1 mole B = 2 mole C

If the number of mole of A, B and C is respectively a, b and c at any time dt then,

The change of A and B,

da/dt = db/dt = k₂c – k₁ab

The change of C,

dc/dt = 2k₁ab – 2k₂c

So we get,

da = db = (k₂c – k₁ab) x dt

and,

dc = (2k₁ab – 2k₂c) x dt

For a given value for dt, the change of A, B and C will be as above. That means we have to sum those values to the previous amount to get the new amount.

Note that, da, db will be always negative and dc will be always positive. Because, A and B are decreasing and C is increasing.

Along with the code we have to represent the data flow in graphics. So, we need to use C++ graphics as well.

Chemical Reactor Simulation Code in C++

Here is the simplest code for chemical reactor simulation in C++. Please read the full article to understand everything.

#include<bits/stdc++.h>
#include<graphics.h>
using namespace std;

int main(){
	int gd = DETECT, gm;
    initgraph(&gd, &gm, "");

    double XZero, YZero, XMin, YMin, XMax, YMax;
    XZero=100, YZero=400;
    XMin=0, YMin=0;
    XMax=800, YMax=600;

    line(XMin,YZero,XMax,YZero);        //X Axis
    line(XZero,YMin,XZero,YMax);        //Y Axis

    double AQuantity, BQuantity, CQuantity;
    cout<<"Quantity of A: ";
    cin>>AQuantity;
    cout<<"Quantity of B: ";
    cin>>BQuantity;
    CQuantity=0;

    double K1, K2, dt, T, dAB, dC;
    cout<<"Reaction Constants (K1, K2): ";
    cin>>K1>>K2;
    dt = 0.01;
    T = 5;
    delay(15000);

    for(double t=0; t<T; t=t+dt){

        putpixel(t/dt+XZero , YZero-AQuantity , BLUE );    //A Quantity
        putpixel(t/dt+XZero , YZero-BQuantity , BLUE );    //B Quantity
        putpixel(t/dt+XZero , YZero-CQuantity , BLUE );    //A & B Initial Quantity
        delay(20);

		dAB = K2*CQuantity - K1*AQuantity*BQuantity;
		dC = 2*K1*AQuantity*BQuantity - 2*K2*CQuantity;
        AQuantity = AQuantity + dAB*dt;
        BQuantity = BQuantity + dAB*dt;
        CQuantity = CQuantity + dC*dt;
    }
    getch();
	closegraph();
    return 0;
}

Sample Input:

Quantity of A: 150
Quantity of B: 100
Reaction Constants (K1, K2): 0.01 0.01

Sample Output:

Explanation of Chemical Reactor Code

Here, I have described some simple steps to discuss details explanation of the above code.

Step-1: Basic of C++ and Graphics

As the basic of C++ you need to declare header functions, namespace, main function etc.

As we need graphics so we also need to write the basic code for the graphics features.

Here are those codes.

#include<bits/stdc++.h>
#include<graphics.h>
using namespace std;

int main(){
    int gd = DETECT, gm;
    initgraph(&gd, &gm, "");


    ...............
    ...............
    ...............
    getch();
    closegraph();
    return 0;
}

Step-2: Draw the Coordinate System

In our previous tutorial we have explained the details about drawing the coordinate system. Please view the previous tutorial here to understand drawing coordinate system with basic of simulation in step-1.

The code for coordinate system is the below fragment:

    double XZero, YZero, XMin, YMin, XMax, YMax;
    XZero=100, YZero=400;
    XMin=0, YMin=0;
    XMax=800, YMax=600;

    line(XMin,YZero,XMax,YZero);        //X Axis
    line(XZero,YMin,XZero,YMax);        //Y Axis

Note that, here we must declare the coordinate values as double because we need to perform mathematical operation with double values. If we perform any operation between integer value and double value then some data loss may occur.

Step-3: Take Input A, B and Initialize C

We have to take input of A and B from console as these two are unknown to us. And the amount of C is 0 at starting point. So, we have to write the code as belo:

    double AQuantity, BQuantity, CQuantity;
    cout<<"Quantity of A: ";
    cin>>AQuantity;
    cout<<"Quantity of B: ";
    cin>>BQuantity;
    CQuantity=0;

Step-4: Declare and Assign Constants

In the reaction we have two constants k₁ and k₂. We need to take these two values from console. The higher the values of these constants the faster the reaction will flow.

We have another constant dt which is the unit on the time domain (X axis). And we need a value T which will represent the highest time in time domain.

Below fragment describe it well.

    double K1, K2, dt, T, dAB, dC;
    cout<<"Reaction Constants (K1, K2): ";
    cin>>K1>>K2;
    dt = 0.01;
    T = 5;
    delay(3000);

Here, the delay gives you some time to switch from console to graphics window.

Step-5: Derive A, B and C for 0 to T and Draw Graph

The below part is the most important segment in the code. Here we have continuously derived the current amount of A, B and C and put them as a pixel in graphics window.

    for(double t=0; t<T; t=t+dt){
        putpixel(t/dt+XZero , YZero-AQuantity , BLUE );    //A Quantity
        putpixel(t/dt+XZero , YZero-BQuantity , BLUE );    //B Quantity
        putpixel(t/dt+XZero , YZero-CQuantity , BLUE );    //A & B Initial Quantity
        delay(20);

        dAB = K2*CQuantity - K1*AQuantity*BQuantity;
		dC = 2*K1*AQuantity*BQuantity - 2*K2*CQuantity;
        AQuantity = AQuantity + dAB*dt;
        BQuantity = BQuantity + dAB*dt;
        CQuantity = CQuantity + dC*dt;
    }

Here, the loop runs from 0 to T by incrementing 0.01 (dt) at each step.

Hope you have understood the problem well. If you have any confusion feel free to comment here.