## Mass Conservation Equations

1. A stirred vessel is used to reduce the concentration of a contaminant in a liquid by adsorption onto particles of activated carbon. Its operation can be described by the following mass conservation equations

dC

dt = –K1(qeq – q)

dq

dt = K2(qeq – q)

Where C (in kg m–3) is the concentration of the contaminant in the liquid and q (in kg kg–1) is the quantity of contaminant adsorbed onto the activated carbon. qeq (in kg kg–1) is the equilibrium concentration adsorbed onto the activated carbon given by a piecewise-defined adsorption isotherm

8 qeq = 0.35C C 5 kg m–3

1 + 0.75C

< qeq = 0.246C0.25 C > 5 kg m–3

:

The model parameters are K1 = 5 10–2 kg m–3 s–1 and K2 = 5 10–4 s–1. At t = 0 s,

C = 50 kg m–3 and q = 0 kg kg–1.

Question 1 Solution

clc, clear all, close all

%% Question 1

K1 = 5e-2;

K2 = 5e-4;

C = 50;

q = 0;

% Define initial condition

C0 = 50;

q0 = 0;

% Define time of somulation and step

tf = 1*3600; % 3600 seconds or 1 hour

N = 10000; % Number of points to use for the time window

tspan = linspace(0, tf, N); % 10000 points between 0 and 3600

[t, x] = ode45(@(t,x)[-K1*(qeq(x(1))-x(2));K2*(qeq(x(1)) - x(2))], tspan, [C0, q0]);

C = x(:,1);

q = x(:,2);

% Plot solution

figure

subplot(1,2,1)

plot(tspan, C)

xlabel('Time (s)')

ylabel('{C(t)}')

grid on

title('{C(t) vs. time}')

subplot(1,2,2)

plot(tspan, q), grid on

xlabel('Time (s)')

ylabel('{q(t)}')

title('{q(t) vs. time}')

% Display the values at t = 1h in COmmand Window

fprintf("The value of C(t) at t = 1 h is: %.4f kg/m^3\n", C(end));

fprintf("The value of q(t) at t = 1 h is: %.4f kg/kg\n", q(end));

a. Develop a function in MATLAB that takes values of C as input argument and returns values of qeq as output according to the adsorption isotherm.

[5 marks]

1(a) Solution

function ret = Question1a_10601875(C)

% This function is for Question 1, Part a)

if C <= 5

ret = 0.35*C/(1+0.75*C);

else % C > 5

ret = 0.246*C^(0.25);

end

end

## MATLAB Graph to Represent Evolution of the Concentration

b. Develop a program in MATLAB to represent graphically the evolution of the concentration of the contaminant in the liquid during the first hour of operation. The plot should have appropriate labels in both axes with SI units. The program should also write in the Command Window the value of the adsorbed concentration at t = 1 h, indicating as well the SI units.

[20 marks]

1(b) Solution

clc, clear all, close all

%% Question 1

K1 = 5e-2;

K2 = 5e-4;

% Define initial condition

C0 = 50;

q0 = 0;

% Define time of somulation and step

tf = 1*3600; % 3600 seconds or 1 hour

N = 10000; % Number of points to use for the time window

tspan = linspace(0, tf, N); % 10000 points between 0 and 3600

[t, x] = ode45(@(t,x)[-K1*(Question1a_10601875(x(1))-x(2));K2*(Question1a_10601875(x(1)) - x(2))], tspan, [C0, q0]);

C = x(:,1);

q = x(:,2);

% Plot solution

figure

subplot(1,2,1)

plot(tspan, C)

xlabel('Time (s)')

ylabel('{C(t)}')

grid on

title('{C(t) vs. time}')

subplot(1,2,2)

plot(tspan, q), grid on

xlabel('Time (s)')

ylabel('{q(t)}')

title('{q(t) vs. time}')

% Display the values at t = 1h in COmmand Window

fprintf("The value of C(t) at t = 1 h is: %.4f kg/m^3\n", C(end));

fprintf("The value of q(t) at t = 1 h is: %.4f kg/kg\n", q(end));

## Van Der Waals Equation of State

2. The van der Waals equation of state is used by chemical engineers to predict properties of fluids. It is given by

p = RT – a

v – b v2

Where p is the pressure (in Pa), T the absolute temperature (in K) and v the molar volume of the fluid (in m3 mol–1). The parameter R = 8.314 J K–1 mol–1 is the ideal gas constant, and a = 3pcvc2 and b = vc/3 are model parameters, where vc and pc are the molar volume and pressure of the fluid at the critical point conditions, respectively.

The p – v curve for methane at T = 128 K is shown in Figure 2.1. For this compound, pc = 4.6 106 Pa and vc = 9.86 10–5 m3 mol–1. The vapour pressure, pvap, is defined as the pressure at which the areas A1 and A2 are equal (see Figure 2.1).

Figure 2.1: van der Waals isotherm for methane at T = 128 K.

a. The van der Waals equation can be expressed in polynomial form as

pv3 – (pb + RT )v2 + av = ab

Develop a MATLAB program to calculate the three possible solutions for the molar volume at a given value of pressure and T = 128 K. The program should return the three values in the Command Window with their SI units.

Hint: The three values can be calculated simultaneously using a particular built-in MATLAB function.

2(a) Solution [5 marks]

clc

clear all

%% Question 2, Part a)

R = 8.314; % ideal gas constant (J/k/mol)

pc = 4.6e6; % Pa

vc = 9.86e-5; % m^3 /mol

T = 128; % in Kelvin

% Calculate a and b

a = 3*pc*vc^2;

b = vc/3;

% Define the value of pressure (pavp)

p = 2.75e6;

% Calculate using roots function

% pv^3 - (pb+RT)v^2 + av - ab

v = roots([p, -(p*b + R*T), a, -a*b]);

fprintf("The three values of volume for T = %.0fK and p = %d Pa are:\n", T, p);

fprintf("v1 = %.6f m^3 /mol\nv2 = %.6f m^3 /mol\nv3 = %.6f m^3 /mol\n", v(1), v(2), v(3));

% Iterate v to get the graph of p vs. v

vspan = 0:1e-6:3e-4;

pspan = R*T./(vspan-b) - a./vspan.^2;

b. Develop a MATLAB function that takes as input argument a guess for the vapour pressure of methane at T = 128 K and returns the difference between the areas A1 and A2.

[10 mark]

Question 2(b) Solution

function ret = Question2b_10601875(p)

% Parameters

R = 8.314; % ideal gas constant (J/k/mol)

pc = 4.6e6; % Pa

vc = 9.86e-5; % m^3 /mol

T = 128; % in Kelvin

% Calculate a and b

a = 3*pc*vc^2;

b = vc/3;

% Calculate roots

v = roots([p, -(p*b + R*T), a, -a*b]);

% sort ascending

v = sort(v);

% Area A1 is between v(1) and v(2), area A2 is between v(2) and v(3)

N = 10000; % N equally-spaced points between intervals

h1 = (v(2)-v(1))/N;

h2 = (v(3)-v(2))/N;

vspan1 = v(1):h1:v(2);

vspan2 = v(2):h2:v(3);

p1 = R*T./(vspan1-b) - a./vspan1.^2;

p2 = R*T./(vspan2-b) - a./vspan2.^2;

% area under curve using trapz

A1 = trapz(vspan1, p1);

A2 = trapz(vspan2, p2);

% Difference in areas

ret = A2-A1;

end

c. Develop a program in MATLAB that uses the function developed in Question 2b to calculate the vapour pressure of methane at T = 128 K. The program should return the value of pvap in the Command Window with its SI units. Limit your search between the values of pressure 1.5 106 and 3 106 Pa.

Note: If you were not able to solve Question 2b, you can still solve Question 2c by assuming that you have a function named diffarea in the Current Folder that takes as input argument a guess for the vapour pressure of methane at T = 128 K and returns the difference between areas A1 and A2.

[10 marks]

2( C) Solution

clc

clear all

% Define limits of pressure

lb = 1.5e6;

ub = 3e6;

% Calculate difference in areas and select points for minimum difference

values = [];

pspan = lb:50:ub;

for p = pspan

Adiff = Question2b_10601875(p); % this is diffarea

values = [values, abs(Adiff)];

end

% Plot difference in areas

figure

plot(pspan, values)

grid on

xlabel('Pressure (Pa)');

ylabel('{|A_2-A_1|}')

title('Difference in areas');

[min_val, loc] = min(values);

fprintf("The pressure where the difference in areas is minimum is: %.0f Pa\n", pspan(loc))