Solved with COMSOL Multiphysics 5.

Process Control Using a PID

Controller
Introduction
In the chemical process industry it is often important to control a specific process. PID
control (proportional-integral-derivative-control) is one way to achieve that, but it can
be difficult to optimize the parameters in the PID algorithm. This example illustrates
how you can implement a PID control algorithm to simulate a process control system
and to find the optimal PID parameters.
This application is a generic example but could resemble the environment in a
combustion chamber where the concentration at the ignition point is crucial. Two gas
streams with different oxygen concentrations are mixed in the combustion chamber.
The concentration is measured at the ignition point before complete mixing of the
streams is reached. The control algorithm alters the inlet velocity of the gas with the
lower oxygen content to achieve the desired total concentration at the ignition point.

Model Definition
The model geometry appears in Figure 1. At the upper inlet, a gas stream with high
oxygen content enters the reactor at a velocity of 10 mm/s, while a gas with a lower
oxygen level enters from the left. The oxygen concentration is measured at a

1 |

PROCESS CONTROL USING A PID CONTROLLER

Solved with COMSOL Multiphysics 5.1

measurement point, and the inlet velocity of the less concentrated stream is altered by
the PID control algorithm to achieve the desired concentration at that point.

Upper inlet

Measurement point
Controlled inlet
Outlet

Figure 1: Model geometry.

The model uses the Laminar Flow interface to describe the fluid flow and the
Transport of Diluted Species interface for the mass balance. The corresponding
equations read (assuming incompressible flow and absence of reactions)

T
u
[ ( u + ( u ) ) ] + u u + p = 0
t
u = 0

----c+ ( D c ) = u c
t
To formulate the boundary conditions for the mass-transport equation, begin by
assuming that you know the two inlet concentrations. In addition, assume that the
reactant transport at the outlet is mainly driven by convection, that is, neglect diffusion
in the main direction of the convective flow. A no-flux boundary condition describes
all walls. The boundary conditions for the mass balance are:

2 |

BOUNDARY

CONSTRAINT

Upper inlet

c = cin,top

Controlled inlet

c = cin,inlet

PROCESS CONTROL USING A PID CONTROLLER

Solved with COMSOL Multiphysics 5.1

BOUNDARY

CONSTRAINT

Outlet

n (Dc) = 0

Walls

Nn=0

Here c is the concentration; cin,top and cin,inlet are the inlet concentrations (mol/m3)
for the upper and controlled inlets, respectively; D is the applied diffusivity (m2/s);
and N is the molar flux (mol/(m2s)).
The model uses the following boundary conditions for the fluid flow:
BOUNDARY

CONSTRAINT

Upper inlet

u = (0, vin,top)

Controlled inlet

u = (uin, 0)

Outlet

p0 = 0

Inlet sections

nu=0

Walls

u=0

Here u is the velocity vector (m/s), vin,top is the inlet velocity at the top inlet, and uin
is the PID controlled velocity. At the outlet, set the pressure to 0. No Slip boundary
conditions describe all walls except the inlet sections where slip conditions apply,
allowing for a smooth transition to a laminar velocity profile.
The PID control algorithm used to calculate uin is
t

u in = k P ( c c set ) + k I ( c c set ) dt + k D ----- ( c c set )

t

(1)

with the following parameters:

PARAMETER

VALUE

cset

0.5 mol/m3

kP

0.5 m4/(mols)

kI

1 m4/(mols2)

kD

10-3 m4/mol

In practice, the derivative constant, kD, is set to 0 in most cases as this parameter can
be difficult to determine. Moreover, the derivative term may increase the fluctuations
in the system because it amplifies noise in the error c cset.

3 |

PROCESS CONTROL USING A PID CONTROLLER

Solved with COMSOL Multiphysics 5.1

Results and Discussion

The two plots in Figure 2 show the oxygen concentration and the velocity stream lines
in the chamber after 0.05 s and 2 s, respectively. The figures show that the measured
concentration depends strongly on the flow field. At start-up, when the inlet velocity
of the stream entering from the left is very low, the sensor is entirely exposed to the
highly concentrated stream, and as the left inlet velocity increases the opposite relation
occurs.

4 |

PROCESS CONTROL USING A PID CONTROLLER

Solved with COMSOL Multiphysics 5.1

Figure 2: Oxygen concentration and velocity streamlines after 0.1 s (top) and 1.5 s
(bottom).
Figure 3 shows the inlet velocity and concentration in the measurement point as a
function of time for two different values for the kP parameter. The solid line represents
the results for a kP value of 0.5 m4/(mols) while the dashed line corresponds to kP

5 |

PROCESS CONTROL USING A PID CONTROLLER

Solved with COMSOL Multiphysics 5.1

equal to 0.1 m4/(mols). The results evaluated for the smaller kP value oscillate more
before stabilizing. Thus, it is clear that for this case the higher kP value yields a more
stable process control.

Figure 3: PID-controlled inlet velocity (top) and concentration in the measurement point
(bottom) as a function of time for kP = 0.5 m4/(mols) (blue) and kP = 0.1 m4/(mols)
(green).

6 |

PROCESS CONTROL USING A PID CONTROLLER

Solved with COMSOL Multiphysics 5.1

Application Library path: COMSOL_Multiphysics/Multiphysics/pid_control

Modeling Instructions
From the File menu, choose New.
NEW

1 In the New window, click Model Wizard.

MODEL WIZARD

1 In the Model Wizard window, click 2D.

2 In the Select physics tree, select Fluid Flow>Single-Phase Flow>Laminar Flow (spf).
3 Click Add.
4 In the Select physics tree, select Chemical Species Transport>Transport of Diluted
Species (tds).
5 Click Add.
6 In the Select physics tree, select Mathematics>ODE and DAE Interfaces>Global ODEs and
DAEs (ge).
7 Click Add.
8 Click Study.
9 In the Select study tree, select Preset Studies for Selected Physics Interfaces>Time
Dependent.
10 Click Done.
GLOBAL DEFINITIONS

Parameters
1 On the Home toolbar, click Parameters.
2 In the Settings window for Parameters, locate the Parameters section.

7 |

PROCESS CONTROL USING A PID CONTROLLER

Solved with COMSOL Multiphysics 5.1

3 In the table, enter the following settings:

Name

Expression

Value

Description

v_in_top

0.01[m/s]

0.01 m/s

Velocity, upper inlet

c_in_top

1[mol/m^3]

1 mol/m

Concentration, upper
inlet

c_in_inle
t

0.2[mol/m^3]

0.2 mol/m

Concentration,
controlled inlet

c00

0.5[mol/m^3]

0.5 mol/m

Initial concentration,
chamber interior

1e-4[m^2/s]

1E-4 m/s

Diffusivity

c_set

0.5[mol/m^3]

0.5 mol/m

Setpoint concentration

k_P_ctrl

0.5[m^4/
(mol*s)]

0.5 m4/(smol)

Proportional parameter

k_I_ctrl

1[m^4/
(mol*s^2)]

1 m4/(smol)

Integral parameter

k_D_ctrl

1e-3[m^4/mol]

0.001 m4/mol

Derivative parameter

GEOMETRY 1

Create the geometry. To simplify this step, insert a prepared geometry sequence.
1 On the Geometry toolbar, click Insert Sequence.
2 Browse to the applications Application Library folder and double-click the file
pid_control.mph.

3 On the Geometry toolbar, click Build All.

MATERIALS

Material 1 (mat1)
1 In the Model Builder window, under Component 1 (comp1) right-click Materials and

choose Blank Material.

2 Select Domain 1 only.
3 In the Settings window for Material, locate the Material Contents section.
4 In the table, enter the following settings:

8 |

Property

Name

Value

Unit

Property group

Density

rho

1.2[kg
/m^3]

kg/m

Basic

Dynamic viscosity

mu

3e-5

Pas

Basic

PROCESS CONTROL USING A PID CONTROLLER

Solved with COMSOL Multiphysics 5.1

LAMINAR FLOW (SPF)

Inlet 1
1 On the Physics toolbar, click Boundaries and choose Inlet.
2 Select Boundary 1 only.
3 In the Settings window for Inlet, locate the Velocity section.
4 In the U0 text field, type u_in_ctrl*(u_in_ctrl>0).

Inlet 2
1 On the Physics toolbar, click Boundaries and choose Inlet.
2 Select Boundary 7 only.
3 In the Settings window for Inlet, locate the Velocity section.
4 In the U0 text field, type v_in_top.

Outlet 1
1 On the Physics toolbar, click Boundaries and choose Outlet.
2 Select Boundary 13 only.

Wall 2
1 On the Physics toolbar, click Boundaries and choose Wall.
2 Select Boundaries 2, 3, 6, and 8 only.
3 In the Settings window for Wall, locate the Boundary Condition section.
4 From the Boundary condition list, choose Slip.
TR A N S P O R T O F D I L U T E D S P E C I E S ( T D S )

Transport Properties 1
1 In the Model Builder window, under Component 1 (comp1)>Transport of Diluted
Species (tds) click Transport Properties 1.
2 In the Settings window for Transport Properties, locate the Diffusion section.
3 In the Dc text field, type D.
4 Locate the Model Inputs section. From the u list, choose Velocity field (spf).

Initial Values 1
1 In the Model Builder window, under Component 1 (comp1)>Transport of Diluted
Species (tds) click Initial Values 1.
2 In the Settings window for Initial Values, locate the Initial Values section.
3 In the c text field, type c00.

9 |

PROCESS CONTROL USING A PID CONTROLLER

Solved with COMSOL Multiphysics 5.1

Inflow 1
1 On the Physics toolbar, click Boundaries and choose Inflow.
2 Select Boundary 1 only.
3 In the Settings window for Inflow, locate the Concentration section.
4 In the c0,c text field, type c_in_inlet.

Inflow 2
1 On the Physics toolbar, click Boundaries and choose Inflow.
2 Select Boundary 7 only.
3 In the Settings window for Inflow, locate the Concentration section.
4 In the c0,c text field, type c_in_top.

Outflow 1
1 On the Physics toolbar, click Boundaries and choose Outflow.
2 Select Boundary 13 only.
DEFINITIONS

Next, add a probe to sample the concentration and its time derivative at the point x =
0, y = -0.002.
1 On the Definitions toolbar, click Probes and choose Domain Point Probe.
2 In the Settings window for Domain Point Probe, locate the Point Selection section.
3 In row Coordinates, set y to -0.002.
4 In the Model Builder window, expand the Component 1 (comp1)>Definitions>Domain
Point Probe 1 node, then click Point Probe Expression 1 (ppb1).
5 In the Settings window for Point Probe Expression, type c_mp in the Variable name

text field.
6 Click Replace Expression in the upper-right corner of the Expression section. From

the menu, choose Component 1 (comp1)>Transport of Diluted Species>c Concentration.

7 In the Model Builder window, under Component 1 (comp1)>Definitions right-click
Domain Point Probe 1 and choose Point Probe Expression.
8 In the Settings window for Point Probe Expression, type ct_mp in the Variable name

text field.
9 Locate the Expression section. In the Expression text field, type ct.

10 |

PROCESS CONTROL USING A PID CONTROLLER

Solved with COMSOL Multiphysics 5.1

Variables 1
1 On the Definitions toolbar, click Local Variables.
2 In the Settings window for Variables, locate the Variables section.
3 In the table, enter the following settings:
Name

Expression

u_in_ctrl

nojac(k_P_ctrl*(c_m
p-c_set)+k_I_ctrl*I
[mol*s/
m^3]+k_D_ctrl*ct_mp
)

Unit

Description
Velocity, controlled
inlet

The nojac operator ensures that the above expression gives no Jacobian
contribution. In practice, this means that the control velocity will always be
evaluated based on the previous time step. This is necessary to avoid evaluation of
an implicit time derivative in the inlet condition, which is not supported in the time
dependent solver.
Moreover, 'I' refers to the time integral in Equation 1, which you define next.
GLOBAL ODES AND DAES (GE)

Global Equations 1
1 In the Model Builder window, expand the Component 1 (comp1)>Global ODEs and
DAEs (ge) node, then click Global Equations 1.
2 In the Settings window for Global Equations, locate the Global Equations section.
3 In the table, enter the following settings:
Name

f(u,ut,utt,t) (1)

Initial value
(u_0) (1)

Initial value
(u_t0) (1/s)

Description

It-(c_mp-c_set)

Time integral
term

MESH 1

1 In the Model Builder window, under Component 1 (comp1) click Mesh 1.

2 In the Settings window for Mesh, locate the Mesh Settings section.
3 From the Element size list, choose Finer.
4 Click the Build All button.

11 |

PROCESS CONTROL USING A PID CONTROLLER

Solved with COMSOL Multiphysics 5.1

STUDY 1

Use a parametric sweep to solve for two different values of the proportional parameter,
k_P.

Parametric Sweep
1 On the Study toolbar, click Parametric Sweep.
2 In the Settings window for Parametric Sweep, locate the Study Settings section.
3 Click Add.
4 In the table, enter the following settings:
Parameter name

Parameter value list

k_P_ctrl

0.1 0.5

Parameter unit

Step 1: Time Dependent

1 In the Model Builder window, under Study 1 click Step 1: Time Dependent.
2 In the Settings window for Time Dependent, locate the Study Settings section.
3 In the Times text field, type range(0,0.05,1) range(1.1,0.1,6).

Solution 1
1 On the Study toolbar, click Show Default Solver.
2 In the Model Builder window, expand the Solution 1 node, then click Time-Dependent
Solver 1.
3 In the Settings window for Time-Dependent Solver, click to expand the Time
stepping section.
4 Locate the Time Stepping section. From the Method list, choose Generalized alpha.
5 From the Steps taken by solver list, choose Intermediate.

This forces the solver to take at least one step in each of the time intervals you
specified.
6 Find the Algebraic variable settings subsection. From the Error estimation list, choose
Exclude algebraic.
7 On the Study toolbar, click Compute.
TABLE

Go to the Table window.

12 |

PROCESS CONTROL USING A PID CONTROLLER

Solved with COMSOL Multiphysics 5.1

RESULTS

Velocity (spf)
The presence of the derivative term leads to a warning message from the solver. As
already mentioned in the introduction, this term is difficult to determine and also
sensitive to noise, so it is often set to 0.

Concentration (tds)
Add a streamline plot of the velocity to the default surface plot that shows the
concentration at the end of the simulated time span (Figure 2). Study the solution at
t = 0.05 s and t = 2 s.
1 In the Model Builder window, under Results right-click Concentration (tds) and

choose Streamline.
2 In the Settings window for Streamline, locate the Streamline Positioning section.
3 From the Positioning list, choose Magnitude controlled.
4 In the Density text field, type 10.
5 On the Concentration (tds) toolbar, click Plot.
6 In the Model Builder window, click Concentration (tds).
7 In the Settings window for 2D Plot Group, locate the Data section.
8 From the Time (s) list, choose 0.05.
9 On the Concentration (tds) toolbar, click Plot.
10 Click the Zoom Extents button on the Graphics toolbar.
11 From the Time (s) list, choose 3.
12 On the Concentration (tds) toolbar, click Plot.

1D Plot Group 4
Plot the PID-controlled inlet velocity (Figure 3).
1 In the Model Builder window, under Results click 1D Plot Group 4.
2 In the Settings window for 1D Plot Group, locate the Plot Settings section.
3 Select the y-axis label check box.
4 In the associated text field, type u_in,ctrl (mm/s).
5 In the Model Builder window, expand the 1D Plot Group 4 node, then click Global 1.
6 In the Settings window for Global, click Replace Expression in the upper-right corner

of the y-axis data section. From the menu, choose Component

1>Definitions>Variables>u_in_ctrl - Velocity, controlled inlet.

13 |

PROCESS CONTROL USING A PID CONTROLLER

Solved with COMSOL Multiphysics 5.1

7 Click to expand the Legends section. Find the Include subsection. Clear the
Description check box.
8 On the 1D Plot Group 4 toolbar, click Plot.

Proceed to plot the concentration at the measurement point as a function of time

(Figure 3).

1D Plot Group 6
1 In the Model Builder window, right-click 1D Plot Group 4 and choose Duplicate.
2 In the Settings window for 1D Plot Group, click to expand the Title section.
3 From the Title type list, choose Manual.
4 In the Title text area, type Concentration, measurement point.
5 Locate the Plot Settings section. In the y-axis label text field, type c_mp
(mol/m³).

6 In the Model Builder window, expand the 1D Plot Group 6 node, then click Global 1.
7 In the Settings window for Global, click Replace Expression in the upper-right corner

of the y-axis data section. From the menu, choose Component 1>Definitions>c_mp Probe variable c_mp.
8 On the 1D Plot Group 6 toolbar, click Plot.

The resulting plot should look like that in the lower panel of Figure 3.

14 |

PROCESS CONTROL USING A PID CONTROLLER