PROCESS MODELING SIMULATION &

OPTIMIZATION

Submitted By
Aalok Nishad(22BT01010)
Vishal Purohit(22BT01017)
CHEMICAL ENGINEERING

Submitted To
Mrs. Priyanka Pandya
GSFC University, Vadodara

Date of Submission: November 7th, 2025

 Question based on Model and Prototype
Simulitude

1. A geometrically similar model of an air duct is built to 1:25 scale and

tested with water which is 50 times more viscous and 800 times more
dense than air. When tested under dynamically similar conditions, the
pressure drop is 2 bar.
Find corresponding pressure drop in full-scale prototype.

Step 1: Given that

Model: water,

prototype: air

Step 2:
lm / lp = 1/25
μm = 50 μp
ρm = 800 ρp

Pressure drop in model = 2 bar

Step 3:
To find = pressure drop in prototype?
Step 4:
For dynamic viscosity similarity between prototype & model,
Reynolds number for both of them should be equal:
𝑉𝑝 𝜇𝑝 𝜌 𝑙𝑚
𝑉𝑚 = � 𝑚 × �
𝑚× 𝑝
�𝑝 �
1 � 1
�

𝑉𝑝 = 1 × × 800 ×
5016 25
=
25
𝐹 𝑚×𝑎 𝜌 ×𝑉×
𝑃 = = 𝑣
𝐴 𝐴
𝑇×𝐴
=
𝜌 × 𝐿3 × 𝑣 𝜌 ×1×
=
𝑇×𝐴 𝑣
= 𝑇×𝐴

=𝜌 𝑣²
Step 5:

Pressure drop in prototype:

𝜌 � 2 = 1 16
𝑝
(25)2 = 5.12 × 10
𝑉
𝜌 800
𝑝 �
2
𝑚 𝑚

Step 6:
Pressure drop in prototype:

= 5.12 × 10−4 × 2 = 10.24 × 10−4 bar

= 102.4 Pa

Step 7: Conclusion

Even though the model shows a pressure drop of 2 bar, the actual full-
scale air system will experience a very small pressure drop (≈102 Pa) due
to the low density of air.
2. A geometrically similar model of water is built at 1/25 scale and tested
with air, which is 50 times less viscous and 1/800 times denser than
water. When tested under dynamically similar conditions, the pressure
drop is 4 × 10⁻² bar in model.
Find corresponding ΔP in full-scale prototype.

Step 1: Given that

Model = air,

Prototype = water

Step 2:
lm / lp = 1/25
μm / μp = 1/50
ρm / ρp = 1/800

Pressure drop in model = 4 × 10⁻² bar = 10.24 × 10⁻⁴ bar

Step 3:
To find: Pressure drop in prototype?
Step 4:
For dynamic similarity between model & prototype:
Reynolds no. for both must be equal:
𝜌𝑉𝐿 𝜌𝑉𝐿
∣𝑝= ∣𝑚
𝜇 𝜇
𝑉𝑝 𝜇 𝜌 𝑙𝑚
×
𝑉𝑚 = × �𝑝
𝑝 𝑚

� 𝑚 �𝑝 1 �
� 1 �
= 50 × × = 0.0025
𝐹 800 𝑚 ×25
𝑃= = 𝑎
𝐴
𝐴
𝜌×𝑣×𝑣
=
𝑇×𝐴
𝜌 × 𝑙3 × 𝑣
= =
𝑇×𝐴

𝜌×𝑙 ×𝑣
𝑇×𝐴
=𝜌 𝑣²

Thus pressure ∝ ρV²

Step 5:
Pressure drop in prototype:
 𝑃𝑝 𝜌 �𝑝2
𝑃𝑚 𝑝 �𝑉 = 800 × (0.0025)2
= 2�
𝜌�
�

𝑃𝑝 = 0.005 𝑃𝑚
2
𝑃𝑝 = 0.02 bar = 0.02 × 105 = 2000 N/m
�

= 5.12 × 10−6 bar

Step 6: Conclusion
Pressure in full-scale prototype, where air is considered as fluid, must be
lower than value found with other case.
3. A geometrically similar model of H₂O is built to 1/40 scale and tested
with air, which is 50 times more viscous and 800 times more denser than
air. When tested under dynamically similar conditions, pressure drop in
model is 2 bar.
Find corresponding ΔP in full-scale prototype.
Step 1: Given

Scale: 𝑙𝑚/𝑙𝑝 = 1/40

𝜇𝑚/𝜇𝑝 = 50 ⇒ 𝜇𝑝/𝜇𝑚 = 1/50

Model fluid = water, Prototype fluid = air

𝜌𝑚/𝜌𝑝 = 800 ⇒ 𝜌𝑝/𝜌𝑚 = 1/800

Model pressure drop: Δ𝑃𝑚 = 2 bar

Step 2:

Dynamic similarity (Reynolds equality)

𝜌𝑉𝐿 𝜌𝑉𝐿
∣𝑝= ∣𝑚
𝜇 𝜇
𝑉𝑝 𝜇 𝜌 𝑙𝑚 1 1
⇒ = ⋅ ⋅
𝑝 𝑚

�𝑚 �𝑚 �𝑝 �𝑝 = ⋅ 800 ⋅
� � � � 50 40
= 0.4

Step 3:
For geometrically/dynamically similar internal flow,

Δ𝑃 ∝ 𝜌𝑉2 Δ 𝜌 �𝑝2 𝜌𝑝 𝑉𝑝 2 1
=( )
⇒ 𝑃𝑝 = 𝑝 �𝑉 )( 80 × (0.4) = 2 × 10 .
2 −4
𝜌 0
=
2
Δ𝑃 𝜌 𝑉
𝑚

𝑚
𝑚 𝑚 𝑚
Step 4:

Δ𝑃𝑝 = 2 × 10−4 × 2 bar = 4 × 10−4 bar = 40 Pa.

Step 5: Conclusion
The calculated pressure drop in the full-scale prototype is:

Δ𝑃𝑝 = 4 × 10−4
bar

This confirms that the prototype pressure drop is far smaller than the
model’s 2 bar, because air is much lighter and less viscous than water.
 Ko

PROCESS MODELLING SIMULATION AND OPTIMIZATION

CEC Assignment

SUBMITTED BY
Dhavalgiri Goswami (22BT01002)
Jaydevsinh Solanki(22BT01022)

SUBMITTED TO
Ms. Priyanka Pandya
 Illustration 4.1

Consider that benzoic acid (BA) is continuously extracted from

toluene using water as the solvent. The two streams are fed into
tank A (mixer) where they are stirred vigorously and the mixture is
then pumped into tank B (settler) where it is allowed to settle into
two layers (extract and raffinate). The upper toluene layer (raffinate)
and the lower water layer (extract) are removed separately. Model
the system to find out what fraction of BA has passed into the
solvent phase.

Solution

Step 1: The first step in modelling of any system is to schematically

represent the given system by understanding the actual process
occurring in the system.

(Fig. 4.1 Single-stage mixer settler)

Step 2: The second step is to see if there is a possibility of

simplifying the schematic representation of the system under
consideration.

The above system can be idealized and simplified as follows by

combining two stages into a single stage (Fig. 4.2), as we know that

Input – Output + Generation – Depletion = Accumulation

 (Fig. 4.2 Idealized single-stage solvent extraction system)

The various streams have been labelled and two material balances
have already been used: (a) conservation of toluene and (b)
conservation of water. These flow rates have been expressed on a
solute-free basis to simplify the analysis. The concentration of BA in
each stream has also been stated and this completes the
mathematical model.

Step 3: The third step is to make reasonable assumptions, and also

justify the assumptions made with reference to retaining accuracy.

So far, it has been assumed that:

(1) All flow rates are steady,

(2) Toluene and water are immiscible (basic criterion for

choosing a solvent), and

(3) The feed concentration c remains constant.

Step 4: Look for additional assumptions to be made in order to

balance the number of unknowns and the number of equations
required.

We assume that the mixer is so efficient that the two streams

leaving a stage are always in equilibrium with one another. This can
be expressed mathematically as:
 Y = mx

Where m is the distribution coefficient.

Mass balance for BA (component balance):

Input of BA = Rc (kg/s) Output

of BA = Rx + Sy (kg/s)

As the system is in steady state, Input = Output:

Rc = Rx + Sy (Eq. 4.2)

Substituting y = mx:

Rc = Rx + Smx

⇒ Rc = (R + mS)x

Therefore,

X = Rc / (R + mS) (Eq. 4.4) Y

= mRc / (R + mS) (Eq. 4.5)

Step 5: Express quantities in dimensionless form.

Fraction of BA extracted:

Sy / Rc = [1 / ((R / (mS)) + 1)] (Eq. 4.7)

Define:

Α = R / (mS)

E = Sy / Rc

Substituting:

E = 1 / (α + 1)

Multiplying numerator and denominator by (α – 1):

E = (α – 1) / (α² – 1)

Step 6: The
 final result should be expressed in generalized form for extension
to multistage systems.

Example 4.1

Compute the fraction of solute that could be extracted in a single-

stage solvent extraction using the numerical values of:

S = 12R, m = 1/8, and c = 0.1 kg/m³

Solution:

If S = 12R, m = 1/8, and c = 0.1,

X = Rc / (R + mS)

= 0.1R / [R (1 + 12/8)]

= 0.04
Therefore,

Y = mx = (1/8) × 0.04 = 0.005

Fraction of BA extracted:

E = (Sy) / (Rc)

= (12R × 0.005) / (0.1R)

= 0.6 = 60%

Alternatively,

E = 1 / (α + 1)

Where α = R / (mS)

= R / [(1/8) × 12R]

= 8/12

So,
 E = 1 / [(8/12) + 1]

= 12 / 20

= 0.6 = 60%

Conclusion:
The steady-state single-stage solvent extraction model
demonstrates how benzoic acid (BA) can be extracted from toluene
using water as a solvent under equilibrium conditions. By applying
mass balance and equilibrium relations, a mathematical model is
developed to determine the fraction of solute extracted into the
solvent phase.

The analysis shows that the extraction efficiency depends primarily

on the dimensionless ratio (α = R / mS) — the ratio of feed flow
rate to solvent flow rate at equilibrium. The derived expression:

E = 1 / (α + 1)

Reveals that increasing the solvent flow rate (thus decreasing α)

improves the extraction efficiency.

From Example 4.1, using given parameters (S = 12R, m = 1/8, c =

0.1 kg/m³), the computed extraction efficiency is 60%, which
validates the model and confirms that a significant portion of solute
can be extracted in a single-stage operation.

Overall, this case illustrates the importance of equilibrium relations

and flow ratios in designing and optimizing solvent extraction
systems.
PMSO CEC: Steady-state Two-stage Solvent Extraction

Submitted by:

Tirth patel (23BT01D032)

Avinash Prajapati
(23BT01D033)

Submitted to:

Ms. Priyanka Pandya

Subject:

Process Modelling, Simulation and Optimization

 Topic: Steady-state Two-stage Solvent Extraction

Illustration 4.2
Develop a model and obtain an expression for the extraction of BA in a steadystate operation with two
stages. The input data are same as in the previous problem. Each stage still consists of two tanks, a
mixer and a settler, with counter-current flow through the stages. The idealized flow system is shown
in Fig. 4.3, where the symbols have the same meaning as in the previous example, and the different
concentrations in a particular phase are distinguished by suffixes.

Solution

Step 1: System Description and Assumptions:-

● The system consists of two tanks (stages), each with its own mixer and settler.
● Streams: The feed (with concentration c) enters stage 1 with flow rate R, and solvent
stream enters with flow rate S.

Assumptions:-
● The feed concentration c remains constant.
● The mixer is so efficient that the two. streams leaving a stage are always in equilibrium with
each other.

Step 2: Applying the Equilibrium Relationship:-

● For each stage : The concentration in the solvent phase leaving that stage (y) is related to the
concentration in the other phase (x) as:
y1 = mx1 (for stage 1)..................(1)
y2 = mx2 ( for stage 2 )...............(2)
where m is the distribution coefficient. It is constant for the given system. A BA mass balance is now
taken for each stage:

Stage 1 Stage 2

Input of BA (kg/s): Rc + Sy2 RX1

Output of BA (kg/s): RX1 + SY1 Rx2 + SY2

Step 3: Setting Up Mass Balances :-

Since the operation is carried out in steady state, the BA must enter and leave a stage at the same rate,
i.e.

Stage 1: Re + Sy1 = Rx1 + Sy1….............(3)

Stage 2: Rx, = Rx, + Sy1….........................(4)

Step 4 : Express y1 in terms of other variables :

To get an expression for y1, in terms of known quantities, substitute x1, and y2 from Eqs (1) and (2)
in Eqs (3) and (4). Equation (3) becomes

Rc + mSx2 = (Ry1 / m) + Sy1…....................(5)

= [ (R + mS) / m ] y1…...................................(6)
And Equation No: 4 becomes
(Ry1 / m) = Rx2 + mSx2….........................(7)
Equation No:8 gives an expression for y1
y1 = ( m / R ) ( Rx2 + mSx2 )…........................(8)

Step 5 : Express x1 in terms of other variables

Since x2, is unknown, it is to be expressed in terms of known quantities. Substituting y1,
from Eq. (8) in Eg. (6), we get

Rc + mSx2 = [ (R + mS) / m ] (m / R) (Rx2 + mSx2)

 R^2c + mRSx2 = R^2x2 + mRSx2 + mRSx2 + m^2S^2x2

 R^2c = x1(R^2 + mRS + m^2S^2)

Equation No: 9 gives an expression for x1.

 x2 = [ R^2c / R^2 + mRS + m^2S^2 ]...................(9)
Step 6 : Substitute x2 into the expression for y1
For equation no : (9) to (8) we get :-

y2 = [ MRS(R + mS) / R^2 + mRS + m^2S^2 ]…............................................(10)

Step 7; Find the proportion of BA extracted

The proportion of BA extracted is given by the ratio Sy1/Rc. Sustitute the express of y1 from eq (10) and
(9).
Again, the proportion of BA extracted is

(Sy1 / Rc) = [ mS(R + mS) / R^2 + mRS + m^2S^2 ]….................................(11)

Multiplying and dividing the right-hand side (rhs) of Eq. (4.18) by m^2S^2, we get

(Sy1 / Rc) = [ (R/mS) + 1 / (R^2/m^2S^2) + (R/mS) + 1 ]..........................(12)

Step 8: Now introducing the dimensionless groups E = Sy1/Rc and a = R/mS, we get E

= (a + 1) / (a^2 + a + 1)..................................................................................(13)

Multiplying and dividing the rhs of the above equation by (a - 1), we get

E = (a^2 – 1) / (a^3 – 1).........................................................................................(14)

Step 9:Generalize for N stages

By the method of induction, the expression for E can be generalized for N stages. E

= (a^N – 1) / (a^N+1 – 1).................................................................................(15)

Example Number 4.2 : Compute the fraction of solute that can be extracted in a two-stage
counter-current solvent extraction using the same numerical values of S= 12R,
= 1/8, and c = 0.1 kg/m3, which were used for single-stage extraction, and
compare the results.

Diagram :-

Solution:

Given Data :-

● Number of stages : 2
● Solvent ratio : S = 12R
● Equilibrium constant : m = 1/8
● Concentration

Step No : 1 Calculate the fraction of solute in the raffinate phase from second stage ( X2 )
The fraction of solute remaining in the raffinate after two stages, x2 is calculated using the given
values for s, m and c .
Formula:-
x2 = [ R^2c / R^2 + mRS + m^2S^2 ]

Substitute the given values S = 12R and m = 1/8 c = 0.1 kg/m^3.

x2 = [ R^2(0.1) ] / [R^2 + (1/8)R(12R) + (1/8)^2 (12R)^2

] x2 = [ R^2(0.1) ] / [R^2 + (12/8)R^2 + (144/64)R^2 ]

Simplify the expression by factoring out

R^2: x2 = [R^2(0.1)] /
R^2[1+(12/8)+(144/64)] x2 = [0.1 /
1+1.5+2.25]
x2 = 0.1 / 4.75
x2 = 0.02105
Step No : 2 Calculate the fraction of solute in the extract from the first stage ( y1 )
The fraction of solute in the extract from the first stage, y1, is calculated using the given
values.
Formula :-
y1 = mRc(R+mS) / [R^2 + mRS + m^2S^2]
Substitute the given values S =12R and m =1/8,C =
0.1
Simplify the numerator and use the simplified denominator from Step 1:

y1 = (1/8)R(0.1)(R+1/8(12R)) / R^2+(1/8)^2 R(12R)^2

y1 = (1/8)R(0.1)(R+1.5R) / R^2(1+(12/8)+(144/64)
y1 = (1/8)R(0.1)(2.5R) / R^2(4.75)
y1 = (1/8)(0.1)(2.5) / 4.75
y1 = 0.00658

Step No : 3 Calculate the extraction factor ( a

) Formula:-
a = R / mS
a = R / (1/8)(12R)
a = R / (12/8)R
a = 1 / 1.5
a = 0.66667

Step No : 4 Calculate the fraction extracted ( E

). Formula:-
E = [a^2-1] / [a^3-1]

Substitute the value of a:

E = (0.66667^2)-1 / (0.66667^3)-1
E = 0.44444-1 / 0.29629-1
E = -0.55555 / -0.70369
E = 0.789

Conclusion:-

The fraction of solute extracted is 78.9%. This is a greater degree of extraction than the
60%achieved with a single stage.
 PROCESS MODELLING SIMULATION AND OPTIMIZATION
CEC ASSIGNMENT

SUBMITTED BY
DHAVAL JOSHI (21BT01019)

VISHAL PARMAR (

22BT01D077)

SUBMITTED TO
MS. PRIYANKA PANDYA

pg.
1
4.3 Steady-state Two-stage Cross-current
Illustration 4.3

Consider the idealized flow system shown in Fig. 4.4 where the symbols have the

Solution
The following assumptions are made:

1. Toluene and water are immiscible.

2. We are considering a pure solvent.

3. The mixer is so efficient that the leaving streams are always in equilibrium.

4. The concentration of BA in the feed is constant.

5. The solvent flow rate is the same in both stages.

6. All the quantities are expressed on a solute-free basis.

7. The composition of a stream leaving a stage is the same as that in the

Step-1

BA balance (mass component balance for each stage):

Stage 1 Stage 2

Input of BA (kg/h): Rc + 0 RX,+O

Output of BA (kg/h): SY, + RX, SY, + 4

In steady state:

Stage 1: Rc = Sy + Rx

Stage 2: Rx, = Sy, + Rx,

pg.
2
Step-2

Since there are four unknowns (x1, y1, x2, y2) and only two equations, therefore two

more equations are required and these are the equilibrium relations given by

Y1= m, ( Eq 4.24)

Y2 = m2 ( Eq 4.25)

The objective in this problem is to find the total fraction of BA that can be extracted

by the cross-current mode of extraction operation under steady-state conditions.

The fraction of BA extracted,

𝑠𝑦2 𝑠(𝑦1 =
𝑆𝑦1 = =𝐸
𝑅𝑐 𝑦2)
+
𝑅𝑐

Step-3
Hence, we have to obtain the expressions for y1 and y2 in terms of known quantities

by eliminating x1 and x2. So, from Eq. (4.24), we have

𝑦1
𝑥1 =
𝑚
step-4
From Eq (4.22) and (4.26),

𝑦
𝑅𝑐 = 1 �+ 𝑆) 𝑦1
�

𝑅 + 𝑆𝑦1 =
𝑚( �
�

Therefore,

𝑅
𝑥1 𝑐 + 𝑚𝑆
=
�
�

From Eq. (4.25), we have

𝑦2
𝑋2
= 𝑚

pg.
3
Substituting x1 and x2 from Eq (4.26) and (4.29) in Eq. (4.23), we get
𝑅
𝑐 𝑅𝑦 𝑚
𝑅 + 𝑚𝑠 2 + 𝑠𝑦2 = 𝑠 ) 𝑦2
� = (𝑅 +
� 𝑚 𝑚

𝑌2 𝑚𝑅2𝑐
=
Now
𝑅+𝑚𝑠
𝑀𝑟2𝑐 𝑚𝑅𝑐 𝑅
𝑚𝑅 𝑠 +𝑚𝑠][1+ +𝑚𝑠]
𝐸 = 𝐸[
𝑐 + 𝑚𝑠 (𝑅+𝑚𝑠) [ �
=
Now 𝑅
𝑅𝑐
+
2
𝑅 𝑅𝑐
�

𝑚𝑠
= 𝑚𝑠) (1
𝑅 + 𝑚𝑠)
+
=(

𝑅
𝑅

On multiplying and dividing both the terms in the parentheses by mS, we get R1
𝑅
1 𝑚𝑠 ]
𝐸 = [1 𝑅
𝑅
𝑚𝑠
++ 𝑚𝑠 + 1
1
We know that
𝑅
𝑚 =𝑎
𝑆
𝑎
1
𝐸 = + 𝑙 [1 𝑎+𝑙]
+
𝑎

1
𝐸 = 2𝑎 +
(𝑎 + 1)2

pg.
4
Example 4.4

Compute the fraction of solute that can be extracted in a two-stage cross-current solvent extraction
using the same numerical values of S = 6R, m = 1/8, and c = 0.1 kg/m3, which were used for single-stage
extraction, and compare the results.

Solution

If S = 6R (since 2s here = S in counter-current), m = 1/8, c = 0.1 kg/m3,

Now

𝑅 8
𝑅 = 6 = 1.3333
𝛼 =𝑚𝑠= 1
( () )
6𝑅
8

1
𝐸 = 2(1.3333) + = 0.6735 = 67.35%
(1.3333 + 1)2

Comparison For counter-current steady-state two-stage extraction,

Here, a = 0.66667 (since S = 12R)

Now
𝑎−1
𝐸
= 𝑎+1

By putting value E= 0.789

Final answer is 𝐸 =

79%

CONCLUSION

The conclusion of this example problem is a comparison of extraction efficiency between a two-stage
cross-current process and a two-stage counter-current process for the same overall solvent usage.

pg.
5
PROCESS MODELLING SIMULATION AND OPTIMIZATION

CEC Assignment

SUBMITTED BY

Nishit Tiwari(23BT01D028)

SUBMITTED TO

Ms. Priyanka Pandya

EXAMPLE 4.3 ANALYSIS FOR CONVEXITY AND CONCAVITY

For each of these functions

(a) f(x) = 3x2

(b) f(x) = 2x

(c) f(x) = -5x2

(d) f(x) = 2x2 - x3

determine if f(x) is convex, concave, strictly convex, strictly concave, all,

or none of these classes in the range - ∞ ≤ x ≤ ∞.

For a multivariate function, the nature of convexity can best be

evaluated by examining the eigenvalues of f(x) as shown in Table 4.1
We have omitted the indefinite case for H, that is when f(x) is neither
convex or concave.

Table (4.1) Relationship between the character of f(x) and the state of H(x)

f(x) H(x) All the eigenvalues

of H(x) are
Strictly convex Positive-definite >0

Convex Positive-semidefinite ≥0

Concave Negative-semidefinite ≤0

Strictly concave Negative-definite <0

Solution :

(a) f "(x) = 6, always positive, hence f(x) is both strictly convex and
convex.

(b) f "(x) = 0 for all values of x, hence f(x) is convex and concave. Note
straight lines are both convex and concave simultaneously.

(c) f "(x) = - 10, always negative, hence f(x) is both strictly concave and
concave.
 (d) f "(x) = 6 - 3x; may be positive or negative depending on the value
of x, hence f(x) is not convex or concave over the entire range of x.

EXAMPLE 4.5 DETERMINATION OF POSITIVE-DEFINITENESS OF A

FUNCTION

Repeat the analysis of Example 4.4 for f(x) = x12 + x1x2 + 2x2 + 4

Solution :

𝟐 𝟏
H(x) = (𝟏 𝟎)

The eigenvalues are 1 + √2 and 1 - √2 , or one positive or one negative

value. Consequently, f(x) does not fall into any of the categories in Table
4.1. We conclude that no unique extremum exists.

EXAMPLE 4.9 CONSTRUCTION OF A CONVEX REGION

Construct the region given by the following inequality constraints; is it

convex?

X1≤ 6 ; x2 ≤ 6; x1≥ 0 ; x1 + x2≤ 6; x2 ≥ 0

Solution :

See Figure E 4.9 for the region delineated by the inequality constraints. By
visual inspection, the region is convex. This set of linear inequality
constraints forms a convex region because all the constraints are
concave. In this case the convex region is closed.

FIGURE E4.9

Diagram of region defined by linear inequality constraints.

Conclusion :

The analysis demonstrates that the curvature properties of functions—

convexity and concavity—are governed by the nature of their second
derivatives or Hessian eigenvalues. A positive-definite Hessian indicates
strict convexity, guaranteeing a unique minimum, while a negative-
definite Hessian implies strict concavity and a unique maximum. In
contrast, functions with mixed-sign eigenvalues lack definite curvature
and no unique extremum exists. Moreover, the set of linear inequality
constraints defines a convex region, forming the foundation for various
optimization and design applications in engineering and economics.
PROCESS MODELLING SIMULATION AND OPTIMIZATION

CEC ASSIGNMENT

SUBMITTED BY

Vishwa Gandhi (21BT01014)

Pranav Parmar (23BT01D031)

SUBMITTED TO
Ms. Priyanka Pandya
 INTRODUCTION

The numerical problem associated with this topic is based on mass

balance and equilibrium relationships in a counter-current absorption
column.

It assumes that:

•The gas enters at the bottom and the solvent enters at the top of the
column.
•The solute transfer between phases follows the equilibrium

relation y^* = Kx.

•Flow rates of gas (G) and liquid (L) remain constant on a solute-
free basis.

Illustration 4.8

𝐺𝑁+1 kg mol/s of a wet gas containing 𝑌𝑁+1 kg mol of solute/kg mol of

solute is to be stripped from the gas by absorption in 𝐿𝑜 kg mol/s of lean

wet gas is fed into the base of a plate absorption column where the

𝑋𝑜 kg mol/kg mol of lean oil and the solute in the exit gas is 𝑌1 kg
oil, which is fed at the top of the column. The solute in the entering oil is

mol/kg mol of wet gas. The equilibrium constant 𝐾𝑚 is given as 𝐾𝑚 =

𝑌𝑀/𝑋𝑚,where 𝑌𝑀
and 𝑋𝑚 are the mole fractions of the solute in the gas and liquid phases,

expressed in terms of the absorption factor A = ( 𝐿𝑜/K𝐺𝑁+1) and the

respectively. Show that the performance of the absorber can be

number of ideal stages (N) by the Kremser-Brown equation (which relates

the fractional solute recovery to the absorption factor):
𝑌𝑁+1 − 𝑌1/𝑌𝑁+1 − 𝑌𝑜 = 𝐴𝑁 + 1 − 𝐴/𝐴𝑁 + 1
 𝐺𝑁+1=𝐺𝑚(1-𝑌𝑚)
Y=Y/(1+Y) ---Y=Y/(1-Y)
Considering Fig. 4.10 and taking material balance over plate m gives

𝐿𝑂(𝑋𝑚- 𝑋𝑚−1) = 𝐺𝑁 + 1 (𝑋𝑚+1-𝑋𝑚)

(The gas and liquid ratios on a solute-free basis are same on every plate,

eliminate 𝑋𝑚and 𝑋𝑚−1 from Eq. (4.86). But

though the actual L and G may be different on each plate.) We have to

𝑦𝑚𝐺𝑚= 𝑋𝑚𝐺𝑁+1

𝑋𝑚𝐿𝑚=𝑋𝑚- 𝑋𝑚𝐿𝑜
We know that the equilibrium relation is

𝑌𝑚 = 𝐾𝑚𝑋𝑚
𝑋𝑚
𝑦𝑚 = 𝐿0
𝐺𝑁+1 𝐺𝑚 𝑚
/ =𝐾
𝐿𝑚

𝐺𝑁+1
𝑋𝑚 = ( ) )𝑦𝑚
𝐾𝑚
/
𝐿0 �
(𝐿
� 𝐺
�
�
Similarly,
𝐺𝑁+1 𝐿𝑚−1
𝑋 = ( − 1)
𝐺𝑚−1
)(
𝐿0 𝐾𝑚 𝑚

Substituting the values of 𝑋𝑚 𝑎𝑛𝑑 𝑋𝑚−1 and from Eq. (4.88) and (4.89)
respectively, in Eq. (4.86) gives

= 𝐺𝑁+1(𝑌𝑚+1 − 𝑌𝑚)(4.90)

We know

𝐿
that
= 𝐴𝑚(absorption factor for plate 𝑚)
𝑚
𝐾𝑚 𝐺𝑚
𝐿𝑚−1
= 𝐴𝑚−1(absorption factor for plate 𝑚 − 1)
𝐾𝑚−1𝐺𝑚
−1

𝑌𝑚[𝐺𝑁+1𝐴𝑚 + 𝐺𝑁+1] = 𝑌𝑚+1𝐺𝑁+1 + 𝑌𝑚−1(𝐺𝑁+1𝐴𝑚−1)(4.92)

𝑌𝑚+1 + 𝐴𝑚−1𝑌𝑚−1
⇒ 𝑌𝑚 = (4.93)
1+
𝐴𝑚

bottom of the column as A, and rearranging Eq. (4.93) in the form of a

Taking the mean value of the absorption factor between the top and the

finite difference equation,

𝑌𝑚+1 − (1 + 𝐴)𝑌𝑚 + 𝐴𝑌𝑚−1 = 0(4.94)
 Putting Eq. (4.94) into the operational form and factorizing,

[E² − E(1 + A) + A]Yₘ = [(E − 1)(E − A)]Yₘ = 0 (4.95)

which has the solution

Yₘ = C₁ + C₂Aᵐ⁻¹ (4.96)

Put

m = 0, Y₀ = C₁ + C₂ (4.97)

and
m = 1, Y₁ = C₁ + C₂A (4.98)

Similarly,

m = N + 1, Yₘ𝗁₁ = C₁ + C₂Aᴺ⁺¹ (4.99)

∴ Yₘ𝗁₁ − Y₀ = C₂ (Aᴺ⁺¹ − 1) (4.100)

and

Yₘ𝗁₁ − Y₁ = C₂ (Aᴺ⁺¹ − A) (4.101)

Dividing Eq. (4.101) by Eq. (4.100), we get

(Yₘ𝗁₁ − Y₁) / (Yₘ𝗁₁ − Y₀) = (Aᴺ⁺¹ − A) / (Aᴺ⁺¹ − 1) (4.102)

Equation (4.102) is the well − known Kremser − Brown equation.

The physical significance of this equation is: The normalized

fractional recovery of the solute in

a multicomponent absorber operated in a counter − current mode of

operation depends on the absorption factor.

 CONCLUSION

From the above calculations, it can be concluded that the absorption

column is capable of achieving the required solute removal within the
estimated number of stages. The results show a clear relation between
the absorption factor and the overall recovery of the solute. As the
absorption factor increases, the number of stages required for a given
separation decreases, which agrees with the theoretical [Link]
numerical also helped in understanding how the Kremser–Brown equation
can be applied to predict column performance without doing complex
experimental trials.
Process Modelling Simulation
and Optimization

CEC Assignment

Submitted By

Ayush Chhodvadiya (23BT01D025)

Yuvraj Rajput (23BT01D035)

Submitted To

Ms. Priyanka Pandya

GSFC University, Vadodara

Illustration 4.9

A continuous flow distribution column is fed with a binary mixture of A and

B. The relative volatility a of the mixture is constant. By considering the
plate above the feed plate, derive the difference equation that relates the
liquid composition to the plate number n if the overall plate efficiency is
100%. Reduce the equation to a linear one with constant coefficients.
Show that the axis translation required involves the points where the
equilibrium and operating lines intersect.

Solution: -

𝑥𝑛
Mole fraction of the more volatile component in the liquid leaving

𝑛
plate

𝑥𝐷
Mole fraction of the more volatile component in the distillate (product
from condenser)

𝑦𝑛∗
with the liquid of composition 𝑥𝑛
Mole fraction of the more volatile component in vapor in equilibrium

𝑦𝑛 Mole fraction of the more volatile component in vapor leaving plate 𝑛

𝐿 Constant molar flow rate of liquid down the column
𝑉 Constant molar flow rate of vapor up the column
𝐷 Molar flow rate of distillate product from the condenser

Fig.4.1 Multistage Distillation

Considering Fig.4.1 and writing the material balance for the more volatile
component around the section which includes the nth plate and the
condenser,

𝑉𝑌𝑛−1 = 𝑂𝑋𝑛 − 𝐷𝑋𝐷 = 0 (1)

This is known as the operating line for the enriching section. The definition
of relative volatility is
𝑦∗ (1−𝑥 )

𝛼 (1−𝑦
= 𝑛∗ )𝑥𝑛 (2)
� 𝑛
�
This is the equation of the equilibrium line. Combining Eqs (1) and (2), we
get
𝑥 𝑥
− =0
𝐷𝑥𝐷

+
[
𝐷𝑥𝑥
+
𝐷𝑛(𝛼−1)−𝛼𝑉
]𝑥
(3)

𝑛 𝛼−1 𝑂(𝛼−1 𝑛− 𝑂(𝛼−1)

𝑛−1 1
)
Let,
1
𝑎= 𝐷𝑥𝐷(𝛼 − 1) − 𝛼𝑉 𝐷𝑥𝐷
𝛼 − 1, 𝑏 = 𝑂(𝛼 − 𝑂(𝛼
, 𝑐 =− 1)
1)
Equation (3) becomes

𝑥𝑛𝑥𝑛−1 + 𝑎𝑥𝑛 + 𝑏𝑥𝑛−1 + 𝑐 = 0 (4)

This is known as Riccati equation. Put 𝑥𝑛 = 𝑥𝑎 + 𝛿

(𝑥𝑛 + 𝛿)(𝑥𝑛−1 + 𝛿) + 𝑎(𝑥𝑛 + 𝛿) + 𝑏(𝑥𝑛−1 + 𝛿) + 𝑐 = 0 (5)

⇒ 𝑥𝑛𝑥𝑛−1 + 𝑥𝑛(𝑎 + 𝛿) + 𝑥𝑛−1(𝑏 + 𝛿) + 𝛿2 + (𝑎 + 𝑏)𝛿 + 𝑐 = 0

If the value of 𝛿 is such that

𝛿2 + (𝑎 + 𝑏)𝛿 + 𝑐 = 0 (6)
we have

𝑥𝑛𝑥𝑛−1 + 𝑥𝑛(𝑎 + 𝛿) + 𝑥𝑛−1(𝑏 + 𝛿) = 0 (7)

⇒1+ + =0
𝑎+𝛿 𝑏+𝛿
(8)
𝑥𝑛− 𝑥𝑛

Put 𝑉𝑛 =
1

1/𝑥𝑛

∴ 𝑉𝑛(𝑏 + 𝛿) + 𝑉𝑛−1(𝑎 + 𝛿) + 1 = 0 (9)

𝑉𝑛(𝑏 + 𝛿) + 𝑉𝑛−1(𝑎 + 𝛿) = −1 (10)

Putting this in the operational form,

[𝐸(𝑏 + 𝛿) + (𝑎 + 𝛿)]𝑉𝑛−1 = −1 (11)

 ∴ 𝐸 = − (𝑏+𝛿)
(𝑎+𝛿)
(12)

The complementary solution is

𝑉𝑛−1 = 𝐾[− ]
(𝑎+𝛿) 𝑛−1

(𝑏+𝛿)
(13)

The particular solution is

𝑉𝑛
= (𝑏+𝛿)+
−1 (14)
−1 (𝑎+𝛿)

𝑉
Hence the complete solution is
= 𝐾[− ] −
(𝑎+𝛿) 𝑛−1
1 (15)
𝑛−1 (𝑏+𝛿) (𝑏+𝛿)+
(𝑎+𝛿)

∴ 𝑉𝑛 (𝑎+𝛿)
= 𝐾[− ]
− (𝑏+𝛿)+(𝑎+𝛿)
(16)
𝑛 1
(𝑏+𝛿)

Since 𝑉𝑛 = 1/(𝑥𝑛 − 1),𝐾 =constant,

𝑥𝑛
=𝛿+
1
(𝑎+𝛿) 𝑛 1
𝐾[− ] −
(17)
(𝑏+𝛿) (𝑏+𝛿)+(𝑎+𝛿)

(𝑥𝑖, 𝑦𝑖):
Let the intersection of the equilibrium and operating lines occur at the point

𝑉𝑦𝑖 − 𝑂𝑥𝑖 − 𝐷𝑥𝐷 = 0 (18)

𝛼 =(1−𝑦
𝑖 𝑦 (1−𝑥 )
𝑖
𝑖)𝑥𝑖 (19)
On eliminating 𝑌𝑖, we have from the above two
equations,

𝑥�2 + (𝑎 + 𝑏)𝑥𝑖 + 𝑐 = 0 (20)

�

With the usual notation for 𝑎, 𝑏,and 𝑐,

𝛿 = 𝑥𝑖
Therefore, the translation of axis required to linearize the Riccati equation
corresponds to the new origin of coordinates at the intersection of the
equilibrium and operating lines.
Conclusion

This derivation demonstrates the mathematical transformation of the

nonlinear vapour-liquid equilibrium relation into a solvable linear form
using relative volatility and coordinate translation. The resulting
expression provides insight into plate-to-plate variation of composition,
enabling precise modelling of the enriching section in distillation design
and simulation.
 PROCESS, MODELLING,
SIMULATION & OPTIMIZATION

CEC ASSIGNMENT

SUBMITTED BY
HARSH CHAUHAN (23BT01D024)
YAGNESH SONI (23BT01D038)

SUBMITTED TO
MS. PRIYANKA PANDYA
MODELS IN HEAT TRANSFER OPERATION

Illustration 5.4

1.26 kg/s of sulphuric acid of heat capacity 1508 J/kg K is to be cooled in a

two-stage counter-current cooler of the following type. Hot acid at 174°C
is fed to a tank where it is stirred well in contact with cooling coils. The
continuous discharge from this tank at 88°C flows to a second stirred tank
and leaves at 45°C. Cooling water at 20°C flows into the coil of the second
tank and then to the coil of the first tank. The water is at 80°C as it leaves
the coil of the hot acid tank. To what temperature would the contents of
each tank rise if, due to trouble in the supply, the cooling water is
suddenly stopped for 1 h? Calculate also the intermediate water
temperature before failure. The capacity of each tank is 4536 kg of acid
and the flow rate of water is 0.975 kg/h.
Solution:
GIVEN THAT:

MA = Mass flow rate of H2SO4 = 1.26 kg/s

CPA = Heat capacity of H2SO4 = 1508 J/kg K

CPW = Heat capacity of Water = 4187 J/kg K

T0A = 174°C

T1A = 88°C

T2A = 45°C

T0W = 20°C

T1W =?

T2W = 80°C

W = Capacity of each tank = 4536 kg acid

MW Flow rate = 0.975 kg/h

SOLUTION

For Tank 1

Considering heat balance around the tank 1:

Rate of heat input by cooling water = Rate of heat output by sulfuric acid

∴ MWCPW (80 – T1W) = MACPA (174 – 88)

∴ 0.975 * 4107 (80 - T1W) = 1.26 * 1508 (86)

∴ 80 - T1W = 40.8

∴ T1W = 39.2 ≈ 40°C

For tank 2

MWCPW ΔTW = MACPA ΔTA

0.975 * 4107 (T1W – 20) = 1.26 * 1504 (88 – 45)

T1W = 40.34°C ≈ 40°C

Before failure: Steady State

After failure: Unsteady state
  Given condition is when cooling water supply is stopped for 1 hour =
3600S
 The rate of cooling water = 0
 Unsteady state condition prevails as the acid is still flowing at the
same rate as before

Now, for this condition, let assume

Outlet acid temperature of tank 1 T1A =?

Outlet acid temperature of tank 2 T1A =?

Where θ = 0 Sec

θ = 3600 Sec

Heat balance for unsteady state condition:

For tank 1
𝑑 T1A
MCPA T0A - MCPA T1A = WCPA
𝑑θ -----------(1)
For tank 2
𝑑 T2A
MCPA T1A - MCPA T2A = WCPA
𝑑θ -----------(2)

From equation 1
𝑑 T1A
MCPA (T0A - T1A) = WCPA
𝑑θ

* 𝑑θ =(𝑇0𝐴−𝑇1𝐴)
𝑀 𝑑 T1A
𝑊

With respect to T1A ∴ -ln (T0A - T1A) If

with respect to T0A ∴ ln (T0A - T1A)

Integrating above equation

𝑀
* θ + C = -ln (T0A - T1A)
𝑊

Now rearrange
* θ + C-------------------------------------------------------------------(3)
𝑀
-ln (T0A - T1A) =
𝑊
Applying boundary conditions

When θ = 0 sec, T1A = 88°C

∴ -ln (T0A - T1A) = *θ +C

𝑀
𝑊

C = -ln (T0A - T1A) or C = -ln (T0A - 88)

Put value of C in 3 equation

* θ -ln (T0A – 88)
𝑀
-ln (T0A - T1A) -
𝑊

*θ
𝑀
or -ln (T0A - T1A) + ln (T0A – 88) =
𝑊

multiply by (-ve)
*θ
𝑀
ln (T0A - T1A) -ln (T0A - 88) = -
𝑊

∴ ln a – ln b = ln 𝑎
𝑏

∴ ln [T0A − T1A 𝑀
*θ
T0A − 88 𝑊
]=-
𝑀
−( ∗
T0A − T1A θ)
𝑒
Or = 𝑊

T0A − 88 −( 𝑀∗ θ)
Or (T0A - T1A) = (T0A - 88) * 𝑒 𝑊

Or T1A = T0A – (T0A – 88) * 𝑒 −(

𝑀
∗
θ)
𝑊
-----------(4)

Put value of T0A

T1A = 174 – (174 – 88) * 𝑒

1.26
−( ∗ θ)
4536

T1A = 174 – 86 𝑒(−2.78∗10^−4 ∗ θ)--------------------------------------------------------------------- (5)

Now θ = 3600 Sec, putting in equation 5

T1A = 174 – 86 𝑒(−2.78∗10^−4 ∗ 3600)

T1A = 142.4°C -----------(6)

From equation 2
𝑑 T2A
MCPA T1A - MCPA T2A = WCPA
𝑑θ
𝑑 T2A
MCPA (T1A - T2A) = WCPA
𝑑θ -----------(7)
𝑑 T2A
𝑀
(T1A - T2A) =
𝑊 𝑑θ

𝑑
T2A
𝑀
T 2A =
𝑀
T1A--------------------------------------------------------------------------(8)
𝑑θ
𝑊
+
𝑊

Putting value of T1A from equation 5 in equation 8

𝑀 [ 174 – 86 𝑒(−2.78∗10^−4) θ ]--------------------------------------------------------------(9)

𝑊
=

Equation 9 is the first order differential equation that can be solved using
method of integral factor

Basic IF order of differential equation

𝑑𝑦
𝑑𝑥
+ p (x) y = Q (x)-----------------------------------------------------------------------------(10)

By the method of IF

y * IF = ∫ 𝑄 (𝑥) IF dx + C--------------------------------------------------------------------(11)

where, IF = 𝑒 ∫ 𝑝 (𝑥) 𝑑𝑥--------------------------------------------------------------------(12)

standard form:
𝑑
𝑑θ T2A + P (θ) T2A = Q (θ)------------------------------------------------------------------(13)

T2A * IF ∫ 𝑄 (θ) IF dθ + C-------------------------------------------------------------------(14)

IF = 𝑒 ∫ 𝑝 (θ) 𝑑θ--------------------------------------------------------------------------------(15)

Since
= 𝑃(θ)}-------------(16)
𝑀
𝑀
* T2A = P (θ) T2A {
𝑊 𝑊

𝑀
𝑊
= P (θ)------------------------------------------------------------------------------------------(17)
IF = 𝑒 ∫ 𝑝 (θ) 𝑑θ -------- (17 i)

=𝑒 ∫ ∗ 𝑑θ from 17 (i)
𝑀

IF = 𝑒 ∫ ∗ 𝑑θ
𝑀

𝑊
𝑀
=𝑒𝑊 ∗
θ
We have, M = 1.26, W = 4536

∴ IF = 𝑒 ∗ θ
1.26
4536

= 𝑒(2.78∗10^−4) ∗ θ

From equation 11 & 14

y * IF = ∫ 𝑄 (𝑥) IF dx + C

(11) T2A * IF ∫ 𝑄 (θ) IF dθ + C

(14)

A * 𝑒(2.78∗10^−4) ∗ θ = (174 – 86 𝑒−2.78∗10^−4 ∗ θ) * (𝑒2.78∗10^−4 ∗ θ)dθ + C

∫
T2𝑀
𝑊

T2A * 𝑒(2.78∗10^−4) ∗ θ [ 174 * 𝑒2.78∗10^−4 ∗ θ dθ – 86 dθ + C]

𝑀
𝑊
=

- ∫ 86 ∗ dθ + C]
𝑀
T 2A =
174 ∗
𝑒2.78∗10^−4 ∗ θ
𝑊 2.78∗10^−4
[

∫ 𝑒3𝑥 = 𝑒3𝑥/3

Integrating
𝑀
[ 174 * 𝑒2.78∗10^−4 ∗ θ dθ – 86 dθ] + C
𝑊

174 ∗ – 86 θ) + C---------------------------------------------------------(18)
𝑀 (𝑒2.78∗10^−4 ∗ θ
𝑊
=
2.78∗10^−4

Put θ = 0 & T2A = 45°C

= 2.78 ∗ 10^ − 4
𝑀
1.26=
𝑊

4536

45 = 2.78 ∗ 10^ − 4 (
174
Putting in
 2.78∗10^−4

- 0) +C
C = -129
Put C in 18th equation

𝑀 174 ∗ – 86 θ) – 129
𝑊 𝑒
( 2.78∗10^−4 ∗ θ
=
2.78∗10^−4

174 ∗ 𝑒2.78∗10^−4 ∗ −
T2A = 2.78 ∗ 10^ −
129 86 θ
θ −
𝑒2.78∗10 −4𝝷
-
𝑒2.78∗10^−4 ∗ 𝑒2.78∗10
4𝝷
θ∗ 𝑒2.78∗10^−4

86 θ∗ 2.78∗10−4 129
−4𝝷 −4𝝷
𝑒2.78∗10 𝑒2.78∗10
= 174 - -

Put θ = 3600 Sec

86 (3600)∗ 2.78∗10−4 129

−4 (3600) −4 (3600)
𝑒2.78∗10 𝑒2.78∗10
= 174 - -

= 174 – 31.61 – 47.41

= 174 – 79.05

= 94.97 °C
TOPIC: Absorption of CO2 From Air using NaoH Solution

Example 4.10: Some experiments are done on the absorption of air

mixture in 2.5N caustic soda. In one experiment at atmospheric pressure,
the results obtained are as follows (Fig. 4.1 1):

Gas rate at the bottom of the tower = 0.34 kg/m2s

Liquid rate at the top of the tower = 3.94 kg/m 2s

The CO2 in the inlet gas is 315 ppm (parts per million) and in the exit gas
it is 31 ppm. Find the mole fraction of CO2 absorbed by NaOH.

FORMULA:

I. GN+1(YN+1-Y1) = L0(XN-X0)

FIGURE:

X0= 0
y1 = 31 x 10-6
L0 = 3.94 kg/m2s

YN+1 =315 x 10-6

GN+1= 0.34 XN=?
kg/m2s
Given that:

Gas rate at the bottom of the tower = 0.34 kg/m2s

Liquid rate at the top of the tower = 3.94 kg/m 2s

At bottom of tower:

yn+1 = 315x10-6

Gn+1 = 0.34 kg/m2s

Gn+1= 0.34/29 = 0.0117 kmol/m2s (Divide by Air molecular weight)

At the top of the tower:

Y1= 31x10-6

x0 = 0

L0= 3.94 kg/m2s

To Calculate: Mole Fraction of Co2 Absorbed by NaoH (xn)

Solution:

Step:1: Formula

GN+1(YN+1-Y1) = L0(XN-X0)

Step:2

Basis: 1000kg

2.5N NaoH Contains 2.5(40) =100 kg/m3 NaoH

Mean Molecular weight of liquid=

(kg of NaoH× Molecular weight of NaoH) + (kg of water× Molecular

weight of Water)

(Total mass of mixture)

= (100×40) + (900×18)
1000

= (4000+16200)
1000

=20.2 kg/kmol

∴ L0 = 3.94 kg/m2s

20.2 kg/kmol

=0.195 kmol/m2s

By overall Material balance:

=GN+1(YN+1-Y1) = L0(XN-X0)

=0.0117(315×10-6-35×10-6) = 0.195(xN-0)

=0.0117(2.84×10-04) = 0.195 XN

=3.3228×10-6= 0.195 XN

= 3.3228×10-6 = XN
0.195

XN = 17.04×10-6 ppm

CONCLUSION:

The mole fraction of CO₂ absorbed by 2.5N NaOH is 17.04 ppm, showing
limited CO₂ absorption under the given operating conditions. The
calculation is based on an overall material balance (input = output)
principle and verified using consistent gas–liquid flow data.
PROCESS MODELLING SIMULATION AND
OPTIMIZATION

CEC Assignment

SUBMITTED BY
Pratham Panchal (23BT01D029)
Keyur Parmar (23BT01D030)

SUBMITTED TO
Ms. Priyanka Pandya
TOPIC : Unsteady-state Steam Heating of a Liquid
(Dynamic Lumped Parameter Rigid Analytical Model)

Illustration 5.2

A closed kettle (Fig. 5.2) of total surface area A m² is heated through

this surface by condensing steam at temperature Tₛ K. The kettle is
charged with M kg of liquid of heat capacity Cₚ J/kg K at a
temperature of T₀ [Link] the process is controlled by a heat-transfer
coefficient h W/m² K, how does the temperature of the liquid vary
with time?

Fig. 5.2 Solution Steam heating of a liquid in a kettle

In general, Cₚ in J/kg K is referred to as specific heat. In fact, the

specific heat of a fluid (liquid or gas) is defined as the ratio of the
heat capacity of that fluid (C) to the heat capacity of the reference
fluid (Cₚ; i.e., water for liquids and air for gases). Thus, specific heat
is a dimensionless quantity, and heat capacity (Cₚ) has the units of
J/kg K. We follow this convention throughout this book.

Consider a time interval δθ :

Input – Output = Accumulation. (1)
Heat input (J) = h A (Tₛ – T) δθ. (2)
Heat output (J) = 0. (3)
Heat accumulation (J) = M Cₚ (dT / dθ) δθ. (4)
Put (2), (3) ,(4) in equations (1)

hA (Tₛ - T) δθ = M Cₚ (dT / dθ) δθ

On separating the variables, we get

( hA dθ / M Cₚ ) = dT / (Tₛ - T) (5)

On integrating, we get

( hA / M Cₚ ) ∫ dθ = ∫ dT/ ( Tₛ– T)

( hA / M Cₚ ) θ + B = - ln (Tₛ– T) (6)
B can be evaluated by using the boundary conditions, i.e., when θ =
0 , T = T₀ Therefore,

B = - ln ( Tₛ– T₀). (7)

Put (7) in equation (6)

( hA / M Cₚ ) θ – ln (Tₛ - T₀) = - ln (Tₛ - T)

( hA / M Cₚ ) θ – ln (Tₛ - T₀) = - ln (Tₛ - T) + ln (Tₛ - T₀)

( hA / M Cₚ ) θ = ln [(Tₛ -T₀) / (Tₛ - T)]

( hA / M Cₚ ) θ = –ln [(Tₛ -T) / (Tₛ -T₀)]

–( hA / M Cₚ ) θ = ln [(Tₛ -T) / (Tₛ -T₀)]

Exp [-(hA / M Cₚ) θ] = (Tₛ - T) / (Tₛ - T₀)

(Tₛ - T) / (Tₛ - T₀) = exp [–( hA / M Cₚ ) θ]

This is the expression for variation of T with (time).

Example 5.3

A closed vessel of total surface area 40 m² is heated through this

surface by condensing steam at a temperature of 100°[Link] vessel
is charged with 600 kg of liquid having a heat capacity of 2512 J/kg
K at a temperature of 25°[Link] the process is controlled by a heat-
transfer coefficient of 142 W/m² K, model the system in
dimensionless form and find out the temperature of the liquid after
1 hour.
Solution
This is the case of unsteady-state steam heating of a liquid. Hence
the modelling is the same
Given:

Tₛ = 100°C

T₀ = 25°C

Cₚ = 2512 J/kg K
A = 40 m²
M = 600 kg
H = 142 W/m² K

θ = 1 h = 3600 s

(100 – T) / (100 – 25) = exp[-(40 × 142 / (600 × 2512)) × 3600]

(100 – T) / 75 = exp (-13.566878)

(100 – T) / 75 = 1.28226 × 10⁻⁶

100 – T = 75 × 1.28226 × 10⁻⁶

100 – T = 9.5 × 10⁻⁵

T = 100°C

Conclusion

The given illustration and analysis explain the unsteady-state

heating of a liquid in a closed vessel using condensing steam as the
heating source. By applying the principle of heat accumulation, the
temperature change of the liquid with time is modeled using a
lumped capacitance approach.
The energy balance shows that the heat transferred from steam to
the liquid depends on the heat-transfer coefficient (h), surface area
(A), mass of liquid (M), and specific heat capacity (Cₚ). The final
temperature-rise equation obtained is:

(Tₛ - T) / (Tₛ - T₀) = exp [–( hA / M Cₚ ) θ]

This formula predicts how the liquid temperature increases over
time. From the worked example, it is shown that when the heating
surface area, heat-transfer coefficient, and heating time are
sufficiently large, the liquid temperature approaches the steam
temperature (Tₛ). Therefore, after enough time, the liquid becomes
nearly equal in temperature to the condensing steam.
 PROCESS, MODELLING,
SIMULATION & OPTIMIZATION

CEC ASSIGNMENT

SUBMITTED BY
MOHAMMED FARHAN (23BT01D036)
DEV SHIROLAWALA (23BT01D037)

SUBMITTED TO
MS. PRIYANKA PANDYA
 Illustration 5.7 — Unsteady-State Heat Transfer in a
Tubular Gas

A supply of hot air is to be obtained by drawing cool air through a heated

cylindrical preheater. The pipe is 0.1 m in diameter and 1.52 m long, and
is maintained at a temperature of Tw = 316°C throughout its length. The
average values of the properties of the air are as follows: thermal
conductivity k = 0.035 W/m°C, specific heat capacity Cp = 1005 J/kg°C,
overall heat-transfer coefficient U = 4 W/m²/°C, density ρ = 0.801 kg/m³,
flow rate m = 7.87 x 10⁻³ m³/s, inlet temperature To = 31°C, and
temperature rise on heating section ΔT = 52.3°C. x is the distance
measured in metres from the pipe inlet.

Considering a control volume and a heat balance within the gas in the
axial direction, by mass flow of the gas in the axial direction, and by the
above- mentioned thermal conductivity, specific heat, and flow rate, we
develop an analytical mathematical model for the temperature
distribution in the axial direction in the form of a differential equation for
the steady-state condition.

Solution

The system is shown in Fig.

[Diagram: Pipe of length L, 4 in. diameter, inlet temp To, outlet temp TL,
wall temp Tw, and distance x from the inlet.]

Fig. 5.8 Gas preheater

The cross-sectional area A = πD²/4 = 0.00785 m². The heat balance around
the incremental length Δx is as follows:

| | Heat input | Heat output |

1. By conduction | -kA (dT/dx) | -kA * [ (dT/dx) + (d/dx * (-

kA( dT/dx)) * Δx ] |

2. By mass flow | uρCp T | uρCp * [ T + (dT/dx) * Δx ] |

3. Wall heat transfer

| h (πD Δx) (Tw - T) | ----------------|

Input - output = 0 (steady state)

-kA (dT/dx) + uρCp T + h(πD Δx) (Tw - T)

= [ -kA * [ (dT/dx) + (d/dx *(-kA (dT/dx)) * Δx ] + uρCp * [ T

+ (dT/dx) * Δx ] ] = 0 -------- 1

Simplifying equation 1 we get

hπDΔx(Tw-T) = [ (d/dx *(-kA (dT/dx)) + uρCp (dT/dx)]

hπD(Tw-T) = [ (d/dx *(-kA (dT/dx)) * Δx + uρCp (dT/dx) Δx]

hπD(Tw-T) = [(d²T/dx²)-kA]+ uρCp (dT/dx)

Divide by KA and On simplifying, we get

(d²T/dx²) - (uρCp/kA) (dT/dx) + (hπD/kA) (Tw - T) = 0

Putting Tw - T = t and simplifying, we get (d²t/dx²)

- (uρCp/kA) (dt/dx) - (h πD/kA) t = 0
Substituting the numerical values and h = 52.3x-1/2, we get

(d²t/dx²)+ (4*52.3x-1/2/0.035*0.1)(t)-

(7.87*10-3*0.801*1005/0.035*0.00785)dt/dx

(d²t/dx²)+(59771) 3x-1/2(t)-(23058)dt/dx

(d²t/dx²) - (2.31 x 10⁴) (dt/dx) - (5.98 x 10⁴) x-1/2 t = 0 ---- A

Putting z = x² and simplifying, we get

(d²t/dz²) = (1/4z2) (d2t/dz2) - (1/4z3 * dt/dz)

Now Eq. (A) becomes

(1/ 4 z2 )(d2t/dz2) - (1/ 4 z3)(dt/dx) - ( (2.31 x 10⁴) / 2z ) (dt/dx) - (

(5.98 x 10⁴) / z ) t = 0 ----------- (B)

Multiplying eq B by 4z4 can be simplified to give

Z2(d²t/dx²) – z dt/dz + (23.92 x 10⁴ z3) (t) –(46116z3)dt/dz =0

Z2(d²t/dx²) - z ( - (1 + 4.62 x 10⁴ z²) (dt/dz) - (23.92 x 10⁴ z3) t = 0

Process Modeling Simulation and
Optimization

CEC Assignment

Submitted By
Jyoti R. Wagh (23BT01D040)

B. Tech Chemical Engineering

Submitted To

Ms. Priyanka Pandya

GSFC UNIVERSITY, Vadodara

1
Topic: Equation Of Motion

According to Newton’s 2nd Law

F = Force µf
a = Acceleration
gc= Convertion Constant =32.2 lbmft/sec2

1/gc * d/dt (MVi) = ∑𝑛𝑗= 𝑓

1
𝑗𝑖
(1)

Vi = Velocity in the ith direction ft/sec

Fji =jth Force acting in ith Direction

Equation (1) says that the time rate of Change of Momentum in ith direction
is Equal to the net sum of the Force Pushing in the i th direction.
This is basically dynamic Force balance or Conservation of Momentum.
The Gravity Flow tank Shown Below:

Application of Equation of Motion to a Microscopic System. Let the Length

of the exit line be L and its Cross-Sectional Area by A p. The Vertical
Cylindrical tank has Area of A.
The Part of this Process that is Described by a Force Balance is the Liquid
Flowing through the Pipe.
It will have the Mass equal to Volume of Pipe (Ap * L) times Density of
Liquid(ρ)
The Mass of Liquid will have a Velocity of V equal to Volumetric Flow rate
divided by Cross-section area at pipe.
Assume: - (1) Plug Flow Conditions.
(2) Incompressible Liquid.
2
 (3) All the Liquid is Moving at the Same Velocity
(4) Turbulent Flow.
M= Ap *L*ρ V=
F/Ap

The amount of Liquid in the Pipe will not Change with time but for Changing
rate of Outflow the Velocity of Liquid must be Changed and to Change the
Velocity or the Momentum of Liquid use a Force on the Liquid. The
Direction of Interact of Horizontal since the Pipe is assume to be
Horizontal. The Force Pushing on the Liquid at the left end of Pipe is the
Hydraulic Pressure Force of the Liquid in the Tank.

i.e. Hydraulic Force = (A*P*ρ*h*g)/(gc)

The Static Pressure in the Tank at the end of the Pipe are same, hence not
included.

The only Force Pushing in the Opposition direction from Right to Left and
Opposing the Flow is the Frictional Force due to Viscosity of Liquid.
Since we have assumed Flow is Turbulent the Frictional Force will Proportional
to the Velocity and the Length of Pipe.

Frictional Force= kp*L*V2 Friction=

µ*F

From equation (1) (1/gc) * (d/dt (MVi))

1/gc * d/dt* A(ρ*v) =(A*ρ*h*g)/(gc) – (kf*N) (ALρ/gc)

*(du/dt) = (ALρh/gcL) * (g - gc) * (kf)*(N/ALρ)

du/dt = (h/L * g) – *(kf *N*gc)/ (Ap*L*ρ)

Conclusion: The equation of motion derived from Newton’s Second Law

establishes a fundamental relationship between force, mass, and
acceleration. In the context of fluid flow through a pipe connected to a
tank, it demonstrates that the rate of change of momentum of the liquid is
governed by the net forces acting on it — primarily the driving hydraulic
pressure and the opposing frictional resistance. By applying appropriate
assumptions such as incompressible and turbulent flow, the dynamic force
balance helps predict velocity changes and flow rates accurately. Thus, the
equation of motion serves as a key tool in analyzing and optimizing fluid
flow behavior in engineering systems.

3
PROCESS MODELLING SIMULATION AND OPTIMIZATION

CEC Assignment

SUBMITTED BY
Mahek Motiyani (22BT01008)
Sujal Mehta (23BT01D039)

SUBMITTED TO
Ms. Priyanka Pandya
 TOPIC: SIMULATION
Example 11.1

Find
a) the incidence matrix,
b) the digraph, and
c) the associated adjacency matrix for a process represented by a
set of the following equations in functional form:

f1 (x1, x2) = 0
f2 (x4)=0,
f3 (x3, x6) = 0
f4 (x4, x5) = 0
f5 (x1, x6) = 0
f6 (x2, x3, x5) = 0

GIVEN DATA:
f1 (x1, x2) = 0
f2 (x4)=0,
f3 (x3, x6) = 0
f4 (x4, x5) = 0
f5 (x1, x6) = 0
f6 (x2, x3, x5) = 0

TO CALCULATE:
(a) the incidence matrix
(b) the digraph
(c) the associated adjacency matrix

SOLUTION:

(a) Incidence matrix

o In the incidence matrix, the process functions and their

relationship with the process variables are represented as
rows and columns of the matrix, respectively, using the binary
coding (i.e., 1's and 0's).

o Eg. f₁ is a function of x, and x₂ only and so, for row f₁ against

the columns x₁ and x2, ones (1's) are put while against the
remaining columns of the variables X3, X4, X5, and X6 zeros
(0's) are put.

o Similarly, all the other functions and their relationship with

the variables are represented with 1's and O's
 o The corresponding incidence matrix for the above system is as
follows:

X1 X2 X3 X4 X5
X6
f1 [1 10 0 0 0]
f2 [0 00 1 0 0]
f3 [0 01 0 0 1]
f4 [0 00 1 1 0]
f5 [1 00 0 0 1]
f6 [0 11 0 0 0]

(b) Digraph

o In a digraph, the flow of information (i.e., the sequence of

calculations) is given.

o From the equation f2 (x4) = 0 the value of x4, can be

calculated.

o Using this x4, and the equation f4(x4, x5) = 0 the value of x5
can be obtained.

o Using this x5, and with an assumed value of x 3, the equation

f6(x2, x4, x5) = 0 gives the value of x2.

o Using this value of x₂ and the equation f1(x1, x2) = 0 the

value of x₁ could be computed.

o Subsequently, using this value of x₁ in combination with the

equation f5(x1, x6) = 0 the value of x6 could be obtained.

o Finally, using this x6, and the equation f3(x 3, x6) = 0 the value
of x3 is recalculated.

o This value should be equal to the earlier value of x3.

o If not, the iterations are continued till both the values are
same. This entire flow of information is represented in the
form of a digraph.
 X4 X5
2 4 6

X2 X3

1 5 3
X1 X6

As in above we have assumed x3, we can also assume x2 then the

diagraph will change

X4 X5
2 4 6

X3 X2

3 5 1
X6 X1

As can be seen from the figure and the earlier description, the values
of x4 and x5 can be obtained independently, while the values of the
remaining variables x1, x2, x3, and x6 are calculated iteratively.

(c) Associated adjacency matrix

o In the above digraph, the equations are numbered as the

nodes. In the adjacency matrix, rows correspond to the
vertices from which the flow originates and columns
correspond to the vertices to which the flow is directed; it can
be obtained directly from the incidence matrix as well.

o Rules of Boolean multiplication and union apply to these

matrices.

o Going to-------Row
o Coming from-------Column

o The corresponding associated adjacency matrix for the system

can be written using the digraph as

Assumed value x3
 1 2 3 4 5 6
1 [0 0 0 0 1 0]
2 [0 0 0 1 0 0]
3 [0 0 0 0 1 0]
4 [0 0 0 0 1 0]
5 [0 0 1 0 0 0]
6 [1 0 0 0 0 0]

Assumed value x2

1 2 3 4 5 6
1 [ 0 0 0 0 0 1]
2 [ 0 0 0 1 0 0]
3 [ 0 0 0 0 1 0]
4 [ 0 0 0 0 0 1]
5 [ 1 0 0 0 0 0]
6 [ 0 0 1 0 0 0]

When each equation feeds information to the next equation in

sequence and the last one feeds to the first, a loop of information
flow exists, requiring the whole system to be solved simultaneously.
If a set of equations is part of a larger set of equations which itself
forms a larger or maximal loop, the smaller subset must be solved
together with the bigger subset.

CONCLUSION

In this example, the process was analyzed by representing the

given set of functional equations in the form of an incidence matrix,
digraph, and adjacency matrix. The incidence matrix established the
relationship between process variables and functions, helping to
identify dependencies. The digraph illustrated the sequence of
calculations and the information flow among variables, highlighting
which variables could be determined independently and which
required iterative computation. Finally, the adjacency matrix
quantitatively represented these interconnections and feedback
loops.
Through this systematic approach, the analysis demonstrates how
process modeling and simulation can effectively represent complex
system relationships and dependencies, providing a structured basis
for process optimization and computation sequence planning.
PROCESS MODELLING SIMULATION AND OPTIMIZATION CEC
ASSIGNMENT

SUBMITTED BY
22BT01011
22BTD0162

SUBMITTED TO

Ms. Priyanka Pandya

 Example 11.1

Find the Reduced form of the associated adjacency matrix and

diagraph for a process represented by a set of the following
equations in functional form:

f1 (x1, x2) = 0
f2 (x4)=0,
f3 (x3, x6) = 0
f4 (x4, x5) = 0
f5 (x1, x6) = 0 and
f6 (x2, x3, x5) = 0

Given data:

f1 (x1, x2) = 0
f2 (x4) =0,
f3 (x3, x6) = 0
f4 (x4, x5) = 0
f5 (x1, x6) = 0 and
f6 (x2, x3, x5) = 0

To calculate:

Reduced form of the associated adjacency matrix and diagraph

Given the equations:

f1(x1, x2) = 0
f2(x4) = 0
f3(x3, x6) = 0
f4(x4, x5) = 0
f5(x1, x6) = 0

f6(x2, x3, x5) = 0

Step 1: Identifying Nodes and Dependencies

Let’s define: Each function fi as a node for equations/operations.

Each variable xj as a node for system variables.
However, the classic approach is to focus on variables, drawing an edge from
xj to fi if xj appears in fi. For engineering process modeling, we often just
focus on the variables and their direct dependencies.

The set of variables is {x1, x2, x3, x4, x5, x6}

Dependencies:
f1:x1, x2
f2: x4
f3: x3, x6
f4: x4, x5
f5:x1, x6

f6: x2, x3, x5

Edges will be drawn from the variables on which a function depends to the
variable solved by that function.
Step 2: Drawing the Initial Diagraph

Assume each equation solves for a primary variable, e.g., f1 for x1 , f2 for
x2 , etc. This means we associate fi with xi , following a typical functional
dependency graph.

The functional dependencies:

x1 x2 (from f1 )

x2 x4 (from f2 )

x3 x6 (from f3 )

x4 x5 (from f4 )

x5 x1 and x6 (from f5 )

x6 x2 , x3 , x5 (from f6 )

Draw an arrow from a variable that a function depends on to the variable

solved for by that function.

Step 3: Constructing the Adjacency Matrix

Let the order of variables be:x1, x2, x3, x4, x5, x6

The matrix entry a{ij} = 1 if xj appears in the function used to solve for xi

x1 x2 x3 x4 x5 x6
f1 [0 1 0 0 0 0]

f2 [0 0 0 1 0 0]

f3 [0 0 0 0 0 1]

f4 [0 0 0 0 1 0]

f5 [1 0 0 0 0 1]

f6 [0 1 1 0 1 0]
This matrix represents the un-reduced (original) form.

Step 4: Reducing the Matrix and Diagraph

Reduction in this context typically means removing indirect dependencies.

That is, if xi depends on xj via a sequence of dependencies (a path in the
diagraph), and there is also a direct link, remove the direct link (to only keep
minimal paths, similar to a transitive reduction).

x5 x1 and x6 .

x6 x2 , x3 , x5 .

x2 x4 .

x1 x2 .

Step 5: Drawing the Reduced Diagraph

For each non-zero entry in the reduced adjacency matrix, draw a directed
arrow in the diagraph from x_j to x_i (i.e., from column variable to row
variable).

List of edges:

X3 X2

3 5 1

X6 X1
Conclusion: The reduced adjacency matrix and digraph effectively capture
all variable-function dependencies in the process, providing a streamlined
representation that is instrumental in process analysis, system simplification,
and troubleshooting. This reduced form not only clarifies the direct
influences among variables and functions but also enables systematic
approaches for solving, monitoring, and maintaining complex systems.
PROCESS MODELLING SIMULATION AND
OPTIMIZATION

CEC Assignment

SUBMITTED BY
Kuldeep Yadav (23bt01d026)
Hardik Vasava (23bt01d042)

SUBMITTED TO
Ms. Priyanka Pandya
Topic:- Illustration 11.1
Consider a set of equations to be solved:
 (a) x₁x₄ + x₆²/x₄ – 4 = 0
 (b) x₂x₅ – 3x₆ = 0
 (c) x₁/x₂ + ln(x₃/x₄) – 2 = 0
 (d) x₃³ + 2x₃² – 2x₃ = 0
 (e) x₂ + x₄ – 3 = 0
 (f) x₃(x₃ + x₆) – 7 = 0

Solution
There are two algorithms that are put to use while solving the
problem:
1. Rudd’s algorithm: It locates those variables that occur in the
fewest equations.
2. Steward’s algorithm: It isolates equations with the fewest variables
to be solved independently.
After the application of these two algorithms, the rest of the
equations are solved using the SWS algorithm.
At first, it is seen that equation (d) can be solved directly for x₃. With
the value of x₃ obtained from equation (d), equation (f) can be
solved for x₆. By using Rudd’s algorithm, it is seen that x₅ is the
variable that occurs in the fewest number of equations. Therefore,
the set of equations can be precedence ordered such that equation
(b) is solved last. The remaining three equations are then solved
using the SWS algorithm.
The maximal set of equations are arranged one by one to express
one of its variables, the output variable, in terms of the others.

For example, equation (e) is well suited because of its linear form,
thus giving the output variable as x₂:
 x₂ = 3 – x₄

All occurrences of x₂ in the remaining equations (a) and (c) are then
substituted by the above expression. In this manner, one reduces
the number of unknown variables in the maximal set. The remaining
equations are as follows:
 x₁x₄ + x₆²/x₄ – 4 = 0
 x₁/(3 – x₄) + ln(x₃/x₄) – 2 = 0

Now the substituted variable x₂ will be evaluated after the

remaining unknowns x₁ and x₄ in the maximal set are determined.
Equation (11.6) can be rearranged to isolate x₁:
 x₁ = (4 – x₆²/x₄) / x₄

where x₁ is its output variable. The occurrence of x₁ in the remaining

Eq. (11.7) is then substituted by Eq. (11.8), thus giving:
 (4 – x₆²/x₄) / [x₄(3 – x₄)] + ln(x₃/x₄) – 2 = 0

Although one may be able to solve a system of non-linear equations

analytically as illustrated in this illustration, for computer
implementation, the equations have to be solved iteratively by
linearizing the equations (a variant of the Newton–Raphson solution
technique).

Table 11.2 Order of the solution for the given set of

equations (a)–(f)

Order of Equation Output Final form

solution variable
1 (d) x₃ x₃³ – 2x₃ – 2 = 0
2 (f) x₆ x₆ = 7/x₃ – x₃
3 (c) x₄ (4 – x₆²/x₄) / [x₄(3 – x₄)] +
ln(x₃/x₄) – 2 = 0
4 (a) x₁ x₁ = (4 – x₆²/x₄) / x₄
5 (e) x₂ x₂ = 3 – x₄
6 (b) x₅ x₅ = 3x₆/x₂

For preparing a sequential modular steady-state flowsheet, we need

(a) a partitioning algorithm and (b) a tearing algorithm and a
numerical scheme for solving linear equations (which can be used
iteratively to solve non-linear algebraic equations).
For preparing an equation-oriented flowsheet we need (a) an output
set assignment followed by partitioning/precedence ordering and a
numerical method to solve algebraic equations.

Conclusion:-

In this illustration, the systematic approach for solving a set of

nonlinear algebraic equations using process modeling techniques
has been demonstrated effectively. The application of Rudd’s and
Steward’s algorithms provided a structured way to identify
dependent and independent variables, simplify the problem, and
establish a logical order of solution. Through the sequential modular
approach, each equation was reduced and solved step-by-step by
substituting known variables, ultimately leading to an efficient
solution sequence.

The exercise highlights the importance of equation ordering,

variable isolation, and iterative numerical methods—such as the
Newton– Raphson technique—for solving complex systems
encountered in process engineering. This method not only simplifies
the computational process but also enhances accuracy and
convergence in simulation and optimization of chemical processes.
PROCESS MODELLING SIMULATION
AND
OPTIMIZATION

CEC Assignment

SUBMITTED BY
Mitesh Padhiyar ( 22BT01011 )

SUBMITTED TO
Ms. Priyanka Pandya
Illustration 11.1
Consider a set of equations to be solved:
 (a) x₁x₄ + x₆²/x₄ – 4 = 0
 (b) x₂x₅ – 3x₆ = 0
 (c) x₁/x₂ + ln(x₃/x₄) – 2 = 0
 (d) x₃³ + 2x₃² – 2x₃ = 0
 (e) x₂ + x₄ – 3 = 0
 (f) x₃(x₃ + x₆) – 7 = 0

Solution
There are two algorithms that are put to use while solving the
problem:
1. Rudd’s algorithm: It locates those variables that occur in the
fewest equations.
2. Steward’s algorithm: It isolates equations with the fewest variables
to be solved independently.
After the application of these two algorithms, the rest of the
equations are solved using the SWS algorithm.
At first, it is seen that equation (d) can be solved directly for x₃.
With the value of x₃ obtained from equation (d), equation (f) can be
solved for x₆. By using Rudd’s algorithm, it is seen that x₅ is the
variable that occurs in the fewest number of equations. Therefore,
the set of equations can be precedence ordered such that equation
(b) is solved last. The remaining three equations are then solved
using the SWS algorithm.
The maximal set of equations are arranged one by one to express
one of its variables, the output variable, in terms of the others. For
example, equation (e) is well suited because of its

linear form, thus giving the output variable as x₂:

 x₂ = 3 – x₄

All occurrences of x₂ in the remaining equations (a) and (c) are then
substituted by the above expression. In this manner, one reduces the
number of unknown variables in the maximal set. The remaining
equations are as follows:
 x₁x₄ + x₆²/x₄ – 4 = 0
 x₁/(3 – x₄) + ln(x₃/x₄) – 2 = 0

Now the substituted variable x₂ will be evaluated after the

remaining unknowns x₁ and x₄ in the maximal set are determined.
Equation (11.6) can be rearranged to isolate x₁:
 x₁ = (4 – x₆²/x₄) / x₄
where x₁ is its output variable. The occurrence of x₁ in the remaining
Eq. (11.7) is then substituted by Eq. (11.8), thus giving:
 (4 – x₆²/x₄) / [x₄(3 – x₄)] + ln(x₃/x₄) – 2 = 0

Although one may be able to solve a system of non-linear equations

analytically as illustrated in this illustration, for computer
implementation, the equations have to be solved iteratively by
linearizing the equations (a variant of the Newton–Raphson solution
technique).

Table 11.2 Order of the solution for the given set of equations (a)– (f)

Order of Equation Output variable Final form

solution
(d) x₃ x₃³ – 2x₃ – 2 =
0
(f) x₆ x₆ = 7/x₃ – x₃
(c) x₄ (4 – x₆²/x₄) / [x₄(3
– x₄)] +
ln(x₃/x₄) – 2 = 0

(a) x₁ x₁ = (4 –
x₆²/x₄) / x₄
(e) x₂ x₂ = 3 – x₄
(b) x₅ x₅ = 3x₆/x₂

For preparing a sequential modular steady-state flowsheet, we need

(a) a partitioning algorithm and (b) a tearing algorithm and a numerical
scheme for solving linear equations (which can be used iteratively
to solve non-linear algebraic equations).
For preparing an equation-oriented flowsheet we need (a) an output
set assignment followed by partitioning/precedence ordering and a
numerical method to solve algebraic equations.
Process Modelling Simulation and
Optimization

CEC Assignment

Submitted By

RUTU PATEL (22BT01019)

PATEL HETARTH

(22BT01013)

Submitted To

Ms. Priyanka Pandya

GSFC University, Vadodara

LINEAR MATERIAL BALANCE MODELS
Problem Statement
A chemical plant utilizes three raw materials (A, B, and C) to manufacture
three main products (E, F, and G) through three separate processes (1, 2,
and 3).
Consider the process diagram in Figure E2.10.

Each product (E, F, G) requires different (stoichiometric) amounts of

reactants according to the following mass balances:
 Raw material Maximum (xi, kg/day) Cost (c/kg)

A 40,000 1.5

B 30,000 2.0

C 25,000 2.5

Reactant
Processing Selling
requirements
Process Product cost (c/kg price
(kg/kg
product) (c/kg)
product)

1 E 0.667 A, 0.333 B 0.5 15.0

2 F 0.667 A, 0.333 B 0.5 16.0

0.5 A, 0.167 B,
3 G 1.0 18.0
0.333 C

Formulate a complete linear programming model to determine the optimal

production mix. The model must include:
1. The definitions of all variables.
2. All material balance and supply constraints.
3. The final objective function to maximize the total daily operating
profit.

Solution
Formula Used
1. Mass Balance:
Total mass of reactant consumed=∑Reactant coefficient * Product
mass flow.
 2. Objective Function: Maximize Z = Total Revenue - Total Operating
Costs
 Total Operating Costs=Raw Material Costs + Processing Costs

We first define the system's variables and the constraints that govern
their relationships.
Given that,
The process involves 12 mass flow variables X1 to X12 and is governed by
the relationships and data provided in the tables above. The operating
profit must be maximized.

1. Variables Definitions
We define three sets of variables: the final product flow rates (our main
decision variables) and the total raw material flow rates (which are
dependent on the product flows) and Other sub category.
 Production Variables:
o X8: Daily production rate of Product E (kg/day)
o X9: Daily production rate of Product F (kg/day)
o X10: Daily production rate of Product G (kg/day)
 Raw Material Variables:
o X11: Total daily consumption of Reactant A (kg/day)
o X12: Total daily consumption of Reactant B (kg/day)
o X7: Total daily consumption of Reactant C (kg/day)

 Other:
o X1 :Mass flow of Reactant A entering Process 1(for E).
o X2 :Mass flow of Reactant A entering Process 2(for F).
o X3:Mass flow of Reactant A entering Process 3(for G).
o X4: Mass flow of Reactant B entering Process 1(for E).
o X5: Mass flow of Reactant B entering Process 2(for E).
o X6: Mass flow of Reactant B entering Process3(for E).
 2. Model Constraints
Next, we establish the constraints that govern the system. These include the
mass balance equations and the supply.
 Mass Balance Equations (Equality Constraints): These link the
reactant consumption to the product production rates, based on the
"Reactant requirements" table.
a) X11= 0.667X1+ 0.667X2+ 0.5X3
b) X12= 0.333X4+ 0.333X5+ 0.167X6
c) X11 + X10 + X7 = X8 + X9 + X10

 Supply Constraints (Inequality Constraints): The total consumption

of each reactant cannot exceed its maximum daily availability.
a) X11 ≤ 40,000(Kg/day)
b) X12 ≤30,000(Kg/day)
c) X7 ≤ 25,000(Kg/day)

 Non-Negativity Constraints: The production of each product cannot

be negative.
a) X8 ≥ 20,000(Kg/day)
b) X9 ≥ 25,000(Kg/day)
c) X10 ≥ 30,000(Kg/day)

3. Maximize Profit
Determining Net Profit Contribution (Objective Function components)
A. Income from Selling Price
Selling Price = 0.04 E + 0.033 F + 0.038 G
B. Total Operating
Costs Raw Material
costs
 CostRM = 0.015 A + 0.02 B + 0.025 C
Processing Costs
 CostProcess=0.015 E + 0.005 F + 0.01 G

TOTAL COST=0.015A + 0.02B + 0.025C + 0.015E + 0.005F + 0.01G

Therefore,
Daily Profit = Selling Price – (Raw Material costs + Processing Costs)
= 0.04 E + 0.033 F + 0.038 G - (0.015A + 0.02B + 0.025C + 0.015E
+ 0.005F + 0.01G)
f(x)=0.025 E + 0.028F + 0.028G – 0.015A – 0.02B – 0.025C
=0.025X8 + 0.028X9 + 0.028X10 – 0.015X11 – 0.02X12 -0.025X7
Note that the six variables in the objective function are constrained through
material balances, namely
X11= 0.667X8 + 0.667X9 + 0.5X10

X12 = 0.333X8 + 0.333X9 + 0.167X10

X7 = 0.333X10

Also,
0 ≤ X11 ≤ 40,000

0≤ X12 ≤ 30,000

0≤ X7 ≤ 25,000

Conclusion
This formulation successfully combines the material balance structure
(determining coefficients based on required inputs) with the economic
objectives (maximizing net profit). The final result is a standard Linear
Programming problem with three variables and three functional
constraints.
PROCESS MODELLING SIMULATION AND
OPTIMIZATION

CEC Assignment

SUBMITTED BY
Prince Patel (22BT01014)
Sattwik Sharma (22BT01020)

SUBMITTED TO
Ms. Priyanka Pandya
Sr. Assistant Professor,
Dept. of Chemical Engineering
GSFC University
Modeling and Optimizing Blast Furnace Operation
Optimizing the operation of the blast furnace is important in every
large- scale steel mill. A relatively large number of important variables
(several of which cannot be measured) interact in this process in a
highly complex manner numerous constraints must be taken into
account, and the age and efficiency of the plant significantly affect the
optimum operating point (Deitz, 1997).
Consequently, a detailed examination of this problem demonstrates the
considerations involved in mathematical modeling of a typical process.
The operation of a blast furnace is semicontinuous. The raw materials
are iron ore containing roughly 20 to 60 percent iron as oxides and a
variety of other metallic and nonmetallic oxides. These materials are
combined with coke, which reacts to form blast furnace gas. Limestone
is a flux that helps separate the impurities from the hot metal by
influencing the pH. Apart from the blast furnace gas, which may serve
as a heating medium in other processes, the output of the furnace
consists of molten iron, which includes some impurities (notably carbon
and phosphorus) that must be removed in the steel making process,
and slag, which contains most of the impurities and is of little value.

Operation of the blast furnace calls for determination of the amount of

each ore, a production rate, and a mode of operation that will maximize
the difference between the product value and the cost of producing the
required quantity and quality of molten iron. Figure shows the flow of
materials in the blast furnace, which itself is part of a much larger mill.
One ton of hot metal requires about 1.7 tons of iron-bearing materials,
0.5 to 0.65 tons of coke and other fuel, 0.25 tons of fluxes, and 1.8 to
2.0 tons of air. In addition, for each ton of hot metal produced, the
process creates 0.2 to 0.4 tons of slag, 0.05 tons or less of flue dust,
and 2.5 to 3.5 tons of blast furnace gases. The final product, hot metal,
is about 93% iron, with other trace ingredients, including sulfur, silicon,
phosphorus, and manganese.

The process variables and conceptual models are identified in Figure

under the column “Process Analysis,” which has categories for the
objective function,
equality constraints, and inequality constraints.

Objective Function:

To formulate the objective function, two categories of costs have to be

considered:
1. Costs associated with the material flows (the input and
output variables), such as the costs of purchased materials.
2. Costs associated with the operations related to the process variables
in the model.
The terms that make up the objective function (to be maximized) are
shown in Figure The profit of the blast furnace can be expressed as:

f = Σ (ci * xi) - Σ (cj * xj)

Equality and Inequality Constraints:

1. Iron ore: Ores of different grades are available in restricted quantities.

Different ores have varying percentages of iron and different types
and amounts of impurities. The proportion of each ore that occurs in the
final hot metal is assumed to be fixed by its composition.
2. Coke: The amount of coke that may be burned in any furnace is
effectively limited by the furnace design, and the hot metal temperature
is controlled by the amount of coke.
3. Slag: For technical reasons, the level of impurities in the slag must
be controlled. There is an upper limit on the percentage of magnesium,
upper and lower limits on the percentage of silicon and aluminum, and
close limits on the “basicity” ratio (CaO + MgO)/ (SiO2 + Al2O3).
4. Phosphorus: All phosphorus in the raw material finds its way into
the molten metal. There is an upper limit on the phosphorus
permitted.

The basicity ratio can be expressed as:

(Σ w2i*xi + Σ w3i*xi) / (Σ w4i*xi + Σ w5i*xi) where:

w2i = weight fraction of CaO in feed

w3i = weight fraction of MgO in feed
w4i = weight fraction of SiO2 in feed
w5i = weight fraction of Al2O3 in feed

Preparation Steps:

1. A set of input and output variables.

2. A set of steady-state input–output material and energy
balances (equality constraints).
3. A set of explicit empirical relations (equality constraints).
4. Inequality constraints for raw materials and process variables.
5. An objective function to maximize the net profit of the furnace.
Material Balance Reconciliation
Question:

Suppose the flow rates entering and leaving a process are measured
periodically. Determine the best value for stream A (in kg/h) for the
process
shown from the three hourly measurements indicated of B and C in Figure,
assuming steady-state operation at a fixed operating point.

The process model is:

MA + MC = MB
where M is the mass per unit time of throughput.

Solution:

We need to set up the objective function first. Let us minimize the sum of
the squares of the deviations between input and output as the criterion
so that the objective function becomes:

f(MA) = (MA + 11.1 - 92.4) ^2 + (MA + 10.8 - 94.3) ^2 + (MA + 11.4 -

93.8) ^2

A sum of squares is used since this guarantees that f > 0 for all values of
MA; a minimum at f
= 0 implies no error.

Given Measurements in Figure:

For C:
(a) 11.1 kg/h
(b) 10.8 kg/h
(c) 11.4 kg/h

For B:
(a) 92.4 kg/h
(b) 94.3 kg/h
(c) 93.8 kg/h
Calculation and Answer:

The optimum value of MA can be found by differentiating f with respect to

MA. This leads to an optimum value for MA of 82.4 kg/h,
which is also obtained by computing the averaged measured values:

MA = MB - MC

Final Answer: MA =
82.4 kg/h
Process Modelling Simulation and Optimization
CEC Component

Submitted By

Manan Pathak (22BT01015)

Mann Chavda (22BT01006)

Submitted To

Ms. Priyanka Pandya, Sr. Assistant Professor, Chemical

Engineering GSFC University, Vadodara

1
 OPTIMIZATION OF A THERMAL CRACKER VIA LINEAR
PROGRAMMING

Reactor systems that can be described by a "yield matrix" are potential

candidates for the application of linear programming. In these situations,
each reactant is known to produce a certain distribution of products.
When multiple reactants are employed, it is desirable to optimize the
amounts of each reactant so that the products satisfy flow and demand
constraints. Linear programming has become widely adopted in
scheduling production in olefin units and catalytic crackers. In this
example, we illustrate the use of linear programming to optimize the
operation of a thermal cracker sketched in Figure:1.

Figure:1 Flow diagram of thermal cracker.

Table:1 shows various feeds and the corresponding product distribution

for a thermal cracker that produces olefins. The possible feeds include
ethane, propane, debutanized natural gasoline (DNG), and gas oil, some
of which may be fed simultaneously. Based on plant data, eight products
are produced in varying proportions according to the following matrix. The
capacity to run gas feeds through the cracker is 200,000 lb/stream hour
(total flow based on an average mixture). Ethane uses the equivalent of
1.1 lb of capacity per pound of ethane; propane 0.9 lb; gas oil 0.9 lb/lb and
DNG 1.0.

2
 Feed
Product Ethane Propane Gas oil DNG
Methane 0.07 0.25 0.10 0.15
Ethane 0.40 0.06 0.04 0.05
Ethylene 0.50 0.35 0.20 0.25
Propane - 0.10 0.01 0.01
Propylene 0.01 0.15 0.15 0.18
Butadiene 0.01 0.02 0.04 0.05
Gasoline 0.01 0.07 0.25 0.30
Fuel Oil - - 0.21 0.01
Table:1 Yield structure (wt. fraction)

Downstream processing limits exist of 50,000 lb/stream hour on the

ethylene and 20,000 lb/stream hour on the propylene. The fuel
requirements to run the cracking system for each feedstock type are as
follows:

Feedstock type Fuel requirement (Btu/lb)

Ethane 8364

Propane 5016

Gas oil 3900

DNG 4553

Methane and fuel oil produced by the cracker are recycled as fuel. All the
ethane and propane produced is recycled as feed. Heating values are as
follows:
Recycled feed Heat produced (Btu/lb)

Natural gas 21,520

Methane 21,520

Fuel oil 18,000

3
Because of heat losses and the energy requirements for pyrolysis, the
fixed fuel requirement is 20.0 * 106 Btu/stream hour. The price structure
on the feeds and products and fuel costs is:

Feeds Price ($/lb)

Ethane 6.55
Propane 9.73
Gas oil 12.50
DNG 10.14

Products Price ($/lb)

Methane 5.38 (fuel value)
Ethylene 17.75
Propylene 13.79
Butadiene 26.64
Gasoline 9.93
Fuel oil 4.50 (fuel value)

 Assume an energy (fuel) cost of $2.50/106 Btu.

4
Formula to be used:

 Objective (profit) (in words):

Profit = Product value − Feed cost − Energy cost

 In equation form (summing dollar contributions per hour):

Profit f = ∑ (price * production) - ∑ (price * fresh feed) - fuel

cost
products feeds

We substitute production as linear functions of the feed flows using the

yield fractions from the yield matrix. The energy (fuel) cost comes from
fuel purchased to meet the cracking energy requirement after accounting
for recycled fuel produced from cracker (methane and fuel oil). The fixed
heat loss (20×106 Btu/h) can be treated as constant (affects final
accounting but not LP decision if it’s truly constant), while the variable
part depends on feed mix and recycled fuel.

Figure:

5
Given that:

1. Decision (unknown) variables (defined below).

2. Cracker capacity (expressed as equivalent lb feed): 200,000

lb/stream*h (total capacity measured using equivalent weights:
ethane = 1.1 lb capacity per lb, propane = 0.9, gas oil = 0.9, DNG =
1.0).

3. Downstream product processing limits:

a. ethylene ≤ 50,000 lb/h.
b. Propylene ≤ 20,000 lb/h.

4. Fuel requirement to run cracking (fixed + feed-dependent): fixed

heat loss = 20.0 × 106 Btu/hr. Fuel consumption per feed (Btu per lb
feed):
a. Ethane: 8364 Btu/lb
b. Propane: 5016 Btu/lb
c. Gas oil: 3900 Btu/lb
d. DNG: 4553 Btu/lb.

5. Heating values (Btu/lb): methane (and natural gas) 21,520, fuel oil
18,000. Methane and fuel oil produced by cracker are recycled as
fuel.

6. Prices (per lb):

a. Feeds — ethane $6.55, propane $9.73, gas oil $12.50, DNG
$10.14.
b. Products — methane $5.38 (fuel value), ethylene $17.75,
propylene $13.79, butadiene $26.64, gasoline $9.93, fuel oil
$4.50 (fuel).
c. Energy (fuel) cost assumed $2.50 per 106 Btu.

6
 7. Yield (weight-fraction) matrix for products from each feed in
fractions are mentioned.

To Calculate:
 Objective function (profit):

f = Product value - Feed cost - Energy cost

 Optimal feed rates: x1, x2, x3, x4, x5, x6, x7.

 Maximum profit per hour ($/h) under the given operational conditions.

 To produce the LP formulation, identify design/decision variables,

write objective and constraints, and state process model (mass &
energy balances).

Solution:
The optimization process involves the following procedure;
 Set up the objective function and constraints to maximize profit while
operating within furnace and downstream process equipment
constraints. The variables to be optimized are the amounts of the
four feeds.

 Solve using linear programming.

 Examine the sensitivity of profits to increases in the ethylene

production rate.

7
Step 1: Formula used:

Profit = Product value − Feed cost − Energy cost

i.e.,

Profit f= ∑ (price * production) - ∑ (price * fresh feed ) - fuel cost

products feeds

Step 2: The design or decision variable:

We define the following variables for the flow rates to and from the
furnace (in lb/h):

x1 = fresh ethane feed

x2 = fresh propane feed
x3 = gas oil feed
x4 = DNG feed
x5 = ethane recycle
x6 = propane recycle
x7 = fuel added
All xi ≥ 0

 Assumptions used in formulating the objective function and

constraints are
I. 20*106 Btu/h fixed fuel requirement (methane) to
compensate for the heat loss.
II. All propane and ethane are recycled with the feed, and all
methane and fuel oil are recycled as fuel.

 A basis of 1 hour is used, and all costs are calculated in cents per
hour.

8
Step 3: Objective function (profit). In words, the profit f is

f = Product value - Feed cost - Energy cost

Product value:
The value for each product (in cents per pound) is as follows:

Ethylene: 17.75 * (0.5x1 + 0.5x5 + 0.35x2 + 0.35x6 + 0.20x3 + 0.25x4)

(a)

Propylene: 13.79 * (0.01x1 + 0.01x5 + 0.15x2 + 0.15x6 + 0.15x3 +

0.18x4) (b)

Butadiene: 26.64 * (0.01x1 + 0.01x5 + 0.02x2 + 0.02x6 + 0.04x3 +

0.05x4) (c)

Gasoline: 9.93 * (0.5x1 + 0.5x5 + 0.35x2 + 0.35x6 + 0.25x3 + 0.30x4) _(d)

Total product sales = 9.39x1 + 9.51x2 + 9.17x3 + 11.23x4 + 9.39x5 +

9.51x6 _(e)

Feed cost:
Feed cost ($/h) = 6.55x1 + 9.73 x2 + 12.50x3 + 10.14x4 (f)

Energy cost:

The fixed heat loss of 20 * l06 Btu/h can be expressed in terms of methane
cost (5.38$/lb) using a heating value of 21,520 Btu/lb for methane. The
fixed heat loss represents a constant cost that is independent of the
variables xi, hence in optimization we can ignore this factor, but in
evaluating the final costs this term must be taken into account. The value

9
for x7 depends on the amount of fuel oil and methane produced in the
cracker (x7 provides for any deficit in products recycled as fuel).

We combine (e) and (f) to get the objective function ($/h)

f = 2.84x1 – 0.22x2 – 3.33x3 + 1.09x4 + 9.39x5 + 9.51x6

(g)

Step 4: Constraints:
1. Cracker capacity of 200,000 lb/h.

1.1(x1 + x5) + 0.9(x2 + x6) + 0.9x3 + 1.0x4 ≤ 200,000 (h)

(or)

1.1x1 + 0.9x2 + 0.9x3 + 1.0x4 + 1.1x5 + 0.9x6 ≤ 200,000

2. Ethylene processing limitation of 100,000 lb/h.

0.5x1 + 0.35x2 + 0.25x3 + 0.25x4 + 0.5x5 + 0.35x6 ≤ 100,000

(i)

3. Propylene processing limitation of 20,000 lb/h.

0.01x1 + 0.15x2 + 0.15x3 + 0.18x4 + 0.01x5 + 0.15x6 ≤ 20,000

(j)

4. Ethane recycle.

x5 = 0.4x1 + 0.4x5 + 0.06x2 + 0.06x6 + 0.04x3 + 0.05x4 (k)

Rearranging, (j) becomes

0.4x1 + 0.06x2 + 0.04x3 + 0.05x4 - 0.6x5 + 0.06x6 = 0 (l)

5. Propane recycle.

X6 = 0.1x2 + 0.1x2 + 0.01x3 + 0.01x4 (m)

10
 Rearranging Equation (m),

0.1x2 + 0.01x3 + 0.01x3 - 0.9x6 = 0 (n)

6. Heat constraint.
The total fuel heating value (THV) (in Btu/h) is given by;

(fuel) (methane from cracker)

THV = 21,520x7 + 21,520 * (0.07x1 + 0.25x2 + 0.10x3 + 0.15x4

- 0.07x5 + 0.25x6) + 18,0000 *(0.21x3 + 0.01x4)

(fuel oil from cracker)

THV = 1506.4x1 + 5380x2 + 5932x3 + 3408x4 + 1506.4x5 +

5308x6 + 21,520x7 (o)

The required fuel for cracking (Btu/h) is

8364(x1 + x5) + 5016(x2 + x6) + 3900x3 + 4553x4
(ethane) (Propane) (Gas oil) (DNG)

= 83644x1 + 5016x2 + 3900x3 + 3408x4 + 8364x5 +

5016x6 (p)

Therefore, the sum of Equation (p) + 20,000,000 Btu/h is equal to the

THV from Equation (o), which gives the constraint;

-6857.6x1 + 364x2 + 2032x3 - 1145x4 – 6857.6x5 + 364x6 +21,520x7 =

20,000,000

(q)

7. Non-negativity: xi ≥ 0 for all i.

All constraints above are linear (equalities or inequalities), so the overall

problem is a Linear Program (LP).

11
Table:2 lists the optimal solution of this problem obtained using the Excel
Solver (case I). Note that the maximum amount of ethylene is produced.
As the ethylene production constraint is relaxed, the objective function
value increases. Once the constraint is raised above 90,909 lb/h, the
objective function remains constant.

Flow rate (lb/h)

Stream Case 1 Case 2
x1 (ethane feed) 60,000 21,770
x2 (propane feed) 0 0
x3 (gas oil feed) 0 0
x4 (DNG feed) 0 107,600
x5 (ethane recycle) 40,000 23,600
x6 (propane recycle) 0 1,195
x7 (fuel added) 32,800 21,090
Ethylene 50,000 50,000
Propylene 1,000 20,000
Butadiene 1,000 5,857
Gasoline 1,000 32,820
Methane (recycled to
7,000 19,610
fuel)
Fuel oil 0 1,076
Objective function
369,560 298,590
($/h)
Table:2 Optimal flow rates for cracking furnace for different restrictions on ethylene
and propylene production

Suppose the inequality constraints on ethylene and propylene production

were changed to equality constraints (ethylene = 50,000; propylene =
20,000). The optimal solution for these conditions is shown as case 2 in
Table:2. This specification forces the use of DNG as well as ethane.

12
Step 5: Process model

 Mass balances (component/product yields):

For each product p (methane, ethane, ethylene, propane,

propylene, butadiene, gasoline, fuel oil) the production flow
(lb/h) = sum over feeds of (yield fraction of product p from
feed j) * (feed flow j).

Some product streams (ethane, propane) are routed as

recycle to feed (so recycle variables appear); methane and
fuel oil are recycled as fuel. These linear mass balances
produce the product expressions used in objective/constraints.

 Energy (heat) balance: the required Btu/hr to vaporize/crack the feed

(feed-dependent + fixed losses) must be supplied by the heating
value of purchased fuel plus heating from recycled methane and
fuel oil. This gives the linear heat constraint above.

Conclusion
 Optimal Strategy: The optimal strategy is to prioritize the use of the
most profitable feeds (Gas Oil and Ethane, which have the highest
objective function coefficients) while completely shutting off the less
profitable feeds (Propane and DNG).

 Binding Constraints: The optimal solution is severely constrained by three

factors, which are simultaneously active:
 Ethane Supply: The maximum available ethane is used.
 Cracker Capacity: The cracker is run at its 100% capacity limit.
 Propylene Demand: The production of propylene is exactly met
at the minimum required amount.

 Economic Interpretation: The shadow prices of the binding constraints

are critical. The solution suggests that increasing the availability of
ethane or increasing the cracker's total capacity would lead to a
significant increase in overall profit. The shadow price on the
Propylene Demand constraint indicates the profit loss incurred by
having to meet that minimum target.

 Methodology: The use of Linear Programming is appropriate here

because the system can be accurately described by a "yield matrix"
with linear objective and constraint functions.

13
PROCESS MODELING SIMULATION &
OPTIMIZATION

Submitted By
HARSH JAYSWAL (22BT01004)
TIRTH UPADHYAY
(22BT01023) CHEMICAL
ENGINEERING

Submitted To
Mrs. Priyanka Pandya
GSFC University, Vadodara

Date of Submission: November 7th, 2025

  Example : 12.4 :

DETERMINATION OF THE OPTIMAL REFLUX RATIO FOR A STAGED-

DISTILLATION COLUMN.

 Solution :

1. Design variables-

i. R- reflux ratio -continuous { R >= R min }

ii. N- total no. of ideal stages { N >= N min }
iii. F – feed stage location
iv. D- distillate flow
-Generally, the optimization is done over R & N, feed stage and distillate
are secondary.

2. Objective function-
-To maximize the profit of the distillation column operation.
Maximize Z = (D . Price ) – Steam cost (R) – Cooling cost (R)
i.e :
f = Propylene sales + Propane sales - Utility costs - Raw material costs

f =(CDXDD + CBXBB) + [C′ (1 − XD)D + C′ (1 − XB)B]

D B

−[𝐶1𝑄𝑅 + 𝐶2𝑄𝐶] − [𝐶𝐹𝑋𝐹𝐹 +�𝐶′ (1 − 𝑋𝐹)𝐹]

� ……eq.(i)

-Where, QR = Reboiler heat requirement

QC = Cooling load
-From eq. (e) we have : XF F = XD D + XB B ….substitute in (eq. i )
-And defining : { -W = CB – CD } & { -U= C’D – C’B },
-Therefore,
f = CDXFF + C’B ( 1 - XF )F – CFXF F – C’F ( 1 - X )F –
C1QR -C2QC - WXBB- U( 1-XD )D
-& QR = QC = Λv

-Therefore,
f1 = ( C1 + C2 ) Λv + WXBB + U(1 – XD )D

3. Constraints-

-Calculate Rm: The formula is given as:

1 𝑋𝐷 (1 − 𝑋𝐷)
𝑅 =
𝑚 [
(𝛼 − 1) 𝑋𝐹 − 𝛼
(1 − ]
𝑋𝐹)
……eq. a

𝑁 − 𝑁𝑚 𝑅 − 𝑅𝑚 0.5668
-Substitute the value of Rm in eq. b to find Nm,

= 0.75[1 − ( ) ]
𝑁+1 𝑅+
1
……eq. b

-Now compute Xb from,

……eq. c

-Equations (a) & (c) comprise equality constraints relating X B and R.

 - where, XB = product purity
R = Recovery

4. Process Model-

The process model is the set of relationships, primarily the mass &
component balances & the design corelations that describe the operation
of the distillation column.

Mass & Component balances:

i. Overall mass balance {F=D+B}

ii. Overall component balance { XFF = XDD + XBB }
iii. Liquid flow rate { L = R.D }
iv. Vapour flow rate { V = (R+1) D }
The (a), (b) & (C) equations of the constraints are the core relationships.
used to calculate dependent variables like XD based on the design Variable
R.
Final result :
i. The value of Rm is 11.338 …….from eq.(a)
ii. The optimum reflux ratio is 17.06 ……from table 12.4 b
iii. The cost f1 is $3870/day ……from table 12.4 b
PMSO CEC: OPTIMAL DESIGN AND OPERATION
OF A CONVENTIONAL STAGED-DISTILLATION
COLUMN

Submitted by:
Mr. Aryarajsinh Gohel (22BT01001) Mr.
Aadarsh Shah (22BT01021)

Submitted To:
Ms. Priyanka Pandya

Subject:
Process Modelling, Simulation and Optimization
Problem 12.1:
OPTIMAL DESIGN AND OPERATION OF A CONVENTIONAL STAGED-
DISTILLATION COLUMN

Theory:
Distillation is the most commonly used industrial separation process, and
its optimization problems vary in complexity.
1. Basic Optimization: The simplest case involves finding optimal
operating conditions such as heat input, cooling rates, and product
flow rates for an existing column to minimize cost or energy use
while meeting specific purity or recovery requirements. This is
typically a nonlinear programming (NLP) problem. When all degrees
of freedom are fixed, it reduces to a standard distillation design
problem.
2. Design Optimization: A more advanced problem includes
determining both the optimal operating conditions and the minimum
number of stages. Since the number of stages is a discrete variable,
this becomes a mixed-integer nonlinear programming (MINLP)
problem, balancing capital costs (increasing with stages) and
operating costs (decreasing with stages).
3. Comprehensive Optimization: The most complex task involves
optimizing the number of stages along with the feed and side
stream locations. However, these locations usually fall within a
narrow practical range, making the objective function relatively
insensitive to small variations.

Question:
Determine the optimal feed stage locations and reboiler heat duty for a
four-stage distillation column that minimizes the cost of operation, which
is directly proportional to the reboiler heat duty Q 1. The column separates
a five-component hydrocarbon feed. The specifications are as follows:

 Total Feed: 100 lb mol/h (liquid)

 Feed Enthalpy, hf : 4000 Btu/lb mol
 Feed Composition:
x1 (C₃H₈) = 0.05
x2 (i-C₄H₁₀) = 0.15 x3
(n-C₄H₁₀) = 0.25 x4
(i-C₅H₁₂) = 0.20 x5
(n-C₅H₁₂) = 0.35
  Top Product Flow Rate: 10 lb mol/h (liquid)
 Product Purity Constraint: Mole fraction of n-C₅H₁₂ in the top product
(x5,4) ≤ 0.07
 The decision variables are the reboiler heat duty Q1 and the feed
flow rates to each of the 4 stages (F1, F2, F3, F4), with the sum of all
feeds being 100 lb mol/h.

Figure:

The example focuses on the design and optimization of a steady-state

staged distillation column. A typical column setup is illustrated in Figure
E12.1, with variables and parameters defined in Table E12.1A. The
notation uses superscript F for the feed, subscripts indicating withdrawal
stages, and superscripts V and L to denote the vapor and liquid phases,
respectively.
The column stages are numbered from the bottom (reboiler) upward, with
k = 1 representing the first stage. At the bottom, vapor flow V 0 and liquid
flow L1 are zero, while at the top (condenser), vapor V n and liquid Ln+1 are
zero.
The formulation of the optimization problem proceeds in three parts:
1. Equality constraints – representing mass and energy balances,
phase equilibria, and summation conditions.
2. Inequality constraints – capturing process limits, specifications, and
operating restrictions.
3. Objective function – typically aimed at minimizing total cost or
energy consumption while maintaining desired separation
performance.

Objective Function:
The main costs of operation are the heating and cooling costs that are
related to Q, and Q,, respectively. We assume all the other values of Q,
are zero. Q, is determined from the energy balance, so that Q, is the
independent variable. The cost of operation per annum is assumed to be
directly proportional to Q, because the maintenance and cooling costs are
relatively small and the capital costs per annum are already fixed.
Consequently, the objective function is relatively simple:
Minimize: Q1

Formulas to be used:

FL + FV + VK-1 + Lk+1 = Vk + Lk + WV + WL -(a)

1. Total Material Balance:

K K K K
x FL + y F + yi,K-1VK-1 + xi,K+1LK+1 = yi,kVk + xi,kLk + yi,kW + xi,kWL -(b)
2. Component Material Balance:
F F V V
i,k K i,k K K K

Qk +hF Fk+ Hk-lVk-l + hk+lLk+l= HkVk + hkLk + HkWV+ hkWL -(c)

3. Energy Balance:

K K K

4. Equilibrium relations for liquid and vapor at each stage (one for each
stage):
Yi,k=Ki,kxi,k -(d)

5. Relation between equilibrium constant and p, T, x, y (one for each

stage):
Ki,k = Ki(pk,Tk,xk,yk) -(e)

6. Relation between enthalpies and p, T, x, y (one for each stage):

hk = h(pk,Tk,xk) -(f)

Hk = H(pk,Tk,yk) -(g)

7. Key Inequality Constraints:

 Q1 ≥ 0
 Q4 ≤ 0 (Condenser duty)
 xi,k ≥ 0, yi,k ≥ 0
 Fk ≥ 0
 x5,4 ≤ 0.07

The preceding classic set of algebraic equations form a well-defined

sparse structure that has been analyzed extensively. Innumerable
techniques of solution have been proposed for problems with 0 degrees of
freedom, that is, the column operating or design variables are completely
specified. Our interest here in posing an optimization problem is to have
one or more degrees of freedom left after prespecifying the values of
most of the independent variables. Frequently, values are given for the
following parameters:
(a) Number of stages
(b) Flow rate, composition, and enthalpy of the feed(s)
(c) Location of the feed(s) and side stream withdrawal(s)
(d) Flow rate of the side strearn(s)
(e) Heat input rate to each stage except one
(f) Stage pressures (based on column detailed design specifications)
Given That:
 Number of stages = 4 (Stage 1: Reboiler, Stage 4: Condenser)
 Feed and product specifications as listed in the question.
 The problem has 5 degrees of freedom (Q1, F1, F2, F3, F4).
 Initial guesses for the variables are required to start the
optimization algorithm.

To Calculate: The optimal values of the decision variables (Q1, F1, F2, F3, F4)
that minimize Q1 while satisfying all constraints.

Solution
Step 1: Problem Formulation
The optimization problem is formally stated as:

Minimize Q1 with subject to

∑4𝑘= 𝐹k =100
1

Distillation Model Equations (a-i):

Q1≥0

Fk≥0, k=1,

…,4 x5,4 ≤0.07

Step 2:
Key Inequality Constraints:

Q1 ≥ 0
Q4 ≤ 0 (Condenser duty)
xi,k ≥ 0, yi,k ≥ 0
Fk ≥ 0
x5,4 ≤ 0.07

Initialization:
A reasonable initial guess, assuming an equal split of the feed, is used to
start the numerical solver:
F1 = 25

F2 = 25

F3 = 25
 F4 = 25

Q1 = 5 * 106

Step 3: Numerical Solution Method: A Sequential Quadratic Programming

(SQP) algorithm is used to solve this NLP problem. The solution involves
nested calculations:

minimize 𝑄1.
 Outer Loop: The SQP algorithm adjusts the decision variables to

 Inner Loop: For each set of variables proposed by the SQP, the full
set of distillation column equations is solved to determine all
internal states (temperatures, flow rates, compositions) and check
constraints. Physical properties (enthalpies, equilibrium constants)
are evaluated in this inner loop.

Step 4: Optimization Results:

The SQP algorithm converges to the following optimal solution:

Variable Description Initial Guess Optimal Value

F1 Feed to Stage 25 lb mol/h 23.7 lb mol/h
1 (Reboiler)
F2 Feed to Stage 25 lb mol/h 23.7 lb mol/h
2
F3 Feed to Stage 25 lb mol/h 23.7 lb mol/h
3
F4 Feed to Stage 25 lb mol/h 23.7 lb mol/h
4 (Condenser)
Q1 Reboiler 5 * 106 BTU 3.38 * 106 BTU
Heat
Duty
x5,4 n-C₅H₁₂ in Top -- 0.07
Product

Final Answer:
The optimal configuration that minimizes operating cost is:

Q1 = 3.38 * 106 BTU

Feed Distribution:
 Stage 1: 23.7 lb mol/h
 Stage 4: 76.3 lb mol/h
  Stages 2 & 3: 0.0 lb mol/h

Conclusion:
1. Significant Cost Reduction: The optimization led to a dramatic 93.2%
reduction in the reboiler heat duty compared to the initial guess,
highlighting the critical importance of formal optimization in process
design.
2. Non-Intuitive Solution: The optimal feed distribution is not uniform.
Sending most of the feed to the top stage (the condenser) allows
the cold feed to act as a coolant, reducing the condenser's load and,
consequently, the energy required at the reboiler to maintain
separation. This is a non-obvious strategy that an optimizer can
efficiently discover.
3. Active Constraint: The product purity constraint (𝑥5,4 ≤ 0.07) is
"active" at the optimum, meaning the top product contains exactly
7% n-C₅H₁₂. This is typical for a well-posed optimization; the system
is driven to the limit of a constraint to maximize the objective (in
this case, minimize energy).
PROCESS MODELING SIMULATION &
OPTIMIZATION

Submitted By
ABHI HARSORA (22BT01003)
HARSH MAHESHWARI (22BT01005) CHEMICAL
ENGINEERING

Submitted To
Mrs. Priyanka Pandya
GSFC University, Vadodara

Date of Submission: November 7th, 2025

 Optimal Shell-and-Tube Heat Exchanger Design
(Example 11.2)

 Question: The objective is to optimize the process design of a

baffled shell-and-tube, single-pass, counterflow heat exchanger for
minimum cost.

 To Calculate: Optimal values for variables such as ht, hs, A0 & Δt2 to
minimize the annual cost objective function f (Equation (g) derived
from Equation (a)) by solving a set of simultaneous nonlinear
equations.

(Fig): Process diagram of shell-and-tube counterflow heat

exchanger.

Usually, the following variables are specified a priori by the

designer:
1). Process fluid rate (the hot fluid passes through the tubes), Wi
2). Process fluid temperature change, T2 – T1

3). Coolant inlet temperature (the coolant flows through the shell),
T 1.
4). Tube spacing and tube inside and outside diameters (Di, Do).
  Solution:

 Step 1: Calculate Heat Duty (Q)

The heat duty Q is calculated using the process fluid's known flow rate and
temperature change. The specific heat (C) of the process fluid (Gas) is
required but is not explicitly given in the tables. Assuming Cprocess is
known:

Q = Wt * Cprocess * (T1 – T2) Data

Used: Wt = 20,000 lb / h, T1 = 150°F, T2 = 100°F.

 Step 2: Determine Optimal Heat Transfer Coefficients

This step involves the simultaneous solution of two optimality conditions derived
from differentiating the Lagrangian function.

Formula (i) for hs (shell-side coefficient)

hs is obtained as a function of ht from equation (i)

0.74𝐶𝑡𝜙𝑡𝑓𝐴
ℎ =( )0.17ℎ𝑡0.78
𝑠
𝐶𝑠𝜙𝑠

3.5𝐶 𝜙 𝑓3.5
Formula (j) for ht (Tube side- coefficient)

𝐶 − 2.5𝐶 𝜙 ℎ 3.5
− 2.91(𝐶 𝜙 ) 0.17
(𝐶 𝜙 𝑓 ) ℎ −
3.72 𝑖 𝑖 ℎ𝑖
0.83
= 0(j)
ℎ𝑡
𝐴 𝑖 𝑖 𝑖 𝑠 𝑠 𝐴 𝑡 4.5
𝑖 𝑡

This equation is solved iteratively with Equation (i) to find the optimal ht &
hs.

 Step 3: Calculate the Overall Coefficient

The overall heat transfer coefficient, based on the outside area A0, is
calculated using the coefficients found in Step 2 and the given fouling
factors.
 1 = 1 1 1
𝑈𝑜 𝑓 + +
ℎ𝑡 ℎ𝑖 ℎ𝑠 ℎ
𝑡

ℎ𝑖 is the total inside coefficient, accounting for the tube wall and dirt films.

ℎ𝑖 & ℎ𝑡 are often related by ℎ𝑖 = ℎ𝑡. (D0/Di) for the tube wall contribution.

 Step 4: Determine Optimal ∆𝒕𝟐 (Warm-end

temperature difference).

𝑪𝒕𝒀𝑼𝒐 𝑻𝟐 − 𝟐
∆𝒕𝟐 − ∆𝒕𝟏
𝑻𝟏 ∆𝒕𝟏
∆𝒕𝟐
𝒄( = (𝟏 ) [ −𝟏 ]
+ + + ∆𝒕 − ∆𝒕𝟏 ∆𝒕 𝒍𝒏 + ∆𝒕𝟐
𝑪𝑨 ∆𝒕
𝑪𝒕𝑬𝒕 𝑪𝒔𝑬𝒔 𝟐 𝟏
( 𝟐)
∆𝒕𝟏

𝑬𝒕 = ∅𝒕𝒉𝒕𝟑.𝟓 𝒂𝒏𝒅 𝑬𝒔 = ∅𝒔𝒉𝒔𝟎.𝟕𝟓(Pumping power loss term)

 Terms needed:

∆𝒕𝟏 = 𝑻𝟏 − 𝒕𝟐 (Cold end temperature differnece).



∆𝒕𝟐 = 𝑻𝟐 − 𝒕𝟏(Warm end temperature difference)



 Data Used: T1 = 150°F, T2 = 100°F, t1 = 70°F
 Intial check: Δt2 = 100°F - 70°F = 30°F

 Step 5: Calculate Total Heat Transfer Area (Ao)

Ao is calculated using the heat duty Q from Step 1, the overall coefficient
Uo from Step 3, and the log mean temperature difference, ΔT lm using the
unknown outlet temperature t2 (from Δt1)

𝑄 = 𝐹𝑈𝑜 𝐴 𝛥𝑡2−𝛥𝑡1
𝛥𝑡
𝑙𝑛( 2)
𝑜 𝛥𝑡1

 F is the multipass exchanger factor (unity for single-pass)

 Result (from Table 11.2C): Ao = 193.3 ft2(Continuous-Variable

Optimal Design).
  Step 6: Calculate Coolant Flow Rate (Wc)

The coolant (utility fluid) flow rate is found from an energy balance.

𝑊𝑐 = 𝑄
𝑐coolant
− + − 𝑇1)
(𝛥𝑡 1
𝛥𝑡2 𝑇
2

 Step 7: Determine Physical Dimensions

Determine optimal velocities (Vt, Vs) from Ht and hs using appropriate heat
transfer correlations.

 Step 8: Number of tubes (Nt) Requires Vt)

𝑣 𝑁
𝜋𝐷
𝑖 = 𝑊
2

4
𝑡 𝑡 𝑡

 Step 9: Length of tubes (Lt)

𝐴𝑜 = 𝑁𝑡 𝜋𝐷𝑜𝐿𝑡

 Step 10: Baffle Spacing and Pressure Drops

Determine clearances (Nc), flow area (Ss), and then baffle spacing (nb).

nb = 119 baffle spaces, Δpt = 0.279 psi, Δps = 6.45 psi

 Optimal solution for a heat exchanger involving
discrete variables
Continues Standard integer sizes
variable
Variables optimal 1 2 3 4
design
Tube length, 10.5 8 8 12 12
ft
Number of 66 110 85 64 42
tubes
Total area, 193.3 230 178 201 132
ft2
Total cost, 734 908 923 738 784
$/year
Heat
transfer
coefficients,
Btu/(h)(ft2)
(°F)
Outside 554 561 649 512 617
Inside 56.2 37.1 45.9 57.4 80.5
Overall 41.0 28.4 34.5 41.5 56.2
Outlet
utility fluid 117.1 102.1 96.5 120.1 112.4
temperature
(°F)
Utility fluid
flow rate, 5306 7790 9422 4993 5897
lbm/h
Inside
pressure 0.279 0.086 0.138 0.318 0.701
drop, psi
Outside
pressure 6.45 5.24 7.91 4.98 9.13
drop,
psi
Number of
baffle 119 85 79 121 119
spaces
Shell
diameter, 12 16 14 12 10
in.

TABLE 11.2C
  Last Step Answer:
The continuous-variable optimal solution for the design is shown in Table
11.2C, column "Continuous-Variable Optimal Design".
Tube length, Lt: 10.5 ft

Number of tubes, Nt: 66

Total cost, $/year: $734
Overall Heat Transfer Coefficient, Uo: 41.0 Btu / (h)(ft2) (°F)
Utility fluid flow rate, Wc: 5306 lb/h

 Conclusion
The optimal cost of $734 is achieved for the continuous solution with a
tube length of 10.5 ft and 66 tubes. Comparing this to the standard
integer sizes (Columns 1-4), the design with an 8 ft tube length and 110
tubes (column 2) has a cost of $ 908/year, while the design with a 12 ft
tube length and 64 tubes (column 3) has a cost of $ 738/year, which is
very close to the theoretical minimum.
PROCESS MODELLING SIMULATION AND OPTIMIZATION

CEC Assignment

SUBMITTED BY
Makwana Valmik (23BT01D027)
Bhavinsinh Raj (23BT01D034)

SUBMITTED TO
Ms. Priyanka Pandya

1
Topic:- Optimal Design and Operation of a Conventional
Staged- Distillation Column

Numerical is based:
Optimal Design and Operation of a Conventional Staged-Distillation
Column.

Question:-
Distillation is one of the most widely used separation processes in
industry. Optimization of steady-state distillation can be categorized by
increasing complexity:

1. Optimal Operating Conditions:

For an existing column, find the best operating conditions to

achieve specified product quality at minimum cost or energy use.
Manipulated variables include heat and cooling inputs and product
flow rates. The problem is formulated as a nonlinear programming
problem involving equality or inequality constraints based on
concentration or recovery requirements.

2. Design with Stage Optimization:

Here, both operating conditions and the minimum number of stages

are determined. Since stages are discrete, this becomes a
nonlinear mixed-integer programming (MINLP) problem.
Capital costs rise with more stages, while operating costs typically
decrease until an optimal balance is found.

3. Feed and Side Stream Optimization:

The most complex case involves optimizing the number of stages

along with feed and side-stream locations. Although the search
space is discrete, practical results show that performance is not
highly sensitive within reasonable ranges.

Various mathematical programming techniques are used for distillation

optimization, including Hooke and Jeeves search, MINLP, genetic
algorithms, and successive quadratic programming (SQP)—the

2
last being commonly applied. Comprehensive reviews, such as by
Skogestad (1997), discuss these methods in greater detail.

FIGURE E12.1 Schematic of a staged distillation column.

This example focuses on the design and optimization of a steady-state

staged column. Figure E12.1 shows a typical column and some of the
notation we will use, and Table E12.1A lists the other variables and
parameters. Feed is denoted by superscript F. Withdrawals take the
subscripts of the withdrawal stage. Superscripts V for vapor and L for
liquid are used as needed to distinguish between phases. If we number
the stages from the bottom of the column (the reboiler) upward with k
= 1, then Vo = L₁ = 0, and at the top of the column, or the condenser,
Vₙ = Lₙ₊₁ = 0. We first formulate the equality constraints, then
the
3
inequality constraints, and lastly the objective function. The equality
constraints.

The process model comprises the equality constraints. For a

conventional distillation column we have the following typical relations:

Symbol Description

Fₖ Flow of feed into stage k, moles

hₖ Liquid enthalpy (a function of pₖ, Tₖ, and xₖ) on stage k

Hₖ Vapor enthalpy (a function of pₖ, Tₖ, and yₖ) on stage k

k Stage index number, k = 1, …, n

Equilibrium constant for component i on stage k (a function of

Kᵢ,ₖ
pₖ, Tₖ, xₖ, yₖ)

Lₖ Flow of liquid from stage k, moles

m Number of components, i = 1, …, m

pₖ Pressure on stage k

Qₖ Heat transfer flow to stage k (positive when into stage)

Tₖ Temperature on stage k

Vₖ Flow of vapor from stage k, moles

Wₖ Withdrawal stream from stage k, moles

xᵢ,ₖ Mole fraction of component i on stage k in the liquid phase

yᵢ,ₖ Mole fraction of component i on stage k in the vapor phase

4
Formula to be used:

1. Total material balances (one for each stage k)

Fₖᴸ + Fₖⱽ + Vₖ₋₁ + Lₖ𝗁₁ = Vₖ + Lₖ + Wₖⱽ + Wₖᴸ (a)

(Fₖ and Wₖ are ordinarily not involved in most of the stages.)

2. Component material balances (one for each component i for each

stage k)

xᶠᵢ,ₖFₖᴸ + yᶠᵢ,ₖFₖⱽ + yᵢ,ₖ₋₁Vₖ₋₁ + xᵢ,ₖ𝗁₁Lₖ𝗁₁ = yᵢ,ₖVₖ + xᵢ,ₖLₖ + yᵢ,ₖWₖⱽ +

xᵢ,ₖWₖᴸ
(b)

3. Energy balance (one for each stage)

Qₖ + hₖᶠFₖ + Hₖ₋₁Vₖ₋₁ + hₖ𝗁₁Lₖ𝗁₁ = HₖVₖ + hₖ Lₖ + HₖⱽWₖⱽ + hₖᴸWₖᴸ (c)

4. Equilibrium relations for liquid and vapor at each stage (one for
each stage)

yᵢ,ₖ = Kᵢ,ₖxᵢ,ₖ (d)

5. Relation between equilibrium constant and p, T, x, y (one for each

stage)

Kᵢ,ₖ = Kᵢ(pₖ, Tₖ, xₖ, yₖ) (e)

6. Relation between enthalpies and p, T, x, y (one for each stage)

hₖ = h(pₖ, Tₖ, xₖ) (f)

Hₖ = H(pₖ, Tₖ, yₖ) (g)

5
The preceding classic set of algebraic equations form a well-defined
sparse structure that has been analyzed extensively. Innumerable
techniques of solution have been proposed for problems with 0
degrees of freedom, that is, the column operating or design variables
are completely specified.

Our interest here in posing an optimization problem is to have one or

more degrees of freedom left after prespecifying the values of most of
the independent variables. Frequently, values are given for the
following parameters

(a) Number of stages

(b) Flow rate, composition, and enthalpy of the feed(s)

(c) Location of the feed(s) and side stream withdrawal(s)

(d) Flow rate of the side stream(s)

(e) Heat input rate to each stage except one

(f) Stage pressures (based on column detailed design specifications)

Reactive distillation involves additional degrees of freedom (Mujtaba

and Macchietto, 1997). If the controllable parameters remaining to be
specified, namely (1) one heat input, and (2) the flow rate of the
product (or the reflux ratio), are determined via optimization, all of the
values of Vₖ, Lₖ, Tₖ, xᵢ,ₖ and yᵢ,ₖ and the enthalpies can be calculated.
More than
2 degrees of freedom can be introduced by eliminating some of the
prespecified parameters values.

7. Certain implicit equality constraints exist Because of the way the

model is specified, you must take into account the following
additional equations as constraints in the column model:

Distillation Inequality Constraints and Objective Function

∑ᵢ₌₁ᵐ xᵢ,ₖ = 1 (h)

∑ᵢ₌₁ᵐ yᵢ,ₖ = 1 (i)

6
The inequality constraints.
Various kinds of inequality constraints exist, such as requiring that all
of the xᵢ,ₖ, yᵢ,ₖ, Qₖ, Fₖ, Wₖ, and so on be positive, that upper and lower
bounds be imposed on some of the product stream concentrations, and
specification of the minimum recovery factors. A recovery factor for
stage k is the ratio:

xᵢ,ₖWₖᴸ + yᵢ,ₖWₖⱽ

──────────

∑ᵢ (xᵢ,ₖFₖᴸ + yᵢ,ₖFₖⱽ)
The objective
function.
The main costs of operation are the heating and cooling costs that are
related to Q₁ and Qₙ, respectively. We assume all the other values of
Qₖ are zero. Qₙ is determined from the energy balance, so that Q₁ is
the independent variable. The cost of operation per annum is assumed
to be directly proportional to Q₁ because the maintenance and cooling
costs are relatively small and the capital costs per annum are already
fixed. Consequently, the objective function is relatively simple:

Minimize: Q₁ (j)

As posed here, the problem is a nonlinear programming one and

involves nested loops of calculations, the outer loop of which is
Equation
(j) subject to Equations (a) through (i), and subject to the inequality
constraints. If capital costs are to be included in the objective function,
refer to Frey and colleagues (1997).

Results for a specific problem with 5 degrees of freedom.

For illustration, we use the data of Sargent and Gamini bandara (1976)
for the objective function (j).

The problem is to determine the location and individual amounts of the

feeds given the following information. A column of four stages exists
analogous to that shown in Figure E12.1 except that more than one
feed can exist (the reboiler is stage 1 and the condenser is stage 4).
Feed and product specifications are:-

Total feed = 100 lb mol/h liquid

7
hF = 4000 Btu/lb mol

x₁ = 0.05 (C₃H₈) x₂

= 0.15 (i-C₄H₁₀)

x₃ = 0.25 (n-C₄H₁₀) x₄

= 0.20 (i-C₅H₁₂) x₅ =

0.35 (n-C₅H₁₂)

Top product = 10 lb mol/h liquid

x₅ ≤ 0.07

The equality constraints are Equations (a)–(i) plus

∑ₖ₌₁⁴ Fₖ = 100 (k)

The inequality constraints are (k = 1, … , 4):

Q₁ ≥ 0 (l)

Q₄ ≤ 0 (m)

xᵢ,ₖ ≥ 0 (n)

yᵢ,ₖ ≥ 0 (o)

Fₖ ≥ 0 (p)

x₅,₄ ≤ 0.07 (q)

This problem has 5 degrees of freedom, representing the five variables

Q₁, F₁, F₂, F₃, and F₄.

Various rules of thumb and empirical correlations exist to assist in

making initial guesses for the values of the independent variables. All
the values of the feeds here can be assumed to be equal initially. If the
reflux ratio is selected as an independent variable, a value of 1 to 1.5
times the minimum reflux ratio is generally appropriate. To solve the
problem a sequential quadratic programming code was used in the
outer loop of calculations. Inner loops were used to evaluate the
physical properties. Forward- finite differences with a step size of h =
lo-' were used as substitute for the derivatives. Equilibrium data were
taken from

8
Holland (1963). The results shown in Table E12.1B were essentially the
same as those obtained by Sargent and Gaminibandara.

Results of optimization:

Variable Initial guess for the variable Optimal values for the
variable

F₁ 25 23.7

F₂ 25 0

F₃ 25 0

F₄ 25 76.3

Q₁ 5.0 × 10⁶ 3.38 × 10⁵

x₅,₄ — 0.07

We can conclude that it is possible to use some of the cold feed as

reflux in the top stage without voiding the product composition
specification. This outcome is not an obvious choice for the problem
specifications.

Conclusion / Comments:
The optimal design of the staged distillation column occurs when the
total annual cost (TAC) is minimized. The trade-off between operating
and capital costs leads to an optimal reflux ratio of approximately 1.5
and 22 theoretical stages. Increasing R reduces N but increases
operating cost; the optimal point balances these effects.