### PGEE11169 Numerical Methods for Chemical Engineers Coursework 1

**Hello, dear friend, you can consult us at any time if you have any questions, add WeChat: daixieit**

**PGEE11169**

**Numerical Methods for Chemical Engineers**

**Coursework 1**

The following coursework is to be submitted to the online dropbox on Learn by **15:00 on Thursday the 13**^{th} **of October 2022**. This is **coursework number 1 of 3**. Each coursework is worth 30% of the final mark.

You should submit the solution to this coursework in the form of a short report on Learn. This report must be typed, must not be longer than 7 pages with at least 2cm margins and a font size of 11pt.

You need to submit also the Python notebook through Noteable. While the report needs to be self-contained, some elements of the implementation are graded in the Python notebook. These elements are clearly stated in the Python notebook. In addition, the Python notbook will be run with the Restart & Run All command and should contain sufficient print and plot statements to check your results.

The marks for this coursework will be awarded for the correctness of the solution as well as for the clarity and logical progression of the solution. Details on laying out mathematical calculations are given on Learn.

1. **Numerical complexity**

Consider the breakthrough simulation from the first lecture. Here, we discretise the column along the z axis with grid spacing of Az. We use three point formulas for the convection and diffusion terms which gives us the following discretisation

dci(t) = d Ci+i — 2ci + Ci_i *3 **c **—* 4ci_i + c《_2 1 — e

*dt* (Az )^{2}^{ }^{V}* 2A**z** **e dt* ，

~~D"~~ = *k* (q*(%) — *q**i* (ci)) , *i **= 1**,...,n*

where *c* and q《are the fluid and adsorbed concentrations, respectively, and e, v, *D *and k are the bed void fraction, flow velocity, diffusion coefficient and mass transfer constant, respectively. The equilibrium adsorbed concentration q* is described by the Langmuir adsorption isotherm.

We arrange the fluid and adsorbed concentrations to the vector x in the following way

^{x} = [^{C}1*, *^{q}*1**, *^{c}*2**, *^{q}*2**, . . . , *^{c}*n**, *^{q}*n* ]

and get the following banded pattern for the right hand side.

(a) In the Python notebook generate the matrix for the right hand side with the [15] create_matrix function for n=4 and with your student ID as seed.

• Print the matrix in the report and describe the structure of the matrix in terms of the non-zero elements.

• Look at the fifth and sixth row of the matrix and describe the non-zero elements in terms of the fluid and adsorbed concentration. Where do they come from?

(b) Use the formula from the lecture slides to calculate the number of elementary [10] operations, i.e. additions/subtractions and multiplications/divisions, required if Gaussian elimination is used to calculate the upper triangular matrix for a range

of sizes up to *n* = 1000. Explain the trend in the required number of elementary operations for the different values of n by referring to the algorithm for naive Gaussian elimination.

(c) Explain how many elementary operations are required if the banded nature is [15] known and used in the calculation of the upper triangular matrix, i.e. only perform operations on non-zero elements. What is the order of the algorithm?

**Remarks:**

• Step through the algorithm column by column and show how many operations are required to eliminate all values below the diagonal in this column.

• Perform the steps for n = 3 and infer the complexity for arbitrary n.

(d) In the Python notebook, write a function that performs Gaussian elimination [10] but takes the banded nature of the problem into account, i.e. you should only operate on non-zero elements. You can use the provided Gaussian elimination function as a template but you need to modify it so that it only operates on non-zero elements.

• Discuss how you are modifying the Gaussian elimination so that you only operate on non-zero elements

• Generate the banded Gaussian elimination function in the Python notebook

(e) In the Python notebook, write a function that generates the matrix and then times how long it takes to calculate the upper triangular matrix with your banded Gauissian implementation as well as with the naive Gaussian implementation for a range of matrix sizes, i.e. matrix size starting at 5 and doubling until you reach 160. The function should return a list of lists with the results.

Print and discuss representative results in the report.

(f) Instead of the single component, we consider a case with two gas phase (c^{1} and c^{2}) and two adsorbed phase (q^{1} and *q** ^{2}*) concentrations. The two species are linked through the multicomponent Langmuir adsorption isotherm

*q**s**c*^{1}

*1* + *b**'**C** ^{1}* + &2C

^{2}

^{ }qsc

^{2}

1 + Sc^{1} + &2C^{2}

How would you arrange the concentrations in the vector x? Explain and discuss your reasoning.

Remark: You don't have to calculate the number of operations but you might want to calculate the bandwidth of the resulting matrix.

2. **Cubic equation of state**

Consider the cubic Redlich-Kwong equation of state

*p RT** **a*

=v - *b* - vTV(vTb)

which relates temperature, pressure and volume of gas. Here the parameters and variables are

• P - gas pressure in [Pa]

• R - gas constant, 8.314 J K^{-1} mol^{-1}

• T - temperature in [K]

• *V* - molar volume [m^{3} mol^{-1}]

• a - constant that corrects for attractive potential of molecules

b - constant that corrects for volume

The constants depend on the critical point of the gas and are calculated by

*b **=* 0.08664 件

P_{c}

Consider the problem of calculating the vapour volume root for water at P = 0.9 bar and *T **= 296* K. The critical temperature and pressure for water are:

*T** _{c}* = 647.1 K

PC = 22.064 MPa

(a) Check that the units are correct and transform the problem into a root finding [15] problem by subtracting the pressure. Write a function in the Python notebook

that defines the root finding problem. Plot the function value for *V* between 0.01 and 0.7 m^{3} mol^{-1}.

(b) In the Python notebook, write a function that performs Newton's method with [15] a finite difference approximation of the derivative. Apply it to the root finding problem starting from two different initial values *V**i* = 0.01 m^{3} mol^{-1} and V2 =

0.7 m^{3} mol^{-1}. Comment and explain the behaviour of Newton's method.

2022-10-10