[Go to previous section: 1.4| Go
to next section: 1.6]

[ Go to CENG301 Homepage
| CENG301 Notes Table of Contents]

- 1.5.0 Intro
- 1.5.1 The Distillation Function: fraco
- 1.5.2 The Separator Function: sept
- 1.5.3 The Reactor Function: reactr

Page 2

- 1.5.4 The Root Finding "Module": ssec2
- 1.5.5 Extraction of a Solute with phase
- 1.5.6 Extraction of a Solute with: phaseg

The root finding functions **ssec1** and
**ssec2** were developed for a course in computer
applications to illustrate the ease of producing a general purpose
function to find a scalar root of one nonlinear equation. These
functions use the secant method discussed in Reklaitis on page 285 to
find a root in a specified interval. The functions will be found to
be very useful in this course since many of our equations require
trial and error solution for a single unknown.

Here is what help shows for **ssec2**:

>>help ssec2function x=ssec2(xs,fx,c,erc,nmax) finds linear approximation of f(x)=0 xs gives an interval to seek a solution to f(x) = 0 in. fx is a character vector that tells how to find f(x). c is a vector of parameters in f(x) erc gives the difference between successive values of x allowed when convergence is deemed to be achieved. nmax is the maximum number of trials, before quitting. If erc and nmax are not given, then the program uses: erc=.0001 and nmax=10 If only two arguments are given, c is set to null. The program uses the function ssec1 for each iteration. Example 1: >> cs=[1 -2 -4.25 7.5]; >> disp(ssec2([2 3],'polyval(c,x)',cs,.0001,10)) Example 2: >> disp(ssec2([3 4],'sin(x)'))

The program uses the function **ssec1** for each
iteration.

Here is the first example suggested in the help listing for
**ssec2**:

>>cs=[1 -2 -4.25 7.5];>>disp(ssec2([2 3],'polyval(c,x)',cs,.0001,10))2.5000<-- One root of the cubic.>>polyval(cs,2.5)ans = 0<-- It checks!

>>disp(ssec2([3 4],'sin(x)'))3.1416<-- Very close to pi

Here is another example that shows how to find multiple roots with
**ssec2**

Example 3: Find all roots of a quadratic:

f(x) = -0.451*x ^{2}+ 0.2836*x + 0.5641

By using the **ssec2** function:

>>ssec2([0 1],'polyval(c,x)',[-.451 .2836 .5641])ans = 1.4761

or by:

>>ssec2([0 1],'polyval([-.451 .2836 .5641],x)')ans = 1.4761

To find the second root, we give it a different interval:

>>ssec2([-1 0],'polyval(c,x)',[-.451 .2836 .5641]> ans = -0.8473

or:

>>ssec2([-1 0],'polyval([-.451 .2836 .5641],x)')ans = -0.8473

Since we are dealing with a polynomial, it is easier to find all the
roots at once with *roots*:

>>help rootsROOTS Find polynomial roots. ROOTS(C) computes the roots of the polynomial whose coefficients are the elements of the vector C. If C has N+1 components, the polynomial is C(1)*X^N + ... + C(N)*X + C(N+1). See also POLY. >>roots([-.451 .2836 .5641])ans = 1.4761 -0.8473

NOTE: Sometimes SSEC2 will not return a satisfactory root,
**but** that **does not **mean that there
are **no** roots. Often the interval chosen does not
bracket a root you want to see or too few trials are done. In either
case, choose a new interval or increase the number of trials. If a
root is still not found, plot the function to see what is wrong or
use another root finding function such as ROOTS, ROOT3, FZERO, or
POLY.

The solution of Example Problem 5.5 in Reklaitis with these functions
is trivial. We want a root of the cubic equation:

f(x) = x - 2(1-x) ^{3}= 0

If we look in the interval: 0 < x < 1, with
**ssec2** we find:

>>ssec2(0:1,'x-2*(1-x)^3')ans = 0.4102<----this is the root

We can easily adapt this procedure to solve many of the trial and
error problems in process calculations. For example, we can get back
to finishing *Example Problem 2.8*
from the text. It is convenient to define a function that will act as
our f(x) for this problem. **Er28** does this:

function er=er28(x) ti=sept(1,x,[.83333,0,.16667],5); sep(ti,5,[4 3]) er=.25-ns(3,1)./ns(3,2);

We need to go through the process of setting *ns* as we did up
to the point where we discovered that we needed a trial and error
solution. All we really need is for the feed stream (5) to be set and
the shape of *ns* to be correct. We will also need to have the
right *cnms* to use **showm**. Recall that the
ratio of benzene (component 1) to the other HC (component 2) in the
extract stream (stream 3) must be 0.25; hence the form of
**er**. If we execute **er28** with our
first assumed value for the mass fraction of benzene in stream 3:
0.1, we get:

>>er28(.1)ans = -0.9470

While an initial guess of 0 gives:

>>er28(0)ans = 0.2500

Obviously the answer must be between 0 and 0.1 so we can use
**ssec2** to find it:

>>ssec2([0 .1],'er28(x)')ans = 0.0231

We also have set our flow matrix to give the desired result:

>>showm(5 ,[4 3],10,2)Compound Inlet | Outlet Stream 5 | 4 3 Total benzene 700.00 | 625.00 75.00 700.00 other hc 300.00 | 0.00 300.00 300.00 SO2 3000.00 | 125.00 2875.00 3000.00 Total 4000.00 | 750.00 3250.00 4000.00

The problem asks for the percent recovery of benzene - the amount of
benzene in the raffinate (stream 4) to the amount of benzene in the
refinery stream (stream 1) is:

>>ns(4,1)/ns(1,1)ans = 0.8929

Therefore, 89.3 percent of the benzene is recovered.

As a final example in the chapter, we will consider extracting a
solute out of an "oil" mixture with "water". One formulation of this
type of problem is in terms of two feed streams containing:

Stream Water Oil Solute 1 W 0.0 S^{1 }2 0.0 O S^{2}

If we mix these streams and separate the result, we can produce
two products: a light stream (of mainly oil and solute) and a heavy
stream (of mainly water and solute.) In the ideal case the "oil"
stream contains no water and vice versa. We may also be able to find
data that determines the distribution of the solute between them. One
approximation to the distribution is given by:

with D a known distribution coefficient. We can actually solve for
all flow rates in such a separation by solving a quadratic equation
involving the amount of solute in one of the phases. We will
demonstrate another way to do this by using our modular functions and
**ssec2**. First set globals and initialize
*cnms* and *ns* with
**start301**,choice** (3) Other** :

You have chosen to not use any compound data If you wish to enter data for a compound, exit (control-C) then execute "create", then run a "new session", choice (1) Enter the number of compounds:3Compound names must be less than 25 characters long! Give the name of compound # 1:waterCompound names must be less than 25 characters long! Give the name of compound # 2:oilCompound names must be less than 25 characters long! Give the name of compound # 3:soluteEnter the number of streams:5

Set the feed streams:

>>ns(1:2,:)=[100 0 2;0 50 10];>>mix(1:2,3)>>showm(1:2,3,10,2)Compound Inlet Outlet Stream 1 2 Total 3 water 100.00 0.00 100.00 100.00 oil 0.00 50.00 50.00 50.00 solute 2.00 10.00 12.00 12.00 Total 102.00 60.00 162.00 162.00

Then write a program to try various fractions of the solute sent to
the oil phase:

function out=phase(dis,x) % Phase: calculates solute distribution in a two-phase separation. % function out=phase(dis,x) % Argument List: % Argument Gives % dis distribution coefficient % x fraction of solute % going into the top phase % phase tries various fractions of solutes sent to a phase % and returns the difference between the calculated and desired % distribution coefficients % See also PHASEG, PICPHASEG % OKB, TYLC global ns out=[]; for k=1:length(x) sep([0 1,x(k)],3,[4 5]);% [No water in 4, All oil in 4, x(k) of Solute in 4] temp=sum(ns(4:5,:)');% total flow in each stream out of separator: row xs=ns(4:5,3)./temp';% fraction solute in streams 4 and 5 subt=xs(1)/xs(2);% ratio of the fractions of solute in 4 to 5 if subt==NaN % check to see if the fractions were both zero subt=0; end out=[out,dis-subt];% we want subt to equal dis end

Note that the first argument of **phase** is the
distribution coefficient and the second argument is the fraction of
the solute going into the top phase. We will assume that in our case
D = 3, and for this value try a couple of values for the fraction to
send to the top:

>>phase(3,.9)<-- D = 3 and the frac. to top = .9ans = -11.9803 >>phase(3,.1)<-- D = 3 and the frac. = .1ans = 2.7595

OK, there must be a solution for D = 3 somewhere between 0.1 and
0.9.

>>ssec2([.1 .9],'phase(3,x)')ans = 0.6226

There it is. We can display our results and check them.

>>showm(3,4:5,11,4)Compound Inlet Outlet Stream 3 4 5 Total water 100.0000 0.0000 100.0000 100.0000 oil 50.0000 50.0000 0.0000 50.0000 solute 12.0000 7.4708 4.5292 12.0000 Total 162.0000 57.4708 104.5292 162.0000 >>x45=[7.4708 4.5292]./[57.4708 104.5292]x45 = 0.1300 0.0433 >>.13/.0433ans = 3.0023<-- Pretty close to the specification.

A more generalized version of **phase** is
**phaseg**. Help on it shows:

>>help phasegPhaseg: Calculates solute distribution in 2-phase separation using modules. function er=phaseg(iw,ioil,isol,dis,in,out,x) Argument Gives iw index of the water in cnms. ioil index of oil in cnms. isol index of the solute in cnms. dis distribution coefficient in index of the inlet stream out indices of the exit streams x fraction of solute going into the top phase er returns the difference between dis and the ratio actually set for the solute in the exit streams. For additional help and picture see PICPHASEG See also PHASE Example: >> ssec2([0.1 0.9],'phaseg(1,2,3,3.0,3,[4 5],x)'); revised help by TYLC

Here is the image of phaseg:

Here is the same simulation that we showed for *phase* with
D=3.0 and a guess of 0.1 for the fraction in the top stream.

>>phaseg(1,2,3,3.0,3,[4 5],0.1)<-- D = 3 and the frac. = .1ans = 2.7595

It may be used in *ssec2* just like *phase*

>>ssec2([.1 .9],'phaseg(1,2,3,3.0,3,4:5,x)')ans = 0.6226

[Go to previous section: 1.4| Go to next section: 1.6]

[ Go to CENG301 Homepage | CENG301 Notes Table of Contents]

Last modified July 31, 1998.