How to model PI flow control with throttling (butterfly) valve in hydraulic centriful pump Simscape/Simulink model?

I want to model PI-controlled butterfly throttling valve to control the flow. The control variable is an early flow measurement. If this value generates a negative value when subtracting this from the setpoint, the valve must open up a bit. If this error is positive, the valve should close some more. The valve orifice has a maximum diameter of 0.19m, and is the inverse of the valve disk diameter as the disk regulates the flow by opening and closing. The control segment of my model is shown below:

For a certain valve I have about as much data as I need:

• Kv value vector over valve opening in degrees, with an approximation by formula . See plot below for normalized values.

• Maximum diameter of orifice, 0.19 m
• Pressure drop can be read from model

I have two initial questions:

1)

What is exactly meant with control member position vector? I thought this is a vector the valve disk moves which goes from 0% (0 m, fully open) to 100% (0.19 m, fully closed) that corresponds to either the orifice area vector, or the pressure drop vector 'dp' and volumetric flowrate table 'q(s,dp). However, in the help guide it says that physical signal at port S is compared to the control member position vecter. I do not understand how this position of the member, which in my case is the diameter of the valve disk, can be compared to an PI-tuned error signal from the process. Furthermore, in the linear option, 'Control member position at closed orifice' is equal to 'Control member travel between closed and open orifice' are both equal to the changing diameter of the disk, so why do I need to input them both?

2)

With the known variables, what would be the best way of accurately inputting the data for an accurate butterfly valve model?

[Sam Chak]

Hi @DB

I don't get the same plot that you have shown. Can you check if I did it correctly?

Update: I guess this is what you meant by normalization. I think this is the orange curve though there is mistake on the legends. If you provide the actual data-pair values, it is possible to apply some curve-fitting techniques.

L = 2.95;
x = linspace(0, 100, 10001);
y = (x/100)./(L + (1 - L)*(x/100));
plot(x, y), grid on
xlabel('Valve Opening'), ylabel('Flow')

[DB]

Hi @Sam Chak, indeed, something went wrong with the color but it should be the orange line. I do not have curve fitting apps and couldn't find a clear way to do it another way. This method comes from literature. The actual Kv values (normalized to 100) are:

Kv_100 = [0    0.1999    2.9188    8.2767   15.5538   25.9496   41.1435   60.6557   82.4870  100.0000]

The L value is found using for loops doing trial and error.

[Sam Chak]

Hi @DB

I have updated the Answer based on your data Kv. Let me know if it works out linearly for your PID control design.

[DB]

@Sam Chak thank you for the quick response. However, as I commented earlier, I still don't understand what you mean with 'Pre-conditioned valve opening' and where you get these values from. How do you know you can use the equation

y2 = (18.15*x2.^2 + 1259*x2 - 7800)./(x2.^2 - 190.8*x2 + 1.207e4);

and where do the parameters 18.15,1259 etc. come from? I have the same question for the equations

u  = (95.4*(x2 + 6.59853))./(x2 - 18.15) - (0.1*sqrt(-296884*x2.^2 + 3.31379e7*x2 + 5.3784e7))./(x2 - 18.15);
y3 = (18.15*u.^2 + 1259*u - 7800)./(u.^2 - 190.8*u + 1.207e4);.

What is the underlying reason or theory that brought you to this?

Answer 1

[Sam Chak]

Hi @DB

Due to the static nonlinearity of the valve characteristics, it is possible to find an inverse function for the valve input that allows the flow output to respond in a linear fashion. Consequently, you should be able to design the PID controller in a linear fashion.

Update: Based the provided data, you can use the fit() function from the Curve Fitting Toolbox to fit a desired curve (math model) or surface to data. Since you used a rational function previously, I also selected the rational function. It seems that a 2nd-degree rational function (rat22) can fit the data.

% Original data Kv
x1 = linspace(0, 100, 10);
y1 = [0 0.1999 2.9188 8.2767 15.5538 25.9496 41.1435 60.6557 82.4870 100.0000];
[fitobject, gof] = fit(x1', y1', "rat22")

Warning: Start point not provided, choosing random start point.

fitobject = 
     General model Rat22:
     fitobject(x) = (p1*x^2 + p2*x + p3) / (x^2 + q1*x + q2)
     Coefficients (with 95% confidence bounds):
       p1 =       18.15  (-5.708, 42)
       p2 =        1259  (132.1, 2386)
       p3 =       -7801  (-2.324e+04, 7641)
       q1 =      -190.8  (-202.9, -178.6)
       q2 =   1.207e+04  (1.097e+04, 1.316e+04)

gof = struct with fields:
           sse: 1.7081
       rsquare: 0.9999
           dfe: 5
    adjrsquare: 0.9997
          rmse: 0.5845

The R-squared (coefficient of determination) and the plotting result confirms that the 2nd-degree rational function fits the 10-point data of very closely.

% Approximated Kv
x2 = linspace(0, 100, 10001);
y2 = (18.15*x2.^2 + 1259*x2 - 7800)./(x2.^2 - 190.8*x2 + 1.207e4);

figure(1)
plot(x1, y1, '.'), hold on
plot(x2, y2), grid on, hold off
xlabel('Valve Opening, %'), ylabel('Flow, %')
legend('data K_{v}', 'approx. K_{v}', 'location', 'southeast')

From the approximated function, you can use the finverse() function to compute the inverse of function , such that . In short, this g(x) is the inverse valve characteristics.

syms x
f(x) = (18.15*x^2 + 1259*x - 7800)/(x^2 - 190.8*x + 1.207e4);
g = finverse(f);
simplify(g)

ans(x) = 

Physically, it implies that if a certain percentage of flow is desired, the we can use the inverse function to calculate required percentage of valve opening that allows the actual flow output to respond in a linear proportion to the reference .

% Inverse valve characteristics
u  = (95.4*(x2 + 6.59853))./(x2 - 18.15) - (0.1*sqrt(-296884*x2.^2 + 3.31379e7*x2 + 5.3784e7))./(x2 - 18.15);
y3 = (18.15*u.^2 + 1259*u - 7800)./(u.^2 - 190.8*u + 1.207e4);

figure(2)
subplot(2, 1, 1)
plot(x2, y2), grid on
xlabel('Valve Opening u, %'), ylabel('Actual Flow K_{v}, %')
subplot(2, 1, 2)
plot(x2, u), grid on, ylim([0 100])
xlabel('Reference Flow K_{v, ref}, %'), ylabel('Required valve opening u, %')

figure(3)
plot(x2, y3), grid on, ylim([0 100])
xlabel('Reference Flow K_{v, ref}, %'), ylabel('Actual Flow K_{v}, %')

[DB]

Hi @Sam Chak thanks for the answer. This looks as if it could work. However, what do you mean with pre-conditioned valve opening? And how did you get to the values 5900,2000 and 39?

[Sam Chak]

Hi @DB

I have updated my Answer with explanations and also made some corrections on the terms, as I'm unsure what professional terms are used in your industrial valve/process field. The suggested solution is purely mathemathetical theory. There may be some losses in practice.

However, I have found that the authors of "Two-Degree-of-Freedom Flow Rate Control for Pneumatic Valves using Fast Response Flow Meter", also applied the same functional inverse approach. The paper is available to download.

[DB]

@Sam Chak

I do not have access to either the Curve Fitting toolbox for fit(), nor the Symbolic Math Toolbox for finverse(). I could avoid curve fitting by using the general function with .

However, I still don't understand how I can use these values for a variable orifice. Firstly, I have to select one the following option for Orifice (IL):

• Linear - Area vs. control member position
• Tabulated data - Area vs. control member position
• Tabulated data - Volumetric flow rate vs. control member position and pressure drop

From our discussion, I only have some indication of the ratio % valve open vs. % Kv. The third option is the only one that involves the flow rate, but I don't know about the pressure drop.

What is the best option? I have tried all of them, but my valve just swings between 0% and 100% opening with all of them, and simulation time is about a minute per second simulation.

[Sam Chak]

Hi @DB

I understand that the Curve Fitting and Symbolic Math Toolboxes are unavailable. Is it necessary to demonstrate the determination of the fitted model and the inverse function using MATLAB?

However, I'm still unsure about how the rational function x/(L + (1 - L)*x) is applied in the context of (1) Area vs. control member position and (2) Volumetric flow rate vs. control member position and pressure drop.

The purpose of using the inverse function is to allow you to easily design the PID controller under the assumption that the valve actuation behaves as a linear perfect transmission.

[DB]

Hi @Sam Chak.

The inverse of the general function can easily be found algebraically to .

However, you are correct that it does not give a good approximation of the 'dead zones' of a butterfly valve, one of which is seen in my plot of the actual Kv between 0-20%. This equation comes from this paper.

I'm using MATLAB because I want to model in Simscape for an accurate simulation of an actual hydraulic system. I'm told the variable orifice is a good block for valve modeling, but it is difficult so far to translate actual manufacturer data to Simscape. Therefore, I'm also contemplating using some MATLAB function instead.

[Sam Chak]

Hi @DB

It is okay to explore things that you want to test with. That's how we learn things until we find a better one.

I'm unfamiliar with Simscape, because I was traditionally trained to apply pure math to describe objects statically and processes dynamically. However, Simscape has a comprehensive of library blocks for a variety of industrial-oriented and engineering stuff. Students tend to learn things relatively easier with Simscape because the components are easily idendified and they can connect the equipments required for a particular process.

When it comes maniputing the process so that the desired behavior is achieved, it is a different level because the knowledge about the equipment and the process are required. For some processes, it is possible to analytically determine the actuation signal required to achieve the desired result. Otherwise, tuning by AI and auto-optimizer are mostly preferred by students if procedure is fast and accurate. Unfortunately, some tuning results are difficult to be interpreted by humans.

[Henry]

Hello DB,

I'm also in research of pressure control by PID throttle valve, and you case is very similar to mine described as below. The logic is clear but building blocks is hard. Can you please share your simulink model with me? That will be very helpful for me to build up and quick start. [e-mail]

Thank you.

Henry

The system is described as below. MFC, mass flow controller input a gas nitrogen flow as 100sccm into a chamber with an volume of 90 liters. The chamber is monitored with a pressure sensor and controlled by a throttle valve with a defined flow conductance (flow rate, Q)/(Pressure, P) vs. Opening, at the downstream of throttle valve, it's an vacuum pump with a defined pumping speed. So the system is consisted with MFC, chamber, pressure sensor, throttle valve, and vacuum valve. The pressure of the chamber is controlled by a defined PID. The scope or results that I need is dynamic change of pressure with valve opening tuned by different PID parameters. Please help build the blocks, diagrams in SIMULINK. the dynamic governing equation is Q_in-Q_out =dP/dt. * volume of chamber. Q_in is flow into chamber input from MFC, Q_out = Q/P(valve opening)*P (Q/P(valve opening) is defined as the flow conductance of throttle valve), P is the chamber pressure.