Simulation Testbed for Liquid Chromatography

(1)

Simulation testbed

for liquid chromatography

David Andersson

Master’s Thesis

Master of Science in Engineering Physics

Spring 2021

(2)

Simulation Testbed for Liquid Chromatography David Andersson, daan0094@student.umu.se

Supervisors: Rickard Sjögren Sartorius

Brandon Corbett Sartorius

Daniel Nilsson Dept. of Physics

Examiner: Åke Brännström Dept. of Mathematics and Mathematical Statistics

Master of Science Thesis in Engineering Physics, 30 ECTS Department of Physics

Umeå University

SE-901 87 Umeå, Sweden

(3)

Abstract

The stlc package is proposed as a tool for simulation of liquid chromatography by implementing several lumped kinetic models which combine diffusive mass transport and adsorption isotherm equations. The purpose of the package is to provide computationally eﬀicient approximations to the general rate model of chromatography. Orthogonal collocation is used to discretize the spatial domain and the resulting system of ordinary differential equations is evaluated by one of several solvers made available in the package. Comparisons between numerical and analytical Laplace domain solutions for values of mass transfer coeﬀicient, 𝑘, ranging from 0 to 1000 and lumped dispersion constant values, 𝐷_𝐿, from 10⁻⁵ to 10⁻²are presented. Analytical results were approximated to an 𝐿¹error in the range 10⁻⁵to 10⁻³with a maximum evaluation time of 0.27𝑠 for 100 grid points. The breakthrough curves of the analytical solution are accurately recreated indicating a correct implementation. Variations in accuracy can be partly attributed to oscillations induced by steep gradients in the solution. The oscillations are reduced by the addition further points to the spatial grid. The package is implemented in Python using minimal dependencies and can produce approximations with short evaluation times.

The Python programming language is dynamically typed and uses automatic memory management, properties which can improve productivity and be beneficial to research applications. The addition of this package to the extensive Python ecosystem of libraries can potentially aid future developments in chromatography.

(4)

1 Chromatography

Chromatography is a technique used to separate, purify and identify chemical components of a mixture. The technique can be split into several modes that rely on the same core process. In chromatography one or more substances are dissolved in an medium called the mobile phase. By letting the mobile phase transport over a solid surface, sub- stance components are separated from each other. Mobile phase components are separated by binding to another compound called the solid phase. The solid phase may take the form of a pellet or a spherical particle, objects which can be made porous to increase their surface area.

Commonly in chromatography several of these particles are packed densely in a cylindrical column. The purpose of the stationary phase is to slow down select components in the mobile phase by means of chemical or physical interactions. When the mobile phase transports over the solid phase, mobile phase components with an aﬀinity for the solid phase will be slowed down. Varying aﬀinity is caused by different types of bonds, physical interactions, component solubility and other molecular level phenomena. When the mobile phase leaves the proximity of the stationary phase, outlet concentration of the transported components will vary with time. The goal often being complete separation of one or more components in the initial mixture from the others.

Chromatography sees wide use in many industries as a separation and purification technique. Chemical analysis is a natural application of chromatography because of the capability to separate components in complex mixtures.

As a purification technique, chromatography is coveted due to the ability to separate products from reagents in synthesis. Active components in a drug can be purified to the level required for medicinal treatment. Liquid chromatography, where the mobile phase consists of aqueous or liquid solution, is an important tool in the production of biotherapeutics. In biotherapeutics separation of biologi- cal components such as proteins from the solution in which they grow makes up a large part of the production process.

The means by which proteins are separated define several modes of liquid chromatography. Some of these modes include hydrophobic interaction, hydrophilic interaction, ion-exchange, aﬀinity and size-exclusion chromatography (Zhu and Chen 2017).

In the biopharmaceutical industry, therapeutic proteins are artificially constructed in laboratories for medicinal applications. Advanced treatments like immunotherapy for cancer fall under the use cases for therapeutic proteins. Here the use of liquid chromatography in the development and production of therapeutic proteins such as monoclonal antibodies are wide. Monoclonal antibody immunotherapy treatment of cancer is highly dependent on the ability to cultivate and purify the critical proteins (Zahavi and Weiner 2020). Production can be done using batch column chromatography. Monoclonal antibodies are recovered by filtering fermentation broth, a mixture of organic and non-organic compounds, from impurities which arise from the protein cultivation. Optimizing production yield and avoiding waste of valuable inputs is both chal- lenging and prohibitively expensive (Liu et al. 2010). The many parameters arising from the complexity of interactions in the production chain, together with many input

components, create unpredictable outputs. If inputs and outputs can be adjusted continuously during the production process it is possible yields can be improved and value created (Meyer et al. 2020).

To build knowledge and predictability of the process, empirical or mechanistic modeling can be performed.

Mechanistic modeling aims to describe the physical and chemical interactions of the components involved. A mechanistic model must be able to replicate the effects of the mass transport of solute molecules, the mobile phase, through a column of packed porous particles, the solid phase. Several phenomena including diffusion, dispersion, adsorption-desorption kinetics, influence the mass transfer. Simultaneously interactions between the phases take place due to differences in hydrophobicity, polarity or net- charge depending on the phases used. In the literature several models of liquid chromatography have been proposed and they can be categorized into equilibrium theory, plate theory and rate models. Out of the existing mechanical models used, the general rate model (GRM), is the most comprehensive. It is not possible to derive a closed form solution to the model, hence numerical methods must be used (Shekhawat and Rathore 2019). Numerical methods for solving such complex problems often have to make a trade-off between computation time and accuracy. The purpose of the present work is to build a software simulation testbed which includes implementations of rate models of liquid chromatography evaluated using the computationally eﬀicient method, orthogonal collocation.

2 The General rate model

The GRM considers all contributions to mass transfer simultaneously in a series of partial differential equations (PDEs). To describe the model, consider a mathematical representation of a cylindrical chromatography column packed with porous resin particles. The length of the column is denoted by 𝐿 in the direction of the axial coordinate 𝑧. The packed particles in the column are represented by an arbitrary number of homogeneous spheres distributed along the length of the column in a single row. Each particles radial coordinate is represented by 𝑟 in a spherical geometry. The boundaries of 𝑧 define inflow and outflow points for the mobile phase in the column. Each particle has flow in the outer boundary of 𝑟 from the column due to spherical symmetry. Figure 1 gives a visual representation of the column and particles. To describe the mass transport variations in time, denoted by 𝑡, mass balance equations are used. Due to the mobile phase being stagnant within the particles it must be calculated separately from the mobile phase outside the particles. Hence two mass balance equations are needed. For the mobile phase in the column, the mass balance is given as:

𝑢𝜕𝐶

𝜕𝑧 +𝜕𝐶

𝜕𝑡 + 𝐹𝜕 ̃𝑞

𝜕𝑡 = 𝐷_𝐿𝜕²𝐶

𝜕𝑧². (2.1) In the equation, 𝑢 is the average flow velocity, 𝐶(𝑡, 𝑧) is the mobile phase solute concentration in the column inter- stitials, 𝐹 = ^1−𝜖_𝜖 is the phase ratio, 𝜖 is the total column porosity, 𝐷_𝐿 is the axial dispersion coeﬀicient. To model the interaction of mobile phase with the stationary phase, 𝑞(𝑡, 𝑧, 𝑟) is defined as the stationary phase concentration and its average by ̃𝑞. Choosing an arbitrary porous particle

1

(6)

z L ^r

Film mass transfer

Porous particle Chromatography column

Inlet Outlet

u

Sorption

Diffusion

Dispersion

Figure 1: The left part of the figure shows the modeled chromatography column with the average mobile phase flow rate 𝑢 and column length 𝐿. The right part shows a close up of one of the porous particles packed in the column. The film mass transfer describes the exchange of concentration between the stagnant fluid inside the porous particles and the free mobile phase in the column. Sorption takes places between the stagnant fluid and the stationary phase coating in the particle. In the column uneven flow causes dispersion effects and inside the particles random molecular motions cause diffusion.

in the column the average stationary phase concentration within that particle can then be defined as:

̃

𝑞 = 3 𝑅_𝑝³∫

𝑅_𝑝

0

𝑟²𝑞 𝑑𝑟. (2.2)

The radius of the porous particle is given by 𝑅_𝑝and from the equation it is implied that ^{𝜕 ̃}_𝜕𝑡^𝑞 denotes the rate of adsorption averaged over the particle.

As the mobile phase transports through the column it encounters the packed spherical particles. The contact between mobile phase and particles will give rise to the effects of film mass transfer. This mass transfer will allow the exchange of mobile phase components in the column with stagnant mobile phase inside the pores of the particles. The mass balance equation for the stagnant mobile phase within the pores of the particles is given by:

𝜖_𝑃𝜕𝐶_𝑝

𝜕𝑡 + (1 − 𝜖_𝑃)𝜕 ̃𝑞

𝜕𝑡 = 𝐷_𝑝(𝜕²𝐶_𝑝

𝜕𝑟² + 2 𝑟

𝜕𝐶_𝑝

𝜕𝑟 ) , (2.3) where 𝜖_𝑃 is the porosity of a particle, 𝐶_𝑝(𝑡, 𝑧, 𝑟) is the concentration of mobile phase solute inside the pores and 𝐷_𝑝 is the diffusion coeﬀicient of solute in particle pores. At the boundary between the particle surface and the interstitial volume, the mass flux is conserved according to:

𝜕 ̃𝑞

𝜕𝑡 = 𝐷_𝑝𝜕𝐶_𝑝

𝜕𝑟 ∣

𝑟=𝑅𝑝

= 𝑘_𝑓(𝐶 − 𝐶_𝑝|_{𝑟=𝑅𝑝}), (2.4)

where 𝑘_𝑓 is the external mass transfer coeﬀicient. The natural boundary condition for the center of a particle is given by:

𝜕𝐶_𝑝

𝜕𝑟 ∣

𝑟=0

= 0 (2.5)

due to spherical symmetry. The concentration of the stationary phase solute, is related to the mobile phase solute in the pores, 𝐶_𝑝 by substituting equation (2.4) into equation (2.1). Hence equation (2.1) becomes:

𝑢𝜕𝐶

𝜕𝑧 +𝜕𝐶

𝜕𝑡 + 𝐹 3

𝑅_𝑝𝑘_𝑓(𝐶 − 𝐶_𝑝|_𝑟=𝑅_𝑝) = 𝐷_𝐿𝜕²𝐶

𝜕𝑧². (2.6) As phases interact, the complexity and size of the protein molecules lead to variation in rates of adsorption and desorption. Fast adsorption between mobile phase solute

concentrations within the porous particles are related to the stationary phase through the equilibrium isotherm.

Assuming infinitely fast interactions, a linear adsorption isotherm can be stated as:

𝑞 = 𝐾_𝑎𝐶_𝑝, (2.7)

where 𝐾_𝑎is the adsorption equilibrium constant. Assum- ing slow adsorption a first order kinetic equation offers the following relation:

𝜕𝑞

𝜕𝑡 = 𝑘_𝑎(𝐶_𝑝− 𝑞

𝐾_𝑎) , (2.8)

where 𝑘_𝑎 is the adsorption rate constant. Second order slow adsorption kinetics are given by:

𝜕𝑞

𝜕𝑡 = 𝑘_𝑎𝐶_𝑝(𝑞_𝑚− 𝑞) − 𝑘_𝑑𝑞, (2.9) where 𝑞_𝑚 is the maximum adsorption capacity, 𝑘_𝑑 is the desorption rate constant (Guiochon et al. 2006). The equations must satisfy the boundary conditions at the column inlet and outlet. First type, or Dirichlet conditions, are given by:

𝐶|_𝑧=0 = 𝐶₀, (2.10)

𝐶|_𝑧=𝐿 = 0, (2.11)

where 𝐶₀denotes the concentration of the solute injected at the inlet. Alternatively using Danckwerts and Neumann boundary conditions allows concentration to fluctuate at the column edges. This reflects a more physical approximation of concentration at the boundaries and within the column. The Danckwerts condition adds the mass flux to the inlet as,

𝐶|_𝑧=0= 𝐶₀+ 𝐷_𝐿 𝑢

𝜕𝐶

𝜕𝑧. (2.12)

At the outlet the Neumann boundary condition is given as:

𝜕𝐶

𝜕𝑧∣

𝑧=𝐿

= 0. (2.13)

Initial conditions for a time dependent, constant rate,

(7)

pulse injection profile is given as:

𝐶(0, 𝑧) = 0, (2.14)

𝐶_𝑝(0, 𝑧, 𝑟) = 0, (2.15)

𝐶(𝑡, 0) = 𝐶₀for 0 ≤ 𝑡 ≤ 𝑡_inj, (2.16) 𝐶(𝑡, 0) = 0 for 𝑡_inj≤ 𝑡, (2.17)

𝑞(0, 𝑧, 𝑟) = 0. (2.18)

Where 𝑡_inj is some arbitrary time under which mobile phase is injected into the column. To adapt the equations to the case of several solute components the variables 𝐶, 𝐶_𝑝, 𝑞 are indexed as 𝐶_𝑛, 𝐶_(𝑝,𝑛), 𝑞_𝑛for 𝑛 = 1, 2, ..., 𝑁_𝑐, where 𝑁_𝑐 is the number of components in the system. This has been omitted in the derivation for the sake of clarity.

2.1 Lumped kinetic models

A family of simplified models, referred to here as lumped kinetic models (LKMs), can be derived from the GRM.

LKMs assume that the column is packed homogeneously, isothermal, neglects radial gradients, under a constant volumetric flow, by lumping internal and external mass transfer and constant dispersion for all components. The section will define five similar models with some common equations. The models are; the transport dispersive model (TDM), transport model, reactive dispersive model (RDM), Thomas model and the equilibrium dispersive model (EDM). The models discard the modeling of the spherical particle axis 𝑟, hence only mobile phase concentration 𝐶(𝑡, 𝑧) and stationary phase concentration 𝑞(𝑡, 𝑧) are used.

The TDM is defined by the following mass balance equation:

𝑢𝜕𝐶

𝜕𝑧 +𝜕𝐶

𝜕𝑡 + 𝐹𝜕𝑞

𝜕𝑡 = 𝐷_𝐿𝜕²𝐶

𝜕𝑧². (2.19) Together with first order adsorption desorption kinetics given by:

𝜕𝑞

𝜕𝑡 = 𝑘_𝑚(𝑞^∗− 𝑞) . (2.20) Here 𝑘_𝑚 is the lumped mass transfer coeﬀicient. The isotherm 𝑞^∗ is the equilibrium relationship between the concentrations of components in the phases. In the model an approximate linear isotherm for small concentrations is given by:

𝑞^∗= 𝑎𝐶, (2.21)

where 𝑎 is the Henry coeﬀicient.

The RDM is defined by combining equation (2.19) with second order adsorption desorption kinetics. For LKMs given by:

𝜕𝑞

𝜕𝑡 = 𝑘_𝑎𝐶(𝑞_𝑚− 𝑞) − 𝑘_𝑑𝑞. (2.22) Two more models are derived by setting the dispersion constant 𝐷_𝐿 to zero. This reduces the mass balance equation to:

𝜕𝐶

𝜕𝑡 + 𝐹𝜕𝑞

𝜕𝑡+ 𝑢𝜕𝐶

𝜕𝑧 = 0. (2.23)

By combining equation (2.23) with equation (2.20) the transport model is obtained. Finally equations (2.23) and (2.22) make up the Thomas model (Golshan-Shirazi and Guiochon 1992).

A limiting case of the LKMs, the EDM, is based on the assumption that equilibrium between mobile and stationary phases is achieved instantaneously. The model uses the mass balance equation,

𝜕𝐶

𝜕𝑡 + 𝐹𝜕𝑞

𝜕𝑡 + 𝑢𝜕𝐶

𝜕𝑧 = 𝐷_𝐴𝜕²𝐶

𝜕𝑧², (2.24) Where 𝐷_𝐴 is the lumped apparent dispersion. The stationary phase concentration is given by equation (2.22).

However, since the assumed instant equilibrium implies that ^𝜕𝑞_𝜕𝑡 = ^𝜕𝑞_𝜕𝑡^∗, equation (2.22) becomes:

𝑞 = 𝑞_𝑚𝐾_𝑅𝐶

1 + 𝐾_𝑅𝐶. (2.25)

where the ratio 𝐾_𝑅= 𝑘_𝑎 𝑘_𝑑.

For the LKMs, Dirichlet boundary conditions are given by:

𝐶|_𝑧=0= 𝐶₀, (2.26)

𝐶|_𝑧=𝐿= 0. (2.27)

Inlet Danckwerts and outlet Neumann boundary conditions are given as:

𝐶|_𝑧=0= 𝐶₀+ 𝐷_𝐿 𝑢

𝜕𝐶

𝜕𝑧, (2.28)

𝜕𝐶

𝜕𝑧|_𝑧=𝐿= 0. (2.29)

Initial conditions for a time dependent pulse injection profile is given as:

𝐶(0, 𝑧) = 0, (2.30)

𝐶(𝑡, 0) = 𝐶₀for 0 ≤ 𝑡 ≤ 𝑡_inj, (2.31) 𝐶(𝑡, 0) = 0 for 𝑡_inj≤ 𝑡, (2.32)

𝑞(0, 𝑧) = 0. (2.33)

2.2 Analytical solutions

Single component analytical solutions in the Laplace domain have been presented by Javeed et al. (2013) for the EDM and LKMs with linear isotherms. For an approximate solution of the TDM the Laplace domain solution is given as:

𝐶(𝑧, 𝑠) = 𝐴𝑒^𝜆¹^𝑧+ 𝐵𝑒^𝜆²^𝑧. (2.34) Where 𝜆_1,2 are given by,

𝜆_1,2= 𝑃 𝑒 2 ∓

√√

√

⎷ (𝑃 𝑒

2 )

2

+ 𝑃 𝑒𝐿 𝑢𝑠 − 𝐿²

𝜖𝐷 𝑘²𝑎 1 − 𝜖 𝑠 + 𝑘 1 − 𝜖

+ 𝐿² 𝜖𝐷𝑎𝑘.

(2.35) The Peclet number, 𝑃 𝑒, follows the relation 𝑃 𝑒 = ^𝐿𝑢_𝐷 and the mass transfer coeﬀicient 𝑘, is related to the lumped mass transfer coeﬀicient by 𝑘_𝑚= _1−𝜖^𝑘 . 𝐴 and 𝐵 for Danck- werts boundary condition are given by:

𝐴 = 𝑐₀𝜆₂𝑒^𝜆² 𝑠 (1 − 𝜆₁

𝑃 𝑒) 𝜆₂𝑒^𝜆²− (1 − 𝜆₂ 𝑃 𝑒) 𝜆₁𝑒^𝜆¹

, (2.36)

𝐵 = −𝑐₀𝜆₁𝑒^𝜆¹ 𝑠 (1 − 𝜆₁

𝑃 𝑒) 𝜆₂𝑒^𝜆²− (1 − 𝜆₂ 𝑃 𝑒) 𝜆₁𝑒^𝜆¹

. (2.37)

The transformation of equation (2.34) back to the time domain can be done using numerical methods. The method proposed by Hoog, Knight, and Stokes (1982), is

3

(8)

used to solve equation (2.34) in the present work. To com- pare results with reference solutions the following equation is used for normalization:

𝜇₁= 𝐿

𝑢(1 + 𝑎𝐹 ). (2.38)

When calculating errors for numerical methods it is beneficial to use the vectors that span the entire solution space instead of using averaging. Large single errors may be- come unnoticeable when using an average based method.

The error norm better reflects the magnitude of individual errors due to the use of the absolute value for every vector element. The 𝐿¹error norm is given by:

𝑁𝑇

∑

𝑖=0

|𝐶_𝑅^𝑖 − 𝐶_𝑁^𝑖|Δ𝑡, (2.39)

where the analytic concentration at any one point is denoted by 𝐶_𝑅^𝑖 and the numeric by 𝐶_𝑁^𝑖, for the time step 𝑖.

Δ𝑡 is the time step length and 𝑁_𝑇 is the number of time steps. The 𝐿² error norm is given by:

𝑁𝑇

∑

𝑖=0

√|𝐶_𝑅^𝑖 − 𝐶_𝑁^𝑖|²Δ𝑡. (2.40)

and the relative error:

∑^𝑁_𝑖=0^𝑇|𝐶_𝑅^𝑖 − 𝐶_𝑁^𝑖|

∑^𝑁_𝑖=0^𝑇|𝐶_𝑅^𝑖| Δ𝑡. (2.41)

3 Discretization

Even if closed form solutions to PDEs can be found, they often concern severely limiting cases. Numerical methods are commonly used to solve PDEs with any practical value. The process of numerically approximating a PDE usually involves discretizing the domain of the problem and creating a solvable system of equations by satisfying boundary conditions. Common methods include the finite difference, finite volume, finite element methods and the methods of weighted residuals. This section presents an approximation scheme using a method of weighted residuals (MWRs) known as orthogonal collocation (OC). The method relies on spatially discretizing the domain in order to reduce the problem to a system of ordinary differential equations (ODEs). The resulting system of ODEs is then solvable by one of the methods proposed in Section 3.2. Finally the section ends by presenting the application of these discretization methods on the LKMs derived in Section 2.1.

3.1 The methods of weighted residuals

MWRs approximate solutions to PDEs using trial functions with a finite number of adjustable parameters. The parameters in these methods are found by setting the weighted average of an equation such that the error, known as the residual in this context, is minimized. In addition to OC the subdomain, least squares, moments and Galerkin methods all belong to the MWRs. The methods are defined by the choice of a weighting function used in the

averaging of the residual. To exemplify consider a differential equation,

𝑑²𝑢

𝑑𝑥² + 𝑔 (𝑥, 𝑢,𝑑𝑢

𝑑𝑥) = 0. (3.1)

where 𝑥 defines a spatial variable, 𝑔(𝑥, 𝑢,^𝑑𝑢_𝑑𝑥) an arbitrary function. Then 𝑢 is the sought solution to the differential equation for 𝑥 ∈ [0, 1], 𝑢(0) = 𝑢₀ and 𝑢(1) = 𝑢₁. The unknown solution 𝑢 can be approximated by a linear combination of trial functions ̃𝑢.

𝑢(𝑥) ≈ ̃𝑢(𝑥) =

𝑛

∑

𝑖=1

𝑎_𝑖𝜓_𝑖(𝑥). (3.2)

In the equation 𝑎_𝑖denotes constant coeﬀicients and 𝜓_𝑖(𝑥) the trial functions. The functions can be based on simple monomials or polynomials. For instance, here the trial functions 𝑃_𝑘represent 𝑘^th degree orthogonal polynomials,

̃

𝑢(𝑥) =

𝑛

∑

𝑘=1

𝑎_𝑘𝑃_𝑘−1(𝑥). (3.3)

To form an expression for the residual 𝑅, equation (3.2) is substituted into (3.1),

𝑅(𝑥, 𝑎_𝑖) =

𝑛

∑

𝑖=1

𝑎_𝑖𝑑²𝜓_𝑖

𝑑𝑥² + 𝑔 (𝑥, ̃𝑢,𝑑 ̃𝑢

𝑑𝑥) ≠ 0. (3.4) As the residual approaches zero the approximation becomes more accurate (Finlayson 1972). To solve for the coeﬀicients 𝑎_𝑖, 𝑛 weighting functions are needed. Then the average residual can be minimized over the domain 𝐷 using the following integral,

∫

𝐷

𝑤_𝑗(𝑥)𝑅(𝑥, 𝑎_𝑖)𝑑𝑥 = 0 for 𝑗 = 1, ..., 𝑛. (3.5) Using equation (3.4) as the residual gives:

∫

1

0

𝑤_𝑗𝑅(𝑥, 𝑎_𝑖)𝑑𝑥 =

𝑛

∑

𝑖=1

𝑎_𝑖∫

1

0

𝑤_𝑗𝑑²𝜓_𝑖

𝑑𝑥² + 𝑤_𝑗𝑔 (𝑥, ̃𝑢,𝑑 ̃𝑢 𝑑𝑥) 𝑑𝑥.

(3.6) Applying the weighting functions will result in system of equations. The system can then be solved for the unknown coeﬀicients 𝑎_𝑖 and the approximate solution ̃𝑢 be found. Depending on the choice of weight functions different systems of equations will be obtained. This paper will focus on the collocation method.

3.1.1 Collocation

The collocation method discretizes a domain into an arbitrary amount of subdomains, where each subdomain is shrunk infinitely thin and can be represented by the Dirac delta function. Hence for collocation the weights in equation (3.5) are given by:

𝑤_𝑗= 𝛿(𝑥 − 𝑥_𝑗) (3.7) for 𝑛 subdomains, or collocation points. A consequence of using the Dirac delta function is that the residual is zero at the chosen collocation points. The weighting function used in collocation doesn’t require integration of the trial functions, making it easier to implement. When choosing the collocation points the spacing between the points can be of importance as shown by Schlömilch et al. (1901).

(9)

When interpolating a function equidistant points will lead to increasing magnitude of oscillations at the boundaries.

Though the application of interpolation is different when approximating a PDE, the effect can cause divergence when increasing the amount of points. To negate these effects Michelsen and J. Villadsen (1972) proposed the roots of orthogonal polynomials be used as collocation points.

3.1.2 Orthogonal Collocation

When the collocation points of the approximation are chosen to be the roots of orthogonal polynomials the method is known as OC. Early works on the subject include pub- lications by Ames (1973), J.V. Villadsen (1995) and more recently thoroughly covered by Young (2019). The collocation points are all given by polynomials belonging to the set of Jacobi polynomials orthogonal on [-1,1] as defined by:

∫

1

−1

(1−𝑥)^𝛼(1−𝑥)^𝛽𝑃𝑛^(𝛼,𝛽)𝑃𝑚^(𝛼,𝛽)(𝑥)𝑑𝑥 = 0, for 𝑛 ≠ 𝑚. (3.8)

In the equation 𝑃_𝑚, 𝑃_𝑛denote polynomials of degree 𝑚, 𝑛, 𝛼 and 𝛽 determine which Jacobi polynomial is defined.

Setting 𝛼 = 𝛽 = −¹₂ gives 1^𝑠𝑡 kind Chebyshev points, 𝛼 = 𝛽 = 0 gives Legendre points, 𝛼 = 𝛽 = ¹₂ gives 2^𝑛𝑑 kind Chebyshev points, 𝛼 = 𝛽 = 1 gives Lobatto points.

Radau points are asymmetric and given by 𝛼 = 0, 𝛽 = 1 or 𝛼 = 1, 𝛽 = 0. In OC the domain is usually shifted to [0,1]

which makes scaling along a spatial axis more convenient.

By deconstructing the polynomial 𝑃_𝑛into 𝑃_𝑛= 𝜌_𝑛𝑝_𝑛, the leading coeﬀicients 𝜌_𝑛 are given by,

𝜌_𝑛= Γ(2𝑛 + 𝛼 + 𝛽 + 1) 2^𝑛𝑛!Γ(𝑛 + 𝛼 + 𝛽 + 1)

where Γ(𝑛) = (𝑛 − 1)! and 𝑝_𝑛 is given by the recurrence relation for monomials,

𝑝_𝑛+1= (𝑥 − ̂𝛼)𝑝_𝑛− ̂𝛽𝑝_𝑛−1. (3.9) The coeﬀicients are given by,

̂

𝛼_𝑛= 𝛽²− 𝛼²

(2𝑛 + 𝛼 + 𝛽)(2𝑛 + 𝛼 + 𝛽 + 2),

̂𝛽_𝑛= 4𝑛(𝑛 + 𝛼)(𝑛 + 𝛽)(𝑛 + 𝛼 + 𝛽) (2𝑛 + 𝛼 + 𝛽)²[(2𝑛 + 𝛼 + 𝛽)²− 1].

To obtain the roots of the chosen Jacobi polynomial numerical differentiation using a Newton-Raphson method can be performed. The method uses a recurrence formula derived from equation (3.8) and finds each sought root by iteration to some arbitrary tolerance level. The method finds the 𝑘^throot by,

𝑥^𝑖+1_𝑘 = 𝑥^𝑚_𝑘 −𝑝_𝑛(𝑥^𝑖_𝑘)

𝑝^′_𝑛(𝑥^𝑖_𝑘), (3.10) where 𝑖 indicates iteration, 𝑝_𝑛 a polynomial given by the recurrence relation in equation (3.9) and 𝑝^′_𝑛 its derivative. When constructing trial functions with polynomials a modal approach can be taken as in equation (3.3). It is more convenient to restate into a nodal approximation by letting the parameters represent the approximate solution

at the collocation points 𝑎_𝑖= 𝑢(𝑥_𝑖), instead. The resulting trial solution in equation (3.2) becomes,

̃

𝑢(𝑥) =

𝑛

∑

𝑖=1

𝑢(𝑥_𝑖)𝑙_𝑖(𝑥). (3.11) In this case Lagrange interpolating polynomials, 𝑙_𝑖, are used as trial functions. The polynomials are given by:

𝑙_𝑖(𝑥) =

𝑛+1

∏

𝑗=0 𝑗≠𝑖

𝑥 − 𝑥_𝑗

𝑥_𝑖− 𝑥_𝑗. (3.12) Here 𝑥₀ and 𝑥_𝑛+1 are added boundary points at 0 and 1. The interior points belong to the roots of orthogonal polynomials derived from equation (3.8). The boundary points are appended as they are not included in the set of Jacobi polynomial roots on the interval [0, 1]. They are however required to deal with boundary conditions im- posed on PDEs. The final nodal trial solution is given as:

𝑢(𝑥) ≈

𝑛+1

∑

𝑖=0

𝑢(𝑥_𝑖)𝑙_𝑖(𝑥). (3.13)

3.1.3 Differentiation matrices

To approximate derivatives in terms of the nodal approximations from Section 3.1.2 the trial functions in equation (3.13) are differentiated. Taking the derivative of the log- arithm of equation (3.12) gives a sum for the diagonal:

𝐴_𝑖𝑗= 𝑑𝑙_𝑗(𝑥) 𝑑𝑥 ∣

𝑥𝑖

=

⎧{ {

⎨{ {⎩

𝑛+1

∑

𝑘=0𝑘≠𝑖

1

𝑥_𝑖− 𝑥_𝑘 for 𝑖 = 𝑗, and

̂

𝑝^′_𝑛(𝑥_𝑖)

(𝑥_𝑖− 𝑥_𝑗) ̂𝑝_𝑛^′(𝑥_𝑗) for 𝑖 ≠ 𝑗.

(3.14)

Here ̂𝑝_𝑛^′ in the non-diagonal elements is given by:

̂

𝑝^′_𝑛(𝑥_𝑖) =

𝑛+1

∏

𝑗=0 𝑗≠𝑖

(𝑥_𝑖− 𝑥_𝑗), (3.15)

where the power of the polynomial corresponds to 𝑛 + 1.

The second derivative is then be calculated by matrix mul- tiplying the first derivative:

𝐵_𝑖𝑗= 𝑑²𝑙_𝑗(𝑥) 𝑑𝑥² ∣

𝑥𝑖

=

𝑛+1

∑

𝑘=0

𝐴_𝑖𝑘𝐴_𝑘𝑗. (3.16)

3.1.4 Orthogonal collocation on finite elements

An extended implementation of OC is known as orthogonal collocation over finite elements (OCFE). The method applies OC with 𝑛 points for every 𝑛_𝑒elements partition- ing the domain. This is beneficial because the amount of grid points can be increased without the need to spend resources to calculate increasingly higher order polynomial roots and differential matrices. Collocation points are mapped to each element using the equation:

𝑔^𝑙(𝑥) = (𝑥 − 𝑥_𝑙)

Δ𝑥_𝑙 , (3.17)

where Δ𝑥_𝑙= 𝑥_𝑙+1− 𝑥_𝑙 and 𝑙 = 1, ..., 𝑛_𝑒elements. Differen- tial operators must be scaled by Δ𝑥_𝑙. The first differential matrix becomes:

𝐴^𝑙_𝑖𝑗= 𝐴_𝑖𝑗

Δ𝑥_𝑙, (3.18)

5

(10)

and the second:

𝐵^𝑙_𝑖𝑗= 𝐵_𝑖𝑗

Δ𝑥²_𝑙. (3.19)

The interior element boundaries requires continuity of the first derivative:

1 Δ𝑥_𝑙

𝑛+1

∑

𝑖=0

𝐴_𝑛+1,𝑖𝑦_𝑖^𝑙− 1 Δ𝑥_𝑙+1

𝑛+1

∑

𝑖=0

𝐴_0,𝑛+1𝑦^𝑙+1_𝑖 = 0, (3.20)

for 𝑙 = 1, ..., 𝑛_𝑒− 1 elements and 𝑦^𝑙_𝑖= 𝑦(𝑔_𝑖^𝑙) = 𝑦(𝑥_𝑙+ 𝑔^𝑙_𝑖Δ𝑥_𝑙) (Carey and Finlayson 1975).

3.2 Systems of ordinary differential equations Once spatially discretized using OC, a time dependent PDE as presented in equation (2.19) is reduced into a system of ODEs. This section introduces a selection of methods that can be used to solve such a system. The methods assume that the ODE is stated as an initial value problem on the form:

𝑦^′(𝑡) = 𝑓(𝑡, 𝑦(𝑡)), (3.21) where 𝑡 is time, 𝑓 represents some function applying the ODE. By letting the dependent variable 𝑦, and 𝑡 take discrete steps the initial step can be stated as 𝑦(𝑡₀) = 𝑦₀. Then by using the forward difference algorithm to approximate the derivative a numerical scheme can be derived:

𝑦^′(𝑡) ≈ 𝑦(𝑡_𝑖+1) − 𝑦(𝑡_𝑖)

ℎ = 𝑓(𝑡, 𝑦(𝑡)). (3.22) Here 𝑡_𝑖+1is given by 𝑡_𝑖+1= 𝑡_𝑖+ ℎ, where ℎ is some discrete time step. A basic numerical method for approximating ODEs, the Euler method, is given by rearranging the equation as:

𝑦(𝑡_𝑖+1) = 𝑦(𝑡_𝑖) + ℎ𝑓(𝑡, 𝑦(𝑡)), (3.23) The iteration for the then done over the consecutive time steps 𝑖 = 1, ..., 𝑡_max, where 𝑡_max is the final step (Landau et al. 2015).

3.2.1 The implicit Runge-Kutta Radau IIA method To solve ODEs the implicit Runge-Kutta (RK) family of methods are commonly used due to their stability properties. The RK methods are multistep methods and are distinguished by varying several coeﬀicients. The implicit RK methods can be generalized by:

𝑦_𝑖+1= 𝑦_𝑖+ ℎ

𝑠

∑

𝑗=1

𝑏_𝑗𝑘_𝑗, (3.24)

where the 𝑘_𝑗 is given by:

𝑘_𝑗= 𝑓 (𝑡_𝑖+ 𝑐_𝑗ℎ, 𝑦_𝑖+ ℎ

𝑠

∑

𝑘=1

𝑎_𝑗𝑘𝑘_𝑗) , for 𝑗 = 1, ..., 𝑠.

(3.25) The here 𝑠 denotes the amount of stages and the coef- ficients 𝑎_(𝑗,𝑘), 𝑐_𝑗, 𝑏_𝑗 define the RK method. The coeﬀi- cients are usually presented in what is known as a Butcher tableau. A Butcher tableau follows the notation in presented in table 1.

Table 1: The table shows the relation between equations (3.24), (3.25) and Butcher tableaus in tables 2, 3, 4, describing the implementation of different RK methods.

𝑐₁ 𝑎_1,1 ⋯ 𝑎_1,𝑠 𝑐₂ 𝑎_2,1 ⋯ 𝑎_2,𝑠

⋮ ⋮ ⋱ ⋮

𝑐_𝑠 𝑎_𝑠,1 ⋯ 𝑎_𝑠,𝑠 𝑏₁ ⋯ 𝑏_𝑠

An implementation of an implicit RK method, the Radau IIA method, is based on the Radau quadrature formulas.

The 5th order Radau IIA method is given by equations (3.24), (3.25) with the constants given by the Butcher tableau in table 2 (Hairer and Wanner 1996).

Table 2: The Butcher tableau defining the con- stants for equations (3.24) and (3.25) to construct the 5th order Radau IIA method.

4−√ 6 10

88−7√ 6 360

296−169√ 6 1800

−2+3√ 6 225 4+√

6 10

296+169√ 6 1800

88+7√ 6 360

−2−3√ 6 225

1 ¹⁶⁻₃₆^√⁶ ¹⁶⁺₃₆^√⁶ ¹₉

16−√ 6

36 16+√

6

36 1

9

3.2.2 Improved Euler

The improved Euler method is a two-stage RK method based on the Euler method. The method is given by the following equations:

̂

𝑦_𝑖+1 = 𝑦_𝑖+ ℎ𝑓(𝑡_𝑖, 𝑦_𝑖), (3.26) 𝑦_𝑖+1 = 𝑦_𝑖+ℎ

2[𝑓(𝑡_𝑖+𝑖, ̂𝑦_𝑖+1) + 𝑓(𝑡_𝑖, 𝑦_𝑖)]. (3.27) Here ̂𝑦_𝑖+1 is evaluated at the same time step 𝑦_𝑖+1(Kong, Siauw, and Bayen 2020). As the improved Euler method belongs to the family of RK methods it can also be de- scribed by the Butcher tableau in table 3.

Table 3: The Butcher tableau defining the con- stants for equations (3.24) and (3.25) to construct the improved Euler method.

0 1 1

1 2

3.2.3 Implicit trapezoidal method

The implicit trapezoidal method is a linear multistep method similar to the improved Euler method (Iserles 2008). The method is second order and given by:

𝑦_𝑖+1= 𝑦_𝑖+ℎ

2[𝑓(𝑡_𝑖+1, 𝑦_𝑖+1) + 𝑓(𝑡_𝑖, 𝑦_𝑖)] , (3.28) or in the form of the Butcher tableau in table 4.

(11)

Table 4: The Butcher tableau defining the con- stants for equations (3.24) and (3.25) to construct the implicit trapezoidal method.

0 1 ¹₂ ¹₂

1 2

3.2.4 Backward differentiation formula

The Backward differentiation formula (BDF) methods belong to a different family than the RK methods. Like the implicit RK the BDFs are suited for stiff problems. A generalized BDF method with a constant step size can be stated as,

𝑘

∑

𝑚=1

1

𝑚∇^𝑚𝑦_𝑖+1− ℎ𝑓(𝑡_𝑖+1, 𝑦_𝑖+1) = 0, (3.29) for order 𝑘, and solved for 𝑦_𝑖+1 to iterate (Shampine and Reichelt 1997).

3.3 Discretized lumped kinetic models

To find approximate solutions to the models in Section 2.1 using OC, trial solutions are formed by spatially discretizing the 𝑧 axis and substituting the derivatives with the differentiation matrices from Section 3.1.3. The resulting system of equations is then solved using one of the ODE schemes presented in Section 3.2. Using equation (3.13) a trial solution for the mobile phase solute concentration 𝐶 can be stated as:

𝐶 =

𝑛+1

∑

𝑖=0

𝑙_𝑖(𝑧)𝐶(𝑧_𝑖), (3.30) where the boundary points 𝑧₀and 𝑧_𝑛+1are included in the interpolation. The interior points are chosen to be the roots of an orthogonal polynomial. The residual is formed by substituting the trial solution into the approximated equation.

3.3.1 Transport dispersive model

Starting with the dispersive mass balance equation (2.19), and substituting the derivatives ^𝜕𝐶_𝜕𝑧,^𝜕_𝜕𝑧²^𝐶2, with the operators in equation (3.14) and equation (3.16) respectively, the following system is obtained:

𝜕𝐶(𝑧_𝑖)

𝜕𝑡 +𝐹𝜕𝑞(𝑧_𝑖)

𝜕𝑡 −

𝑛

∑

𝑖=1

𝐷_𝐿𝐶(𝑧_𝑖)𝜕²𝑙_𝑖

𝜕𝑧²∣

𝑧𝑗

+𝑢𝐶(𝑧_𝑖)𝜕²𝑙_𝑖

𝜕𝑧²∣

𝑧𝑗

= 0.

(3.31) Where 𝐶(𝑧_𝑖) are the nodal values for 𝐶. The collocation method implies that the residual is zero from the weighting equation (3.7) in the interior points 𝑧_𝑗for 𝑗 = 1, 2, ...𝑛.

Letting 𝐶_𝑖= 𝐶(𝑧_𝑖), 𝑞_𝑖= 𝑞(𝑧_𝑖), substituting the derivatives 𝐴_𝑗𝑖= ^𝑑𝑙_𝑑𝑧^𝑖|_𝑧𝑗and 𝐵_𝑗𝑖= ^𝑑_𝑑𝑧²^𝑙2^𝑖|_𝑧𝑗the expression simplifies into:

𝜕𝐶_𝑖

𝜕𝑡 + 𝐹𝜕𝑞_𝑖

𝜕𝑡 −

𝑛

∑

𝑖=1

(𝐷_𝐿𝐵_𝑗𝑖+ 𝑢𝐴_𝑗𝑖)𝐶_𝑖= 0. (3.32) The first order kinetics in equation (2.20) are restated as:

𝜕𝑞_𝑖

𝜕𝑡 = 𝑘_𝑚(𝐶_𝑖𝑎 − 𝑞_𝑖) . (3.33)

The boundary conditions are implemented from equations (2.28) and (2.29). The inlet boundary as:

𝐶∣𝑧=0= 𝐶_inj+𝐷_𝐿 𝑢

𝑛+1

∑

𝑖=0

𝐴_(0,𝑖)𝐶_𝑖, (3.34)

where 𝐶_inj is the injection concentration. The outlet boundary condition is given by:

0 = 𝑢

𝑛+1

∑

𝑖=0

𝐴_{(𝑛+1,𝑖)}𝐶_𝑖. (3.35)

3.3.2 Reactive dispersive model

The RDM combines the dispersive mass balance equation (2.19), with second order kinetics stated in equation (2.22).

The first equation converts into equation (3.32) as before and the kinetic equation is given by:

𝜕𝑞_𝑖

𝜕𝑡 = 𝑘_𝑎𝐶_𝑖(𝑞_𝑚− 𝑞_𝑖) + 𝑘_𝑑𝑞_𝑖. (3.36) The boundary conditions are the same as for the TDM and are given by equations (3.34) and (3.35).

3.3.3 Transport model

The transport model uses the non-dispersive mass balance equation (2.23). By substituting the derivatives as before, the following system is obtained:

𝜕𝐶_𝑖

𝜕𝑡 −

𝑛

∑

𝑖=1

𝑢𝐴_𝑗𝑖𝐶_𝑖= 0. (3.37)

The model uses the first order kinetics in equation (3.33).

The Danckwerts inlet boundary condition without dispersion constant is used:

𝐶|_𝑧=0= 𝐶_inj+ 1 𝑢

𝑛+1

∑

𝑖=0

𝐴_(0,𝑖)𝐶_𝑖 (3.38)

and the Neumann outlet condition is given by equation (3.35).

3.3.4 Thomas model

The Thomas model is constructed using the converted non-dispersive mass balance equation (3.37), second order kinetics given by equation (3.36) and the boundary conditions from equations (3.38), (3.35).

3.3.5 Equilibrium dispersive model

To create the EDM collocation is applied to the dispersive mass balance equation (2.24):

𝜕𝐶_𝑖

𝜕𝑡 −

𝑛

∑

𝑖=1

(𝐷_𝐴𝐵_𝑗𝑖+ 𝑢𝐴_𝑗𝑖)𝐶_𝑖= 0. (3.39)

The kinetic model in equation (2.25) is adapted as:

𝑞_𝑖= 𝑞_𝑚𝐾_𝑟𝐶_𝑖

1 + 𝐾_𝑟𝐶_𝑖. (3.40) Finally the same boundary conditions are used as for previous dispersive models given by equations (3.34) and (3.35).

7

(12)

3.4 Injection profiles

The equations in this section are used to define variations in component concentrations at the inlet of a modeled chromatography column. The equations are referred to as injection profiles in this text. A regular step function is defined as:

𝑓_step(𝑥) =

⎧{

⎨{

⎩

0 𝑥 ≤ 𝑥_𝑒₁, 1 𝑥_𝑒₁ < 𝑥 < 𝑥_𝑒₂, 0 𝑥_𝑒₂ ≤ 𝑥,

, (3.41)

where 𝑥_𝑒₁and 𝑥_𝑒₂define the boundaries of the step. Ebert et al. (2002), suggests a smooth interpolated step function which can be used here reduce the stiffness of an injection:

𝑓_smooth(𝑥) =

⎧{

⎨{

⎩

0 𝑥 ≤ 0,

6𝑥⁵− 15𝑥⁴+ 10𝑥³ 0 ≤ 𝑥 ≤ 1,

1 1 ≤ 𝑥.

(3.42)

The maximum value of the equation is limited to [0, 1] by the clamping function:

𝑓_clamp(𝑥) =

⎧{

⎨{

⎩

0 if 𝑥 ≤ 0, 𝑥 if 0 < 𝑥 < 1, 1 if 1 ≤ 𝑥.

(3.43)

A Gaussian profile can be obtained by letting 𝑏 define the peak of the curve, and 𝑐 the axial spread:

𝑓_gauss(𝑥) = exp [−(𝑥 − 𝑏)²

2𝑐² ]. (3.44) A linear slope concentration profile is generated by:

𝑓_slope(𝑥) = {

𝑥 𝑘𝑠 if _𝑘^𝑥

𝑠 ≤ 1, 0 if _𝑘^𝑥

𝑠 > 1, (3.45) where the slope angle is determined by the constant 𝑘_𝑠. To generate a reverse slope the previous equation is trans- formed into:

𝑓_revslope(𝑥) = {1 −_𝑘^𝑥

𝑠 if _𝑘^𝑥

𝑠 ≤ 1, 0 if _𝑘^𝑥

𝑠 > 1. (3.46)

4 The stlc package

The GRM is in many ways a much more complicated model than any of the LKMs presented in Section 2.1. The consequence is that the GRM requires more computational power to calculate approximate solutions in the context of liquid chromatography. Investigations by Golshan-Shirazi and Guiochon (1992) have shown that accurate approximations of the GRM can be achieved by the LKMs, ex- cept in the case of extremely small pore diffusivity. The stlc package is proposed as a collection of models with the aim to achieve approximate results of the GRM in a more eﬀicient and accessible manner. To achieve this the package contains discretized versions of the LKMs presented in Section 2.1, together with a collection of ODE solvers, implemented using the Python language. The goal of the package is to create a simulation testbed that delivers approximate solutions that are as close as possible that obtained by the GRM without sacrificing evaluation time.

4.1 Dependencies

To make the package as portable as possible the use of external dependencies has been minimized. However the benefit of using some third-party packages heavily out- weighs the cost. The stlc package relies mainly on the well established matrix library numpy (Harris et al. 2020). The implicit RK and BDF ODE solvers, equation (3.24), (3.29) respectively, are sourced from the scipy package (Virtanen et al. 2020). The analytic solution in equation (2.34) uses the inverse Laplace transform function from the package mpmath (Johansson et al. 2013). The included plotting functions make use of matplotlib (Hunter 2007) and plotly (P. T. Inc. 2015).

4.2 Development

To support the development process there is software that can aid a programmer in various ways. Such tools can automate repetitive tasks, help with adherence to Python ecosystem or software development standards and speed up development. In addition to the use of git as ver- sion control system the software mentioned in this section makes up some of the tooling used in the development of the package.

4.2.1 Unit testing

The use of unit testing to produce and confirm expected results from functions is a way to verify code is working as intended. The more thorough tests are designed the more security in functionality can be expected. A way to quantify the amount of testing done on a program is to measure code coverage. If every part of the code is tested then coverage will reach 100%. Unit tests have been implemented using the Python unit testing framework and coverage checked with coverage.py (Gareth Rees 2021).

4.2.2 Type checking

Python is a dynamically typed language which allows a simple and flexible syntax but loses the benefits of type safety. The Python enhancement proposal (PEP) 484 (Guido van Rossum and Langa 2014), introduced optional type hinting which allows programmers access to some of the benefits of typing if needed. Type hinting can aid in error discovery, code readability, improve linting and add an extra layer of communication in a function decla- ration. The program mypy (Jukka Lehtosalo 2021), intro- duces optional static type checking into Python according to the Python enhancement proposal PEP 484 and has been used throughout the development of the package in combination with type hinting.

4.2.3 Linting and formatting

Linting is a process that continuously checks code for se- mantic and syntax errors as well as code style. A linter can be integrated into the development environment via a text editors language server. The Pyright static type checker has been used for linting throughout the development. It implements type checking according to PEP 484 and additional PEP. There is overlap with mypy but there the tools complement each other. Furthermore continuous

(13)

code formatting has been automated with YAPF (G. Inc.

2021), configured to follow the PEP 8 style guide which defines standards for formatting (Rossum, Warsaw, and Coghlan 2001).

4.2.4 Documentation

To improve the usability of the package code it has been commented using docstrings adhering to the PEP 257 con- vention (David Goodger 2001). Furthermore a configura- tion file for the Sphinx Python documentation generator has been included in docs\source\conf.py. Sphinx allows the automatic generation of documentation in various out- put formats from docstrings in source code (Brandl 2021).

The documentation for the stlc package is included in ap- pendix A.

4.3 Modules

The stlc package is built up of several modules, all pack- age code is located in the src directory. The folder src/stlc contains implemented models from Section 2.1.

The src/occ contains the package authored by Young (2019) which performs the OC discretization using several functions including implementations of equations (3.10), (3.14) and (3.16). The folder src/examples contains several examples which demo the capabilities of the package. Unit tests are contained in the folder src/tests.

The src/utilities folder contains the implemented analytic solution from equation (2.34) in laplace.py.

plot_functions.py implements some basic plot functionality used in examples and some models. The package file tree is visualized in Section 4.3.1.

4.3.1 Package file tree

The schematic in Figure 2 shows the distribution of the package files from the perspective of the package root directory.

root doc

source conf.py src

stlc lkm.py flkm.py grm.py

boundary_conditions.py injection_profiles.py ode.py

ocfe ocfe.py occ

arrayprint.py jacobi.py occ.py utilities

laplace.py plot_functions.py tests

examples

Figure 2: The stlc package file tree. The package documentation can be found in appendix A.

4.3.2 LKM

The module in lkm.py implements the LKMs discretized in Section 3.3. The models are constructed using multi- ple inheritance. Unique classes are defined for each of the two mass balance equations (3.32), (3.37) and the two kinetic equations (3.33), (3.36). A parent class containing common methods and parameters inherits one of the mass balance classes and one kinetic equation class resulting in a complete model. The model can be one of the TDM, Transport model, RDM, Thomas model or the Equilib- rium dispersive models. Upon creation, the model object applies collocation and solves the resulting system of ODEs using one of the integration methods made available.

4.3.3 Functional LKM

The module in flkm.py implements the TDM, Transport model, RDM and Thomas model with a different interface than the lkm module. Here a model is created by instan- tiating a model specific parameter class. The parameter class is a data structure containing the parameters for one of the TDM, Transport model, RDM or Thomas model.

The parameter data structure is passed together with more general options into the create_model function to generate a model object. The model object combines one of the equations from the mass balance equations (3.32), (3.37) and one of the kinetic equations (3.33), (3.36). The model object can be passed into the solve function which also requires the choice of integration method together with inte- gration specific parameters. Compared to the lkm module flkm is inspired by the functional paradigm and may be more appropriate for certain use cases or for other reasons preferred by some users.

4.3.4 Injection profiles

The injection profiles given by equations (3.41) through (3.46) in Section 3.4 have all been implemented in injection_profiles.py. One of these functions can be passed into a model instance as is, or by way of a lambda function if function parameters other than the current time step are required.

4.3.5 Boundary conditions

In boundary_conditions.py, both Dirichlet, equation (2.26), and Danckwerts inlet boundary conditions have been implemented for both the dispersive mass transport equation (3.34), and non-dispersive in equation (3.38).

The outlet boundary conditions can be either Dirichlet, equation (2.27), or Neumann, equation (3.35) which are used when the Danckwerts boundary conditions are passed into a model.

4.3.6 ODE solvers

The functions contained in ode.py implement ODE solvers based on equations (3.24) through (3.29). The improved Euler and trapezoidal methods are implemented natively in the module, whereas the implicit BDF and RK methods are wrapped scipy implementations. A solver is chosen by passing a string to the solve method in the lkm module or the solve function in flkm.

9

Simulation Testbed for Liquid Chromatography