Portal de Periódicos da CAPES

Application of MATLAB Functions for Time Domain Simulation of Systems With Lines With Fluid Transients

2005; ASM International; Volume: 127; Issue: 1 Linguagem: Inglês

10.1115/1.1852488

1528-901X

Patompong Wongputorn, David A. Hullender, Robert L. Woods, John King,

Fluid Dynamics Simulations and Interactions

Contributed by the Fluids Engineering Division for publication in the JOURNAL OF FLUIDS ENGINEERING. Manuscript received by the Fluids Engineering Division July 15, 2003; revised manuscript received September 13, 2004. Review conducted by Y. Tsujimoto. Hydraulic and pneumatic fluid transmission lines are often internal components within a total dynamic system to be analyzed by time domain simulation. This paper introduces a new and simple approach utilizing MATLAB® computational tools for including the model for the line fluid transients in the model for the overall system. The fluid transients in the lines are based on a laminar flow distributed parameter model, which includes nonlinear frequency dependent viscous friction terms as well as heat transfer effects in gas lines.Time domain simulation of a system containing fluid line transients is quite complex if the transfer functions for the fluid transients consist of hyperbolic Bessel functions, which is the case for the “Dissipative” model 1234567891011121314151617181920212223242526. The total system simulation model formulation procedure is to first combine the equations for the system, including the transfer functions for the fluid transients, into a single transfer function for the system output variable of interest. By obtaining the output variable transfer function, it is a one step MATLAB® process to obtain the simulation for any system input of interest. For most systems, the “solve” algorithm in MATLAB® is needed to obtain this transfer function in terms of symbolic variables representing the system variables and the nonlinear transfer functions for the fluid transients. A least-squares curve fit algorithm named ‘invfreqs’ in the MATLAB®Signal Processing Toolbox is used to formulate a rational polynomial approximation for this complex transfer function in the frequency domain 272829303132. The approximation is very accurate over a designated frequency range corresponding to the bandwidth of the components of the total system and to the frequency content of the input disturbances. As will be demonstrated, the advantage of this new approach is the ability to include the distributed parameter model for the fluid transients without questionable assumptions associated with friction approximations and lumped parameter modeling techniques associated with fluid transients. The model formulation and approximation algorithm is demonstrated for a simple line terminating into a fluid tank and then for a complex hydraulic braking system containing several hydraulic lines 33. Consider the schematic of a fluid transmission line shown in Fig. 1. It is assumed that the line is a tube with uniform and rigid circular cross section. Also, it is assumed that the fluid media (liquid or gas) flows though the line in one dimension with mean speed that is relatively slow compared to the speed of sound in the fluid. In Fig. 1, Pa and Qa denote the pressure and volumetric flow at the left end of the line. Also, Pb and Qb denote the pressure and volumetric flow at the right end of the line. The direction of the arrows represents the direction of the flows when positive. Details of the derivation of the so-called “dissipative” distributed-parameter model for this transmission line are presented in numerous publications 11 and most recently in the work of Nursilo 19. A particularly useful input–output expression representing the solution for this distributed-parameter model is presented in Eqs. (1) and (2), i.e., (1)Qas=C1sPas−C2sPbs(2)Qbs=C2sPas−C1sPbswhere (3)C1s=cosh ΓsZssinh Γs(4)C2s=1Zssinh Γs(5)Γs=Dns r2ν 1+γ−1Bσ1−B1/2(6)Zs=Z01−B1+γ−1BσThe Bessel function ratios, Bσ and B are as follows: (7)Bσs=2J1jσsr2/νjσsr2/νJ0jσsr2/ν(8)Bs=2J1jsr2/νjsr2/νJ0jsr2/νwhere, J0 and J1 are zero and first-order Bessel functions of the first kind. Note that the Bessel functions are functions of complex variables. Also, it is important to note that for a liquid, the specific heat ratio γ is one and thus, Bσ drops out of Eqs. (5) and (6). The next sections of this paper demonstrate ways to incorporate these equations for line fluid transients into a model for a total system. Consider the hydraulic fluid line connected to a tank (fluid capacitance) shown in Fig. 2. There is a valve represented by a linear fluid resistance R at the connection to the tank. The input to this system is the pressure Pa and the output of interest is the pressure in the tank Pt. This first example represents a problem where the transfer function relationship between the input and the output can easily be obtained by manually combining the equations. The equation for a lumped fluid capacitance 34 for the tank is (9)Qbs=CtsPtsThe equation for a lumped linear fluid resistance 34 is (10)Pbs=RQbs+PtsUsing Eq. (2) for the fluid transients in the line gives (11)Qbs=C2sPas−C1sPbsEquations (91011) represent three equations and three unknowns Qbs,Pbs, and Pts. Solving for Pts gives (12)Pts=HsPaswhere (13)Hs=C2s[1+C1sR]Cts+C1sThis transfer function is a very nonlinear function of the Laplace variable s. It is desired to approximate Hs by a simple ratio of rational polynomials in s representing a linear differential equation. In general, the approach to achieving a transfer function approximation is to match the frequency response of the original transfer function, over a frequency range of interest, with the frequency response of a rational polynomial transfer function, which represents a linear ordinary differential equation 32. This approach utilizes a least-squares curve fit algorithm named “invfreqs,” provided in MATLAB® Signal Processing Toolbox, to obtain the polynomial coefficients associated with the transfer function frequency response curve fits. This m-file, “invfreqs,” employs the least squares technique of Levi 27 along with an iteration algorithm, which utilizes the damped Gauss–Newton method. Before generating a rational polynomial representation for Hs, it necessary to determine the frequency range of interest for the system dynamics. The lower frequency is usually close to zero to preserve the dc gain of the transfer function. The upper bound on the frequency range is determined by the largest of either the frequency content of the input Pat, by c/L, or by estimating the inverse of the time constant associated with pressurizing the tank, i.e. (14)1τ=1CtRsf+RUsing MIL-H-5060 hydraulic fluid at 26.7°C (80°F), the approximate fluid properties 34 are β=1.896×109 N/m2,ν=2.0×10−5 m2/s, and ρ=875 kg/m3. If the tank volume, V, is 0.00142 m3, the line length, L, is 5 m, the inside radius, r, of the line is 0.005 m, and the valve resistance, R, is 1.776e8 Ns/m5, then Eq. (14) gives 0.0004 s which corresponds to a frequency of 2500 rad/s. An estimate of the order of the transfer function needed to accurately approximate Hs over the desired frequency range is determined from an examination of the frequency response of Hs. At least two orders are needed for each apparent second order mode. Examination of the frequency response reveals that there are at least three second order modes over this frequency range; thus, the order of the approximation should be at least 8 and possibly 10. A tenth order transfer function was chosen for this example. The procedure of achieving an accurate curve fit is a significant contribution of this paper; so, the MATLAB® m-file, which uses “invfreqs” to generate the approximation for Hs, is provided in the Appendix and on the Author’s web site 33. It is important to note that the number of frequency points and spacing must be determined by trial and error. Since all of the computations are in MATLAB®, the order is probably not a major issue since the transfer function will be stored in memory for further simulation and analysis. For this case, the result is a tenth order rational polynomial transfer function representing a linear ordinary differential equation, i.e. Hs≈−3.398e6s9+1.041e10s8−⋯+1.21e33s+1.238e35s10+2.762e9s9+⋯+7.631e30s2+1.281e33s+1.235e35As shown in Figs. 3 and 4, the inverse frequency algorithm generates a very accurate linear model out to the specified maximum frequency of interest, 2500 rad/s. The time response corresponding to a step pressure input at the left end of the line is shown in Fig. 5; note the pure time delay at the beginning of the response associated with the speed of sound and the time required for the first pressure changes to reach the tank. For a second example, consider the schematic of a hydraulic brake actuation circuit shown in Fig. 6 for a front-end loader. When the brake valve is opened by pressing on the brake pedal, hydraulic fluid passes through line “a” of length La to the “T” connector which diverts fluid to the front and rear brake cylinders. As shown, there are seven different hydraulic lines, which must be modeled in order to analyze the time domain response of the pressures in the brake cylinders. A typical reason for analyzing this system might be to determine how the time required to pressurize the brake cylinders is affected by the diameter of the lines and/or by the temperature of the oil. Using Eqs. (1) and (2) for each of the line sections in Fig. 6, results in the following equations: Line “a”: (15)Q2=C2aP1−C1aP2Line “b”: (16)Q3=C1bP2−C2bP3Q4=C2bP2−C1bP3Line “c” front: (17)Q5=C1cP3−C2cP4Q6=C2cP3−C1cP4Line “c” rear: (18)Q9=C1cP5−C2cP6Q10=C2cP5−C1cP6Line “d”: (19)Q7=C1dP2−C2dP5Q8=C2dP2−C1dP5For the “T” connectors, the equations are (20)Q2=Q3+Q4(21)Q4=2Q5(22)Q8=2Q9For the brake cylinders, the equations are (23)Q6=CbsP4(24)Q10=CbsP6where the lumped parameter capacitance is given by (25)Cb=Vβ+A2KEquations (15161718192021222324) represent 14 equations with 14 unknowns: P2,P3,P4,P5,P6,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9, and Q10. The input is P1. Because of the complexity and number of equations, as compared to the first example, where the solution was easily obtained manually, it is advantageous to use symbolic math and the “solve” algorithm in MATLAB® to get the solution. The MATLAB® m-file which uses “solve” and “invfreqs” for solving these equations and then generating the rational polynomial transfer function for the total system is provided in the Appendix. When reviewing this m-file, it is important to note the symbolic variable substitution procedure for the hyperbolic Bessel functions and parameters. Also, the frequency increments will not always be obvious and must be determined by trial and error attempts to get reasonable curve fits at both the upper and lower frequency bounds. Using the total line length from the brake valve to the rear brake cylinder divided into the speed of sound gives a frequency of 134 rad/s. The inverse of the time constant based on the steady flow resistance and brake cylinder capacitance is 45 rad/s. The frequency range of interest is selected to be 1 rad/sec to 300 rad/sec and the order of the transfer function approximation is chosen to be 10. For the designated frequency range and system order, the magnitude and phase frequency response comparisons are shown in Figs. 7 and 8. The curve fit is observed to be extremely close to the exact out to 300 rad/s. This transfer function between the input pressure P1 and the rear brake cylinder pressure P6 is shown in Eqs. (26) and (27). Note that this transfer function represents a linear ordinary differential equation that is good for any assumed input P1s.(26)P6s=GsP1swhere (27)Gs≈5.987e13s9−8.70e14s8+⋯+3.674e35s+5.45e36s10+1.35e17s9+⋯+3.86e35s+5.45e36Typical simulation results using this transfer function are shown in Figs. 9 and 10 for a step input. As shown in Fig. 9, the performance of the brake system is very sensitive to the temperature of the oil with 0.00635 m (1/4 in.) diameter lines. The viscosity of ISO 46 oil at 10°C (50°F) is so high that it takes almost 4 s for the brake cylinder pressure to reach the desired level. For 79.4°C (175°F), the desired pressure is achieved in less than 0.2 s. As shown in Fig. 10, the performance of the brake system is sufficiently fast, even at the lower temperature, for a line diameter of 0.0127 m (1/2 in.). This second example demonstrates the simplicity of using the MATLAB® Symbolic Toolbox for solving very complex simultaneous equations. Because of the nonlinear properties of the model for the fluid line transients as represented by C1 and C2, it would be essentially impossible to formulate a simulation model for this system without having to simplify our friction model in the equations for the line fluid transients or resort to lumped parameter models for the lines or limit our analysis to methods requiring periodic inputs. For problems where a fluid transmission line is an internal component within a total system to be simulated in the time domain, it is often desirable to be able to accurately model the fluid transients in the line with partial differential equations that include frequency dependent friction terms and heat transfer effects for gases. Thus, accuracy issues associated with friction model approximations and associated with lumped parameter modeling errors are avoided. An algorithm has been introduced that generates very accurate rational polynomial transfer functions for a total system in which one or more lines are internal components. The transfer functions are generated in the MATLAB® working environment using symbolic variable solution algorithms. The contribution of this paper is the procedure for using MATLAB® to simulate dynamic systems with fluid lines in which the fluid transients are significant in the frequency range of interest without having to approximate the fluid transient friction model. A= Brake cylinder piston area, m2c= Speed of sound, m/s Ct= Tank capacitance, V/β 34Dn= Dimensionless dissipative number, νL/cr2J= −1K= Brake cylinder piston return spring constant, N/m L= Line length, m Pt= Tank pressure r= Internal radius of the line, m R= Linear approximation for valve flow resistance, ∂P/∂QRsf= Steady laminar flow line resistance, 8ρνL/πr434S= Laplace variable t= time, s τ = time constant, s V= Tank volume, m3ZO= Characteristic impedance constant, ρc/πr2β = Bulk modulus of fluid, N/m2γ = Ratio of specific heats ρ = Mass density of fluid, kg σ = Prandtl number for the fluid ν = Kinematic viscosity of fluid, m2/s ω = Frequency, rad/s function [Ht]=Example1 % Fluid line terminating into a tank (fluid capacitance). warning off Vis=2.0e−5; r=0.005; Be=1.896e9; Den=875; c=sqrtBe/Den;Zo=Den*c/pi*r∧2; V=0.00142; L=5; Ct=V/Be;R=1.776e8; syms s C1 C2 B=2*besse1j1,j*sqrtr∧2*s/Vis/j*sqrtr∧2*s/Vis*besse1j0,j*sqrtr∧2*s/Vis;Z=Zo/sqrt1−B;Dn=Vis*L/c*r∧2;G=Dn*r∧2*s/Vis/sqrt1−B;H=C2/1+C1*R*Ct*s+C1;%Trans. fun. to be curve fitted. H=subsH,‘C1’,coshG/Z* sinhG;H=subsH,‘C2’,1/Z* sinhG; % Generate frequency data for the curve fit. Note, the last value is the range of interest. w=[10:1:195,200:10:980,1000:20:2500];N=lengthw; % Generate a greater frequency range for an accuracy comparison wc=[10:1:195,200:10:980,1000:20:6000];NC=lengthwc; for k=1:Nsw=j*wck;TFk=subsH,s,sw;TFck=TFk; end N1=N+1; for k=N1:NC %Additional frequencies for accuracy comparisons. sw=j*wck;TFck=subsH,s,sw; % Additional data for the accuracy comparison. end [Num,Den]=invfreqsTF,w,9,10,[ ],100;Ht=tfNum,DenstepHt,0.03 % Generate Freq. Resp. plots to determine the accuracy of curve fit. figure Hc=freqsNum,Den,wc;MH=20* log 10absHc;AH=angleHc*180/pi; MTF=20* log 10(abs(TFc)); ATF=angle(TFc)*180/pi; semilogxwc,MH,‘r’,wc,MTF,‘b’; title(‘Magnitude Comparison Plots’); xlabel(‘Frequency rad/sec’); ylabel(‘Decibels’); figure semilogx(wc,AH,‘r’,wc,ATF,‘b’); title(‘Phase Comparison Plots’); xlabel(‘Frequency rad/sec’); ylabel(‘Degrees’); function [G]=example2 % Hydraulic brake system for a front-end loader. Vis=10e−6; Den=820; % ISO 46 Oil at 175°F % Vis=300e−6; Den=880; % ISO 46 Oil at 50°F warning off; r=0.0032;Be=1.7235e8;c=sqrtBe/Den;Zo=Den*c/pi*r∧2;V=2.46e−6;A=0.0111;K=3.78e6; La=.559; Dna=Vis*La/c*r∧2; Ld=2.134; Dnd=Vis*Ld/c*r∧2; Lb=0.762; Dnb=Vis*Lb/c*r∧2; Lc=0.7112; Dnc=Vis*Lc/c*r∧2; syms s B=2*besse1j1,j*sqrtr∧2*s/Vis/(j*sqrtr∧2*s/Vis*B=2*besse1j1,j*sqrtr∧2*s/Vis/j*sqrtr∧2*s/Vis*besse1j0,j*sqrtr∧2*s/Vis;Z=Zo/sqrt1−B;Ga=Dna*r∧2*s/Vis/sqrt1−B;Gd=Dnd*r∧2*s/Vis/sqrt1−B;Gc=Dnc*r∧2*s/Vis/sqrt1−B;Gb=Dnb*r∧2*s/Vis/sqrt1−B;Sol=solve‘Q2=C2a*1−C1a*P2’,‘Q3=C1b*P2−C2b*P3’,‘Q4=C2b*P2−C1b*P3’,‘Q5=C1c*P3−C2c*P4’,…‘Q6=C2c*P3−C1c*P4’,‘Q7=C1d*P2−C2d*P5’,‘Q8=C2d*P2−C1d*P5’,‘Q9=C1c*P5−C2c*P6’,…‘Q10=C2c*P5−C1c*P6’,‘Q2=Q3+Q4’,‘Q4=2*Q5’,‘Q8=2*Q9’,‘Q6=Cb*s*P4’,‘Q10=Cb*s*P6’,…‘P2,P3,P4,P5,P6,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9,Q10’;H=Sol.P6;collectH,s;H=subsH,‘C1a’,coshGa/Z* sinhGa;H=subsH,‘C2a’,1/Z* sinhGa;H=subsH,‘C1d’,coshGd/Z* sinhGd;H=subsH,‘C2d’,1/Z* sinhGd;H=subsH,‘C1c’,coshGc/Z* sinhGc;H=subsH,‘C2c’,1/Z* sinhGc;H=subsH,‘C1b’,coshGb/Z* sinhGb;H=subsH,‘C2b’,1/Z* sinhGb;H=subsH,‘Cb’,V/Be+A∧2/K;w=[1:.1:9.9,10:.2:19.8,20:1:99,100:5:195,200:20:300];N=lengthw;wc=[1:.1:9.9,10:.2:19.8,20:1:99,100:5:195,200:20:1000];NC=lengthwc; for k=1:N;sw=j*wck;TFk=subsH,s,sw;TFck=TFk;end N1=N+1; %Additional frequencies generated for accuracy comparisons. for k=N1:NC;sw=j*wck;TFck=subsH,s,sw; end %Create a transfer function by curve fitting. [Num,Den]=invfreqsTF,w,9,10,[ ],100;G=tfNum,Den;stepG figure Hc=freqsNum,Den,wc;MH=20* log 10absHc;AH=angleHc*180/pi;MTF=20* log 10absTFc;ATF=angleTFc*180/pi;semilogxwc,MH,‘r’,wc,MTF,‘b’; title(‘Magnitude Comparison Plots’); xlabel(‘Frequency rad/sec’); ylabel(‘Decibels’); figure semilogxwc,AH,‘r’,wc,ATF,‘b’ title(‘Phase Comparison Plots’); xlabel(‘Frequency rad/sec’); ylabel(‘Degrees’).