Develop a Matlab code implementing the boundary value design method for simple t

ID: 1023211 • Letter: D

Question

Develop a Matlab code implementing the boundary value design method for simple ternary distillation using bottom-up and top-down tray-by-tray calculations. For simplicity, assume CMO and constant relative volatilities. In that case, the governing equilibrium relation is Using the developed code, design a ternary column with moderate (2) and difficult (1.2) adjacent component relative volatilities effecting a direct and an indirect split. The key component impurities in the product streams are 1 mol% each. Also assume an equimolar feed at its bubble point. Submit a short report (no more than 5 pages) containing: Step-by-step design method algorithm. Ternary column profiles at the design reflux and the minimum. reflux for each design. Matlab code (as an appendix) Based on your experience with the simulations, comment on the effect of the third product stream mol fraction on the column design.

Explanation / Answer

The BVP solver is designed to handle systems of ordinary differential equations

y = f(x, y)

where x is the independent variable, y is the dependent variable, and y represents the derivative of y with respect to x dy/dx.

See Choose an ODE Solver for general information about ODEs.

Boundary Conditions

In a boundary value problem, the solution of interest satisfies certain boundary conditions. These conditions specify a relationship between the values of the solution at more than one x. In its basic syntax, bvp4c is designed to solve two-point BVPs, i.e., problems where the solution sought on an interval [a, b] must satisfy the boundary conditions

g(y(a), y(b)) = 0

Unlike initial value problems, a boundary value problem may not have a solution, may have a finite number of solutions, or may have infinitely many solutions. As an integral part of the process of solving a BVP, you need to provide a guess for the required solution. The quality of this guess can be critical for the solver performance and even for a successful computation.

There may be other difficulties when solving BVPs, such as problems imposed on infinite intervals or problems that involve singular coefficients. Often BVPs involve unknown parameters p that have to be determined as part of solving the problem

y = f(x, y, p)

g(y(a), y(b), p) = 0

In this case, the boundary conditions must suffice to determine the value of p.

The function bvp4c solves two-point boundary value problems for ordinary differential equations (ODEs). It integrates a system of first-order ordinary differential equations

y = f(x, y)

on the interval [a, b], subject to general two-point boundary conditions

bc(y(a),y(b)) = 0

It can also accommodate other types of BVP problems, such as those that have any of the following:

Unknown parameters

Singularities in the solutions

Multipoint conditions

In this case, the number of boundary conditions must be sufficient to determine the solution and the unknown parameters.

bvp4c produces a solution that is continuous on [a,b] and has a continuous first derivative there. You can use the function deval and the output of bvp4c to evaluate the solution at specific points on the interval of integration.

bvp4c is a finite difference code that implements the 3-stage Lobatto IIIa formula. This is a collocation formula and the collocation polynomial provides a C1-continuous solution that is fourth-order accurate uniformly in the interval of integration. Mesh selection and error control are based on the residual of the continuous solution.

The collocation technique uses a mesh of points to divide the interval of integration into subintervals. The solver determines a numerical solution by solving a global system of algebraic equations resulting from the boundary conditions, and the collocation conditions imposed on all the subintervals. The solver then estimates the error of the numerical solution on each subinterval. If the solution does not satisfy the tolerance criteria, the solver adapts the mesh and repeats the process. The user must provide the points of the initial mesh as well as an initial approximation of the solution at the mesh points.

BVP Solver Syntax

The basic syntax of the BVP solver is

The input arguments are

odefun

A function handle that evaluates the differential equations. It has the basic form

where x is a scalar, and dydx and y are column vectors. odefun can also accept a vector of unknown parameters and a variable number of known parameters, (see BVP Solver Options).

bcfun

Handle to a function that evaluates the residual in the boundary conditions. It has the basic form

where ya and yb are column vectors representing y(a) and y(b), and res is a column vector of the residual in satisfying the boundary conditions. bcfuncan also accept a vector of unknown parameters and a variable number of known parameters, (see BVP Solver Options).

solinit

Structure with fields x and y:

Ordered nodes of the initial mesh. Boundary conditions are imposed at a = solinit.x(1) and b = solinit.x(end).

Initial guess for the solution with solinit.y(:,i) a guess for the solution at the node solinit.x(i).

The structure can have any name, but the fields must be named x and y. It can also contain a vector that provides an initial guess for unknown parameters. You can form solinit with the helper function bvpinit. See the bvpinit reference page for details.

The output argument sol is a structure created by the solver. In the basic case the structure has fields x, y, yp, and solver.

sol.x

Nodes of the mesh selected by bvp4c

sol.y

Approximation to y(x) at the mesh points of sol.x

sol.yp

Approximation to y(x) at the mesh points of sol.x

sol.solver

'bvp4c'

The structure sol returned by bvp4c contains an additional field if the problem involves unknown parameters:

sol.parameters

Value of unknown parameters, if present, found by the solver.

The function deval uses the output structure sol to evaluate the numerical solution at any point from [a,b].

BVP Solver Options

For more advanced applications, you can specify solver options by passing an input argument options.

options

Structure of optional parameters that change the default integration properties. This is the fourth input argument.

You can create the structure options using the function bvpset. The bvpset reference page describes the properties you can specify.

Integrator Options

The default integration properties in the BVP solver bvp4c are selected to handle common problems. In some cases, you can improve solver performance by overriding these defaults. You do this by supplying bvp4c with an options structure that specifies one or more property values.

For example, to change the value of the relative error tolerance of bvp4c from the default value of 1e-3 to 1e-4,

Create an options structure using the function bvpset by entering

Pass the options structure to bvp4c as follows:

Code the ODE and boundary condition functions. Code the differential equation and the boundary conditions as functions that bvp4c can use:

The code below represents the differential equation and the boundary conditions in the functions shockODE and shockBC. Note that shockODE is vectorized to improve solver performance. The additional parameter is represented by e and is shared with the outer function.

Provide analytical partial derivatives. For this problem, the solver benefits from using analytical partial derivatives. The code below represents the derivatives in functionsshockJac and shockBCJac.

shockJac shares e with the outer function.

Tell bvp4c to use these functions to evaluate the partial derivatives by setting the options FJacobian and BCJacobian. Also set 'Vectorized' to 'on' to indicate that the differential equation function shockODE is vectorized.

Create an initial guess. You must provide bvp4c with a guess structure that contains an initial mesh and a guess for values of the solution at the mesh points. A constant guess of y(x) 1 and y(x) 0, and a mesh of five equally spaced points on [–1 1] suffice to solve the problem for = 10–2. Use bvpinit to form the guess structure.

Use continuation to solve the problem. To obtain the solution for the parameter = 10–4, the example uses continuation by solving a sequence of problems for = 10–2, 10–3, 10–4. The solver bvp4c does not perform continuation automatically, but the code's user interface has been designed to make continuation easy. The code uses the output solthat bvp4c produces for one value of e as the guess in the next iteration.

View the results. Complete the example by displaying the final solution