Process Modeling methods and tools Lecture 9 Ville Alopaeus

1
Process Modeling methods and tools Lecture
9Ville Alopaeus
2
Outline

• Need for partial differential equations
• Some numerical methods
• Finite differences
• Polynomial approximations
• Finite volume methods
• Spectral methods
• Method order and other numeric issues

3
Partial differential equations

• For a function
• changing x leads to changes in f as

f is a function of one variable only. These
result in ordinary differential equation models
4
Partial differential equations

• For a function
• changes in variables lead to changes in f as

by the chain rule. PDEs result in.
5
PDE

• More generally, independent variables x and y,
  dependent c,

Ordinary differential equation
Partial differential equation, differentiation of
the field with respect to several variables
6

• We live in three spatial and one temporal
  dimensions. Therefore models describing real life
  should be partial differential equations
• Usually partial differential equations are more
  complex to be solved than ordinary differential
  equations

7

• If none of the dimensions affect our modeled
  property, we can reduce our model into an
  algebraic model equation (in case of one
  dependent variable)
• If one dimension affects the modeled property but
  the others do not, we end up with an ordinary
  differential equation

8

• If more than one dimension affect our modeled
  property, then we generally end up with a partial
  differential equation

if
and x and y are comparable, the last term can be
cancelled -gt ODE
9
Example tubular reactor
radial profiles (concentration, temperature etc.)
Steady state concentration at any location
within the reactor and any time depend only on
length coordinate ? ordinary differential
equation model (may be several dependent
variables, e.g. components)
10
Example tubular reactor
r
z
radial profiles (concentration, temperature etc.)
Steady state, but concentration at any location
within the reactor and any time depend on axial
and radial locations
PDE
11
Dynamic plug flow
With dimensionless length coordinate
If ?res small compared to changes in c (or
boundary conditions), and since rr(z,t)
Then
ODE
12
A warning

• In many cases, only one spatial dimension is
  important but be careful

A long and narrow cable
A long and narrow wood chip or fibre
Diffusion flux
Heat flux
Diffusion flux
Diffusion along the fibre length cannot be
neglected due to anisotropy!
End effects can be neglected
13
Heat transfer example
qz
y
qx
qx?x
x
z
qy
Heat balance for the cube
14
Heat transfer example
Fouriers law
in each direction
15
Heat transfer
??? differential heat balance for the cube
Insert Fourier law with constant properties
16
Heat transfer
expanded in spatial directions
Second order parabolic partial differential
equation.
17
One-dimensional transient heat transfer
Heat transfer in a slab
Other surfaces insulated or symmetry (variations
in other directions negligible)
Cold surface, Temperature specified
Hot surface, Temperature specified
Transient heat conduction in x-direction only
18
Some solution methods

• Usually partial differential equations are
  discretized in some way to transform them to a
  set of ordinary differential equations
• Similar methods can be used to solve ordinary
  differential equations (transform them into a set
  of algebraic equations)
• There are numerous methods for this.

19
Finite differences

• Perhaps the most easily comprehendable method
• Derivatives are replaced by differences

20
Discretized one-dimensional transient heat
transfer
T1
T2
Transient heat conduction in x-direction only
21
PDE in weather forecast

• Lewis Fry Richardson (1881-1953)
• Proposed in 1922 that the differential equations
  describing atmosphere could be solved numerically
  (!) by using finite differences
• Before that, case-based reasoning was used

22
Richardson weather forecast factory
23
Richardson weather forecast factory

• Richardson estimated that 64 000 people would be
  needed to keep pace with the atmosphere
• Currently computer capacity in European weather
  forecast is 1010 times faster than the
  Richardson factory, corresponding to 1015 people
  doing hand calculations
• Also numerical methods have evolved considerably

24
The same Richardson...

• Extrapolation from two grid sizes

Method order has to be known
25
Method order

• Error is proportional to the grid size to the
  power of the method order

Alternative notations
For first order methods, halving the grid size
halves the error Second order methods ½ grid
size, ¼ error etc...
26
Method order

• Usually low order methods are inaccurate
  (numerical diffusion)
• High order method are often less stable or are
  susceptible to oscillations
• High order does not necessarily mean more
  accurate
• Low order almost necessarily mean less accurate

27
Method order

• Method order can be estimated from the remaining
  term in a Taylor expansion

Error is proportional to the first term that has
been left out eO(hn)
28
Apparent method order

• Calculate a solution with two grid sizes. Compare
  it to the exact solution (analytical or with a
  very dense grid).

Grid 1
Grid 2
29
Apparent method order

• Calculate a solution with two grid sizes. Compare
  it to the exact solution (analytical or with a
  very dense grid).

30
High order approximate solutions to PDEs

• Usually the method order is approximately
  proportional to the number of connections between
  independent variables in a numerical scheme

Two-point difference
eO(h)
or
31
High order approximate solutions to PDEs
Three-point difference
eO(h2)
Where is the third point?
32
High order approximate solutions to PDEs
Taylor series expansions
33
High order approximate solutions to PDEs
Subtract second from the first
Divide by 2h
34
High order approximate solutions to PDEs
Approximation for the derivative
Approximation for the error eO(h2)
The third point is actually fi, but it has a
coefficient 0 on the right hand side
35
High order approximate solutions to PDEs
Five-point difference
eO(h4)
etc. Note that discretization of the boundary
conditions affect the method order. Inaccurate
treatment of the boundaries may ruin the whole
solution accuracy!
36
Advection equation

• Concentration is advected with a flowing fluid
• Hyperbolic partial differential equation
• If v constant ? inital concentration distribution
  just moves without shape change

37
Advection equation

• Discretized with interval 0,1, v1 m/s, t0.5 s

Concentration initially everywhere zero but on
the interval 0.1,0.2
38
Advection equation

• Right hand side is discretized with various
  methods.
• Three point central difference

39

• Three point central difference
• 100 discretization points

40

• First order upwind discretization
• 100 discretization points

41

• First order upwind discretization
• 1000 discretization points

42

• Second order upwind discretization
• 100 discretization points

43
CFL

• In principle, increasing the number of grid
  points we get more and more accurate solution
• This increases number of variables to be solved
  at each time step
• According to the Courant-Friedrichs-Lewy
  stability (or convergence) criterion, maximum
  step size is limited by the grid size

44
CFL

• CFL condition states that time step should be
  such that the flow does not cross several
  discretization points during one time step
• Courant number
• Maximum Courant number value depends on the
  equation to be solved, but usually 1

45
Polynomial approximations

• Usually polynomial approximations are constructed
  so that a trial polynomial
• is formulted in such a manner (its coefficients
  are calculated so) that it satisfies the
  differential equation in best possible way

46
Polynomial approximations

• Tis are trial functions. Here domain always
  scaled to 0,1
• T0 satisfies the boundary conditions of the given
  differential equaiton.
• Each trial function Ti satisfies homogeneous
  boundary conditions

47
Polynomial approximations

• Here only one spatial dimension is considered
  with time-dependent solution
• Only spatial dimension is discretized.
• This method can be used to bring a ODE into a set
  of algebraic equations or a PDE into a set of
  initial value ODEs

48
Polynomial approximations

• There are several criteria how to define the
  optimal solution. For a differential equation

we define a residual function as
49
Various polynomial approximations
1. Collocation method
The trial function satisfies the differential
equation at specific points
50
Various polynomial approximations
1. Collocation method
It can be shown that the optimal collocation is
based on zeros of Jacobi polynomials (the
differential equation is satisfied at those
points) This is called Orthogonal collocation.
51
Various polynomial approximations
Orthogonal collocation
Jacobi polynomials are defined by an
orthogonality condition
52
Various polynomial approximations
Zeros of Jacobi polynomials
E.g. for ?,? 0
These are also used in numerical integration
(Gaussian quadrature)
53
Various polynomial approximations
2. Subdomain method
Integral over each subdomain is set to zero
54
Various polynomial approximations
3. Method of moments
55
Various polynomial approximations
4. Galerkin method
The Finite Element Method is an application of
the same principle
56
Various polynomial approximations
5. Least Squares method
An alternative formulation
57
OCFE

• Orthogonal collocation (and other methods as
  well) can be applied piece-wise. This increases
  accuracy and/or decreases oscillations.

58
Finite volume methods

• The basic idea is to formulate balances for a
  finite volume instead of differential volume
  element
• Also called as integral form of balances
• Often used in Computational Fluid Dynamics (CFD)

59
FVM

• Balance for a finite volume, fluxes calculated at
  the boundaries

qz
qx
qx?x
qy
60
FVM

• Control volume

Control volume boundary
61
FVM
Flux at the boundary must be calculated from the
cell values

• Central difference,
• Upwind (1st order, 2nd order...)

62
FVM

• Higher than first order methods tend to oscillate
  (Godunovs theorem)
• One possible solution is to use high order in
  smooth regions and limit fluxes at the boundaries
  in cases of steep gradients (this prevents
  oscillations)
• These problems appear mainly as related to the
  convection problem (hyperbolic part)

63
Some flux limiters
64
High order FVM
Fluxes at the interfaces are calculated from high
order polynomials. Possibly weighting from
several polynomials. ENO (essentially
non-oscillating) WENO (weighted ENO) etc
65
Spectral methods
Basis functions are waves for which amplitudes
(coefficients) are calculated. From moments
66
Spectral methods
Global weather forecast Each predicted field is
expanded in series of spherical harmonics, e.g.
Trial functions are functions of spatial location
only, and coefficients are functions of time only
67
Moment-based dynamic plug flow reactor model

• A plug flow reactor is modeled with polynomial
  approximations
• Property distributions (concentration profiles
  along reactor length) are modeled based on
  moments
• Source terms for moment time derivatives are
  obtained by solving the reactor model.

68
Moment-based dynamic plug flow reactor model
Reactor model

• Polynomials along reactor length (or part of it)

69
Moment-based dynamic plug flow reactor model
Chebyshev polynomials
70
Moment-based dynamic plug flow reactor model
Reactor model
Boundary conditions can also be solved
approximately, i.e. it is not necessary to
formulate the problem so that one basis function
satisfies the b.c. and others satisfy homogeneous
b.c.
71
Moment-based dynamic plug flow reactor model
Integrating by parts gives boundary conditions
First order approximation (linear profiles) is
actually equivalent to a series of stirred tanks
-model
72
Example Cromatographic separation

• Polynomial profile model is applied to a
  cromatographic separation
• At one end of a separation column, a short pulse
  of two components is introduced
• The two components have different adsorption
  coefficients to the column material

73
Constant concentrations (CSTRs in series), 40
elements
40 variables
74
Linear profile within each element, 20 elements
40 variables
75
3. degree polynomial in each element, 10 elements
40 variables
76
Conclusions

• Partial differential equations are needed in
  cases where several dimensions affect the desired
  results
• Finite differences are perhaps the first and
  simplest approach for PDEs
• High order methods converge rapidly to the exact
  solution as the number of variables is increased,
  but may predict oscillatory solutions

77
Conclusions

• Various polynomial approximations are developed
  for boundary value ODEs and PDEs
• They are based on setting a residual equation to
  zero (various ways ? various methods)
• Finite volume methods are formulated directly for
  integral balances