Numerical simulation of chemical

engineering system

Prepared by:
Bikash Kumar Mahato
Assistant Professor
IOE, Pulchowk Campus
 Overview of Numerical Method
• Differential equations play a dominant role in mathematical modeling.
In practical engineering applications, only a very limited number of
them can be solved analytically.
• The purpose of this unit is to give an introduction to the numerical
methods needed to solve differential equations, and to explain how
solution accuracy can be controlled and how stability can be ensured
by selecting the appropriate methods.
• The mathematical framework needed to solve both ordinary and
partial differential equations is presented.
 ODEs
• A characteristic of a differential equation is that it involves an unknown
function and one or more of the function’s derivatives. If the unknown
function depends on only one independent variable, it is classified as an
ordinary differential equation (ODE).
• The order of the differential equation is simply the order of the highest
derivative that appears in the equations.
• Consequently, a first-order ODE contains only first derivatives, while a
second-order ODE may contain both second and first derivatives.
• The ODEs can also be classified as linear or non-linear.
• An ODE has an infinite number of solutions, but with the appropriate
conditions that describe systems, i.e. the initial value or the boundary
value, the solutions can be determined uniquely.
 Classification of ODEs
• Ordinary Differential Equations (ODEs) can also be classified based on the conditions
provided for their solutions. The two primary types of problems in this context are
Initial Value Problems (IVPs) and Boundary Value Problems (BVPs).
• Initial Value Problems (IVPs)
• An Initial Value Problem specifies the value of the unknown function and possibly its
derivatives at a single point. This type of problem is often used to describe dynamic
systems where the state of the system at a particular time is known, and the goal is to
determine the state of the system at future times.
• General statement of the problem
• An IVP for a first-order ODE has the form: dy/dx = f(x,y), y(x0) = y0
• where y(x0) = y0 is the initial condition.
• IVPs are often solved using methods like Euler's method, Runge-Kutta methods, or
other numerical integration techniques.
 Conti.
• Boundary Value Problems (BVPs)
• A differential equation that has data given at more than one value of the independent
variable is a boundary-value problem (BVP). Consequently, the differential equation
must be of at least second order.
• This type of problem is commonly used in situations where the state of the system is
known at the boundaries of the domain, and the goal is to determine the behavior of the
system within those boundaries.
• BVP for a second-order ODE has the form:

• BVPs are solved using methods like the finite difference method, shooting method,
or finite element method.
 PDEs
• Many phenomena in chemical engineering depend, in complex ways, on space and
time, and often a mathematical model requires more than one independent variable to
characterize the state of a system, i.e. the systems need to be described using partial
differential equations, PDEs.
• Examples of such phenomena include chemical reactions, heat transfer, fluid flow etc. For practical
engineering applications, analytic solutions do not exist, and numerical methods need to be applied.
• This section is not intended to give a complete discussion of PDEs nor of solution
methods. Instead, the aim is to introduce the terminology and some issues involved in
solving PDEs.
• The discussion will be limited to linear PDEs that have two independent variables, e.g.
space and time, or two space variables for a steady-state problem, with the form
 Conti.

• where the subscripts denote the partial derivatives of the dependent

variable, u, and the focus is on the finite difference method.
 Classification of PDEs
• PDEs can be classified in different ways. The classification is important because the
solution methods often apply only to a specified class of PDEs.
• To start with, PDEs can be classified by the number of variables, e.g. ut = uxx, which
contain two independent variables, t and x. The order of a PDE is the order of the
highest order derivative that appears in the PDE.
• For example, ut = ux is a first-order PDE, whereas ut = uxx is a second-order PDE.
• In addition, it is important to make a distinction between nonlinear and linear PDEs.
• An example of a well-known non-linear PDE is the Navier–Stokes equation, which
describes the motion of fluids.
• In a linear PDE, the dependent variable and its derivatives appear in a linear fashion.
The linear second-order PDE in above Equation can be classified as
• parabolic, if B2–4AC = 0;
• hyperbolic, if B2–4AC > 0;
• elliptic, if B2–4AC < 0.
 Conti.
• The B2–4AC term is referred to as the discriminant of the solution, and the
behavior of the solution of Equation depends on its sign. The smoothness of a
solution is affected by the type of equation.
• The solutions to elliptic PDEs are smooth. In addition, boundary conditions at
any point affect the solution at all points in the computational domain.
• Elliptic PDEs often describe the steady-state condition of a variable,
• for example a diffusion process that has reached equilibrium, or a steady-state
temperature distribution.
• In hyperbolic PDEs, the smoothness of a solution depends on the initial and
boundary conditions.
• For instance, if there is a jump in the data at the start or at the boundaries, this jump will
propagate in the solution.
• Parabolic PDEs arise in time-dependent diffusion problems,
• for example the transient flow of heat.
 Conti.
• PDEs can be solved using following methods:
• Finite difference solution of parabolic equations
• Forward difference method
• Backward difference method
• Crank–Nicolson method
 Sensitivity analysis concepts.
• Sensitivity analysis is a method used to determine how the output of a model or system
is affected by changes in its input parameters. This concept is essential in various fields
such as engineering, economics, finance, and environmental science, as it helps to
identify which variables have the most significant impact on the model's outcomes and
to understand the robustness and reliability of the model.
• Key Concepts of Sensitivity Analysis
• Purpose:
• Identify Influential Parameters: Determine which input variables most influence the
output.
• Assess Model Robustness: Understand how uncertainties in inputs affect the
model’s predictions.
• Prioritize Data Collection: Identify which parameters need more accurate
measurement or further study.
• Support Decision Making: Provide insights for making informed decisions based
on the model.
 Conti.
• Types of Sensitivity Analysis:
• Local Sensitivity Analysis: Examines the effect of small changes in input parameters on the
output, typically using partial derivatives or one-at-a-time (OAT) approaches.
• Global Sensitivity Analysis: Investigates the effect of varying multiple input parameters
simultaneously across their entire range, often using statistical methods or variance-based
techniques.
• Methods:
• Derivative-Based Methods:
• Derivative-based sensitivity analysis (DBSA) is a local sensitivity analysis method that involves calculating
the partial derivatives of a model's output with respect to its input parameters.
• These partial derivatives indicate how sensitive the output is to small changes in each input parameter.
• Screening Methods:
• Ford and Flynn (2005) suggest Pearson correlation in order to conduct quick sensitivity analysis, called
screening.
• In this method, correlation coefficients between the output and each parameter are calculated and plotted
against simulation time
• Parameters that have high correlation with output variable are concluded to be the high sensitivity ones
 Conti.
• Variance-Based Methods:
• Variance-based sensitivity analysis is a global sensitivity analysis method that decomposes the
output variance of a model into contributions from each input parameter and their interactions.
• This approach helps to understand the impact of input uncertainty on the output variability and
identifies which parameters are most influential in driving the model's behavior..
• Monte Carlo Simulation:
• Monte Carlo simulation is a powerful and flexible method for sensitivity analysis, particularly
useful for complex models where analytical solutions are difficult to obtain.
• The basic idea is to use random sampling to explore the behavior of the model under different
input conditions and quantify how variations in inputs affect the outputs.
• Steps in Conducting Sensitivity Analysis:
• Define the Model: Establish the mathematical or computational model to be analyzed.
• Select Input Parameters: Identify which parameters to vary and their ranges.
• Perform Simulations: Run the model multiple times with different sets of input values.
• Analyze Results: Evaluate how changes in inputs affect the output, using statistical or graphical
methods.
• Interpret Findings: Determine which parameters are most significant and how sensitive the
model is to changes in these parameters.
 Conti.
• Applications of Sensitivity Analysis
• Engineering: Optimize design parameters and ensure system reliability under
different operating conditions.
• Economics and Finance: Assess how changes in economic variables (e.g., interest
rates, inflation) affect financial models or investment portfolios.
• Environmental Science: Evaluate the impact of uncertainties in environmental
parameters (e.g., emission rates, climate sensitivity) on ecological models.
• Healthcare: Understand the effect of different treatment protocols or patient
characteristics on health outcomes.
 Solution of Navier-Strokes equation using FDM
• In today’s class we will study a popular finite difference method that is
known as, marker-and-cell method and we will solve two dimensional
unsteady Navier-Stokes equations.
• Continuity equation:- ∂u/∂x + ∂v/∂y = 0
• X-momentum equation:-

• Y-momentum equation:-
 Conti.
• So, if you know the fluid properties, then rho and nu are known. Then there
will be unknown u, v and p. So, we have three equations, one is continuity
equation and two momentum equations.
• So, three equations and three unknowns. So, we will be able to solve these
unknowns. But, there is no obvious pressure equation available. So, we need to
derive this special equation from the continuity equation.
• Here we will use marker and cell method. So, this marker and cell technique
was first proposed by Harlow and Welch in 1965 as a method for solving
free-surface problems. The staggered grid is used; which we will discuss
below.
• The u, v, p system is known as primitive variable approach and we solve these
equations in primitive variable approach.
• The differencing scheme for the momentum equations is basically just the
forward-time and centered-space FTCS method. Both for convection terms and
the diffusion terms and the pressure in the u, v, p system should be solved from
a Poisson equation and that we will derive from the continuity equation.
 Conti.
• There are two types of grids based on
the locations of the flow variables,
collocated grid and staggered grid.
• What is collocated grid? In collocated
grid all the solution variables u, v, p and
any scalars like temperature and spaces
are solved at the same node point;
• whereas in staggered grid we solve
pressure and scalars like temperature
and spaces in a particular node but
solve the velocities in staggered way.
• You can see this in above figure also.
Collocated grid sometime it is also
referred as non-staggered grid.
 Conti.
• What is the problem in collocated grid? So, in collocated grid there will be
pressure and velocity decoupling. So, to avoid this pressure and velocity
decoupling, we use staggered grid.
• So, in staggered grid, we solve u velocity in the horizontal face and v velocity in
the vertical. So, you can see these velocities are solved in staggered way and
pressure we solve at center of the cell. So, using staggered grid you can avoid
the pressure and velocity decoupling.
• So, you can see in the figure, all the horizontal arrows wherever it is there. We
solved the u velocities, you can see all the places. But when we solve v velocity,
it is in the vertical arrow. So, in vertical arrow the v velocity is noted.
• Let us now discretize the pressure gradient terms in collocated grid. The point
under consideration i.e., cell center or our node point is P and the Left, Right,
Top and Bottom cell centers are labelled as East (E), West (W), North (N) and
South(S) respectively.
• On similar terms the Left, Right, Top and Bottom face centers are labelled as
East (e), West (w), North (n) and South(s) respectively.
 Conti.
• Let us say that the grid distance is delta x and delta y and is same as the distance
between these face centers.

• The pressure values at the faces are not known. So, we have to interpolate and we take
the average value.
 Conti.
• You can see when we are solving x momentum equation at point P pressures involved
are east and west similarly when we are solving y momentum equation pressures
involved are at north and south not at P.
• So, you can see that although we are solving the velocities at a particular cell we are not
incorporating the pressure at that point in the pressure gradient. In this way, pressure
and velocity may decouple.
• Now let us invoke the continuity equation at the same cell.
 Conti.
• So, if you see this equation carefully, we are taking nodes at east face,
west face north face and south face. But, essentially we are solving the
continuity equation for the control volume having P as the center and not
other neighbouring points as center.
• So, obviously this approach is not correct as no velocities appear from the
main control volume. There may be some oscillations in the problem away
from our node point which may be unnecessarily incorporated in the
solution.
• Thus this is the drawback of collocated grid structure. The problem is also
known as Checker-board pressure distribution or Checker-board
velocity distribution.
 Conti.
• So, in this staggered grid you can see the general notation is same as
discussed above. Writing the continuity equation for cell containing point
P and its neighbours:

• Now, you can see that, the velocities are already available at those points.
Similarly, discretising the pressure gradient terms we get:
 Conti.
• We are solving the pressure at the face center of the main cell. So,
obviously we do not need any interpolation as the pressures are
available at these points already.
• Hence there will be no pressure and velocity decoupling and
essentially you are using the nearest grid point velocities and it is
second order accurate as you are using central difference. So, this is
the advantage of using staggered grid.
 Advantages and Disadvantages of collocated and
staggered grid.
• Obviously, collocated grid is easy to code,
• as we are solving all the variables u, v, p, temperature or spaces at the same
node.
• So it’s data structure is easy because you are having the same index for the u,
v, w.
• But, when you are solving using staggered grid,
• velocities are in staggered way and keeping this data structure is very difficult
because you have separate cell for u, v and p and accordingly you need to keep
the indices.
• But, the advantage is that in staggered grid pressure and velocity decoupling
will not be there.
 Conti.
• In previous class, we discussed why we should use staggered grid over collocated
grid and we have already discussed that in MAC (Marker and Cell) algorithm we use
primitive variable approach and staggered grid.
• So, today, using this MAC algorithm we will discretize the unsteady two-dimensional
Navier-Stokes equations and we will write the discretized equation. Finally, we will
discuss about the algorithm to solve this governing equations.
 Conti.

• Putting all the discretized values back into the momentum equation we
get:
 Conti.
• In the left-hand side, in the temporal term, you have 1 unknown velocity at n+1 level
and all other velocities are known. But if you see the pressure term it is defined at n+1
level which is unknown.
• So, we have to now simplify it in some way. For that we will calculate some
provisional velocity or predicted velocity u tilde assuming the pressure from the
previous time level n.
 Conti.
• So, now you can see that equation (3) and (4), we can uP tilde and vP tilde if we know
the pressure correction, then from equation (5) and (6) we can calculate velocities at
n+1 time level.
• So, now we need to find pressure correction and there is no separate pressure
equation thus, we need to derive it from the continuity equation.
• Now we will satisfy the continuity equation in the main cell of pressure P and we will
derive the equation for pressure correction.
 Conti.

• So, that means provisional velocity must satisfy the continuity equation. So, initially
provisional velocity, it will not satisfy the continuity equation but at convergence, the
divergence of the provisional velocity must tend to zero.
• This equation is difficult to solve because you have a pentagonal matrix. So MAC algorithm
suggests that, you assume the corrections of neighboring terms as 0 i.e., pressure correction in
the neighboring cells equals 0.
Conti.