COMSOL 6.4 - Process Control Using a PID Controller

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. Since an increased flow of that gas will decrease the concentration, the PID coefficients will have a negative sign for this PID controller.

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 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.

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)

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:

Table 1: Mass-balance Boundary conditions.
Boundary	Constraint
Upper inlet	c = cin,top
Controlled inlet	c = cin,inlet
Outlet	n · (−D∇c) = 0
Walls	N · n = 0

Here c is the concentration; cin,top and cin,inlet are the inlet concentrations (SI unit: mol/m3) for the upper and controlled inlets, respectively; D is the applied diffusivity (SI unit: m2/s); and N is the molar flux (SI unit: mol/(m2·s)).

The model uses the following boundary conditions for the fluid flow:

Table 2: Fluid-flow boundary conditions.
Boundary	Constraint
Upper inlet	u = (0, −vin,top)
Controlled inlet	u = (uin, 0)
Outlet	p0 = 0
Inlet sections	n · u = 0
Walls	u = 0

Here u is the velocity vector (SI unit: 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 the following, which is the standard PID control algorithm available in the PID Controller Setup add-in:

(1)

with the following parameters:

Table 3: PID controller parameters.
parameter	value
cset	0.5 mol/m3
kP	-0.5 m4/(mol·s)
kI	-1 m4/(mol·s2)
kD	-10-3 m4/mol

As mentioned, the negative values of the coefficients reflect the fact that cset is lower than the concentration at the upper inlet, cin,top, and that the purpose of the gas stream at the controlled inlet thus is to reduce the concentration. 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 measurement c. The PID Controller Setup add-in includes the option to filter the derivative part, which reduced the influence of high-frequency noise and therefore makes the derivative part more useful.

Results and Discussion

The two plots in Figure 2 show the oxygen concentration and the velocity streamlines 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 startup, 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.

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 green line represents the results for a kP value of 0.5 m4/(mol·s) while the blue line corresponds to kP equal to 0.1 m4/(mol·s). 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/(mol·s) (green) and kP = 0.1 m4/(mol·s) (blue).

COMSOL_Multiphysics/Multiphysics/pid_control

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select > > .

Click .

In the tree, select > .

Click .

Click

In the tree, select > .

Click

Global Definitions

Parameters 1

In the window, under click .

In the window for , locate the section.

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_inlet	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
D	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/(s·mol)	Proportional parameter
k_I_ctrl	-1[m^4/(mol*s^2)]	-1 m4/(s²·mol)	Integral parameter
k_D_ctrl	-1e-3[m^4/mol]	-0.001 m4/mol	Derivative parameter

As mentioned in the Model Definition section, the PID parameters are negative because the setpoint concentration is lower than that at the upper inlet.

Geometry 1

Create the geometry. To simplify this step, insert a prepared geometry sequence.

In the toolbar, click and choose .

Browse to the model’s Application Libraries folder and double-click the file pid_control_geom_sequence.mph.

In the toolbar, click

Click the

button in the toolbar.

In the window, under click .

Definitions

Next, add a probe to sample the concentration and its time derivative at the point x = 0, y = −0.002.

Domain Point Probe 1

In the toolbar, click

and choose .

In the window for , locate the section.

In row , set to -0.002.

Point Probe Expression 1 (ppb1)

In the window, expand the node, then click .

In the window for , type c_mp in the text field.

Click in the upper-right corner of the section. From the menu, choose > > > .

Point Probe Expression 2 (ppb2)

In the window, right-click and choose .

In the window for , type ct_mp in the text field.

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

Proceed to import the PID Controller add-in and set up the control algorithm.

In the toolbar, click

and choose .

Add-in Libraries

In the window, in the tree, select the checkbox for the node > (if it is not already selected).

Click

to load the add-in and close the window.

Developer

In the toolbar, click

and choose > .

Global Definitions

PID Controller Setup 1

In the window, under click .

In the window for , locate the section.

Enter the following settings:


Parameter	Value
Proportional gain	k_P_ctrl
Integral gain	k_I_ctrl
Derivative gain	k_D_ctrl
Reference value	c_set

Keep the default values for the remaining parameters.

Click .

Materials

For the physics setup, you need to specify the density, ρ, and the dynamic viscosity, μ, of the fluid. For this purpose, define a material.

Fluid

In the window, under right-click and choose .

In the window for , type Fluid in the text field.

Select Domain 1 only.

Locate the section. In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Density	rho	1.2[kg/m^3]	kg/m³	Basic
Dynamic viscosity	mu	3e-5[Pa*s]	Pa·s	Basic

You are now ready to set up the physics of the model.

Laminar Flow (spf)

Inlet 1

In the toolbar, click

and choose .

Select Boundary 1 only.

In the window for , locate the section.

In the U0 text field, type pid1.u_in_ctrl. This is the inlet velocity defined by the PID controller.

Inlet 2

In the toolbar, click

and choose .

Select Boundary 7 only.

In the window for , locate the section.

In the U0 text field, type v_in_top.

Outlet 1

In the toolbar, click

and choose .

Select Boundary 13 only.

Wall 2

In the toolbar, click

and choose .

Select Boundaries 2, 3, 6, and 8 only.

In the window for , locate the section.

From the list, choose .

Transport of Diluted Species (tds)

Fluid 1

In the window, under > click .

In the window for , locate the section.

From the u list, choose .

Locate the section. In the Dc text field, type D.

Initial Values 1

In the window, click .

In the window for , locate the section.

In the c text field, type c00.

Inflow 1

In the toolbar, click

and choose .

Select Boundary 1 only.

In the window for , locate the section.

In the c0,c text field, type c_in_inlet.

Inflow 2

In the toolbar, click

and choose .

Select Boundary 7 only.

In the window for , locate the section.

In the c0,c text field, type c_in_top.

Outflow 1

In the toolbar, click

and choose .

Select Boundary 13 only.

Mesh 1

In the window, under click .

In the window for , locate the section.

From the list, choose .

Click

Study 1

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

Parametric Sweep

In the toolbar, click

In the window for , locate the section.

Click

In the table, enter the following settings:


Parameter name	Parameter value list
k_P_ctrl (Proportional parameter)	-0.1 -0.5

Step 1: Time Dependent

In the window, click .

In the window for , locate the section.

In the text field, type range(0,0.05,1) range(1.1,0.1,6).

Solution 1 (sol1)

In the toolbar, click

In the window, expand the node, then click .

In the window for , click to expand the section.

From the list, choose .

This forces the solver to take at least one step in each of the time intervals you specified.

In the toolbar, click

Results

Concentration (tds)

The default plot group contains a surface plot that shows the concentration at the end of the simulated time span, as well as an arrow surface of the flux. Study the solution at t = 0.05 s and t = 2 s.

First, however, add streamlines and adjust the streamline positioning so that the streamline density reflects the velocity magnitude.

Arrow Surface 1

In the window, expand the node.

Right-click and choose .

Concentration (tds)

In the toolbar, click