- Email: [email protected]
MODELING OF PIEZOELECTRICALLY ACTUATED FUEL INJECTORS
MODELING OF PIEZOELECTRICALLY ACTUATED FUEL INJECTORS Dirk Mehlfeldt ∗,1 Hartmut Weckenmann ∗∗ Dr.-Ing. habil. G¨ unter St¨ ohr ∗∗∗ ∗ DaimlerChrysler AG, Stuttgart, Germany. Phone: +49 (0)711 17-52196. [email protected] ∗∗ DaimlerChrysler AG, Stuttgart, Germany. Phone: +49 (0)711 17-58068. [email protected] ∗∗∗ Universit¨ at Siegen, Germany. Phone: +49 (0)271 740-2235. [email protected]
Abstract: In this paper, four mathematical modeling approaches for the electromechanical behaviour of piezoelectrically actuated fuel injectors are discussed. Models with lumped and distributed parameters are investigated together with linear and nonlinear piezoelectric submodels. For all four resulting sets of nonlinear differential equations, numerical integration methods are mentioned. The simulation results of the parameterized models are compared to measurements of a real injection valve. Keywords: Hysteresis, Mathematical models, Numerical simulation, Partial differential equations, Time-varying systems Classification: A (Methods and tools for multi-disciplinary systems development)
1. INTRODUCTION Together with a growing competition in fuel economy aggravated emission restrictions force the automotive engineers to enhance consecutively the thermodynamic processes of internal combustion engines. One of the most recent developments is the gasoline engine with direct injection and stratified combustion. Within this engine concept, the injector is a key component. Figure 1 outlines such a piezoelectrically actuated injector (Gottlieb et al., 2001). The piezoelectric stack actuator with its electrical connections is located in the center of the injector. On the left hand side, it is coupled to a metal element which is in contact to the housing. To the right, the actuator acts against the needle which controls the hydraulic valve at the very right 1
Supported by Conti Temic microelectronic GmbH
Fig. 1. Piezoelectrically actuated injector end. Whenever the actuator is being charged, the induced strain pushes the needle to the right so that the valve is opened. If the actuator is being discharged, a spring pushes the needle back, and the valve closes. Due to the stiff couplings and due to the small inertia of the involved masses, it takes less than 200µs to switch from close to open or back. The common mathematical modeling approach for injectors is an ordinary linear differential equation (Kasper et al., 1997). This approach neither reflects the switching behaviour of the valve nor
433
x˙ i = Ai x + Bi u
the nonlinearity of piezoelectric actuators or the distributed character of the real mechanics. This article introduces and compares four new and improved mathematical modeling approaches. Each mathematical model is parameterized, and the simulation results are compared to measurements. 2. LUMPED PARAMETERS WITHOUT HYSTERESIS At first, the relevant injector mechanics are considered as lumped masses with discrete interactions. Choosing two lumped masses m1 and m2 for both ends of the actuator, together with the consecutive elements and one lumped mass m3 for the valve part of the needle, the injector mechanics become an arrangement as shown in figure 2. m1 k1
FP k2
r1
hN
m2
r3
Fig. 2. Injector mechanics considered as lumped FP is the piezoelectrically induced stress within the actuator. FN is a static force which closes the valve securely when the actuator is discharged. The discrete elements k1 , k2 and k3 represent the stiffness of the housing, the actuator and the needle, respectively. Elements labeled with r stand for velocity proportional dampings. It turned out that the investigated fuel valve has a relatively stiff housing. In this special case, m1 never changes position in the modeled arrangement. Thus, it is omitted in the further examinations. The switching behaviour of the valve requires a distinction of the two cases for the two hydraulic states. Let hN denote the valve clearance. If the valve is closed (hN ≤ 0), the needle is in contact to the valve seat, so that only m2 is free to move. If the valve is opened (hN > 0), the positions of m2 and m3 are variable. Table 1 compares the two resulting simplified arrangements. Table 1. Simplified arrangements Valve state Valve closed hN ≤ 0 (i = 1) Valve opened hN > 0 (i = 2)
Arrangement FP
m2
k2 r2 FP k2 r2
k3 m2
m3 F N k3 r3
Both simplified arrangements can be described separately by sets of ordinary linear differential equations
(1)
i
y = Cx +Du
(2)
with index i ∈ [1, 2] denoting the respective valve state. The input vector u = (IP FN )T is composed by the actuator current IP and the closing force on T the needle FN . The output vector y = (UP hN ) consists of the computed actuator voltage UP and the calculated valve opening hN . As can be seen in table 2, the elements of the state vector x are the actuator charge and the kinetic variables of the lumped masses. Table 2. Elements of the state vector x Element xi1
Physical meaning Electrical actuator charge QP Actuator position hP Actuator velocity h˙
xi2 xi3
P
xi4
m3 F N
k3 r2
i
Valve clearance hN Needle velocity h˙ N
xi5
With these assumptions, the matrices in the equations 1 and 2 become ⎛ 0 0 0 0 0 ⎞ 0 0 1 0 0 r2 ⎜ 1 ∂FP 1 ∂FP − k − k ⎟ 1 A =⎝ − 0 0 ⎠ 2 3 m2 ∂Q m2 ∂h m2 ⎛ 2
A =
0 0
0 0 1 ∂FP
⎜ ⎜ ⎜m2 ⎝
∂QP 0
P
P
1 m2
∂FP ∂hP
1
B
1 0 0 0
0 0
− k2 − k3
⎛1
⎞
0
2
B
−
⎜0 =⎝ 0 0 −
0 0 0 0 1 m3
0 1 r2 m2 0 0
m3
=⎝ 0 0 ⎠ 0 0 0 0
0 0
0 k3
0
⎛
0 0
⎞ ⎟ ⎠
⎛ ⎜ ⎜ C=⎜ ⎝
0 0
0 0 k3
0 0
0 0
0 m2 0 1 k3 r3 − − m3 m3
∂UP ∂QP ∂UP ∂hP 0 0 0
0
⎞ ⎟ ⎟ ⎟ ⎠
⎞T
⎟ ⎟ ⎟ 0⎠ 0
D=0 .
1 0
In this section, a linearized piezoelectricity is considered. Thus, the partial derivatives in the matrices above are to be determined using the linear relations 1 1 s lP d QP − h 2 εs − d2 A NP NP εs − d2 P P d AP 1 ε Q − h FP = NP εs − d2 P εs − d2 lP P
UP =
(3) (4)
between the voltage UP , the charge QP , the force FP and the displacement hP of the actuator (Voigt, 1966). Thereby lP and AP are the length and the cross sectional area of the actuator, and NP is the number of piezoelectric stack layers. The compliance s, the dielectric permittivity ε and the piezoelectric charge constant d are material parameters. Assuming a piece-wise constant input vector u(k) = u(t)|t∈[tk ,tk +T [ for time spans of length T , both linear submodels can be integrated analyti T i i cally. Using Gi = eA T and Hi = 0 eA t Bi dt, equation 1 results in
434
xi (k) = Gi x(k − 1) + Hi u(k − 1)
.
(5)
This equation allows a conclusion from the old state vector at tk−1 to the new vectors at tk = tk−1 + T for each valve state. The appropriate state vector x(k) together wit the output vector y(k) are selected by using the case distinction ⎧ T 1 1 1 ⎪ 0 ⎪ ⎪ x1 (k) x2 (k) x3 (k) 0 ⎪ ⎨ for hN (k) ≤ 0 x(k) =
2 T ⎪ ⎪ x1 (k) x22 (k) x23 (k) x24 (k) x25 (k) ⎪ ⎪ ⎩
elements in the Extended DHM are the operators F Merk and xMerk which are stack memories for past relevant turnaround point coordinates within the (x, F )-plane. Using the operator Γ−1 to denote the numerical solution x = Γ−1 (F, t, θ) of equation 8, the polarization P of the actuator can be given by
(6)
with θ4 = (7)
hN (norm.)
P
U (norm.)
with hN = x24 . The parameterized equations 5 to 7 form a time discrete model for the inspected injector.
0.5
1 t [ms]
10 0 −10 10 0 −10 0
0.5
1 t [ms]
1.5
1 0.5 0 0
0.5
1 t [ms]
0
0.5
1 t [ms]
Table 3. Elements of the state vector x Element xi1
10 0 −10 −20 10 0 −10 −20
3. LUMPED PARAMETERS WITH HYSTERESIS
xi2
Actuator polarization P
xi3
Actuator position hP Actuator velocity h˙ P Valve clearance hN Needle velocity h˙ N T
Using the input vector u = (IP FN ) and output vector y = (UP hN )T again, the matrices for the differential equations in 1 and 2 now are ⎛0 0 0 0 0 0 ⎞ ⎜0 ⎜0 1 A =⎜ ⎜0 ⎝ ⎛0 ⎜0 ⎜0 ⎜ 2 A =⎜ ⎜0 ⎜ ⎝0 ⎛
dx dt
0
+ θ4
describes a friction force F for a displacement x as a nonlinear differential equation. Here, t is the time, θ = (θ1 θ2 θ3 θ4 ) is a parameter vector and sign(·) is the signum function. The very specific
1
B =
∂F
1
m2 ∂P 0 0 0 0
m2
0
0 1 m2
∂F
0 0
⎞
0
P −k −k 2 3 0 k3
⎛
1 ∂P
0 −
∂hP 1 r2 − m2 0 0 −
m3
⎜ ⎟ 2 ⎜∂Q0 P ⎟ ⎜ B = 0 ⎜ 0 ⎠ 0 ⎝ 0 0
∂hP 1 r2 − m2 0 0 0 ∂P
∂hP
0
0
0 0
0
0
⎜∂QP ⎜ 0 ⎝ 0
P −k −k 2 3
∂hP
0
m2 ∂P 0
1 ∂P
0
0
0 1 ∂FP
∂P
0
0 1 ∂FP
0
Based on the Dynamic Hysteresis Model DHM (Dahl and Wilder, 1985), the Extended DHM (Mehlfeldt, n.d.) (8) θ3 Merk Merk −θ4 x−x F −F sign
Physical meaning Electrical actuator charge QP
xi4 xi5 xi6
1.5
The normalized simulation results are compared to measurements in figure 3. This is done for two different injections: In the upper injection, the actuator is charged and discharged in 100µs each, whereas in the lower one, charging and discharging takes 150µs each. On the left hand side, the calculated voltage UP shows errors of up to 11%. On the right hand side, the calculated valve opening hN deviates by up to 27% from the measured values.
θ2
(11)
with the new material properties εP , k σ , k E and sP (Mehlfeldt, n.d.).
1.5
Fig. 3. Lumped masses, no hysteresis: Simulations (solid) and measurements (dashed)
dF = θ1 2− dx
(10)
Since the polarization is given by a differential equation of first order, it is usefull to take P into the state vector x. Thus, x now is composed by the six elements listed in table 3.
1 0.5 0
1.5
hN error [%]
P
U error [%]
0
UP and FP are now given by
UP =
y(k) = C x(k) + D u(k)
1 0.5 0
kσ kE sP .
(9)
lP l 1 1 QP − P P P 2 A εP NP N ε P P kσ AP 1 h F P = AP P P − s lP sP P
for hN (k) > 0
1 0.5 0
1 kE 1 1 Q + h , t, θ NP A P ε P P lP sP P
P = Γ−1
0 0 0 0 0 1 m3
⎞ ⎟ ⎟ ⎟ ⎟ ⎠
⎛∂U
0
0
0
0
0
0
0 0
0 0
0
0
0
0
0
0
0
0
0 k3
1 r3
m3
−
⎞T
⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠
m3
P 0
P ⎜∂Q ⎜∂UP C = ⎜ ∂P ⎜ 0 ⎝ 0 0 0
⎟ ⎟ ⎟ 0⎟ 0⎠ 0
D=0 .
1 0
Contrary to the previous section, the matrices A1 , A2 , B1 and B2 are not constant here. Because the ∂P ∂P partial derivatives ∂h and ∂Q are defined by the P P nonlinear, hard switching differential equation 9, the two sets of differential equations for i=1 and i=2 are also nonlinear and hard switching.
435
Since analytical solutions are not available for this type, the differential equations have to be integrated numerically, e.g. by the forward Euler method. (12) i i i x (k) ≈ T A x(k − 1) x(k − 1) + T B x(k − 1) u(k − 1)
Again, a case distinction selects one single state vector x out of x1 and x2 evaluating hN = x25 . ⎧ T ⎪ x11 (k) x12 (k) x13 (k) x14 (k) 0 0 ⎪ ⎪ ⎪ ⎨ for hN (k) ≤ 0 x(k) =
2 T (13) ⎪ 2 2 2 2 2 ⎪ x (k) x2 (k) x3 (k) x4 (k) x5 (k) x5 (k) ⎪ ⎪ ⎩ 1 for hN (k) > 0
y(k) = C x(k) + D u(k)
(14)
1 0.5 0
hN (norm.)
P
U (norm.)
For figure 4, the time discrete model from the parameterized equations 12 to 14 is simulated and the normalized simulation results are compared to the measurements.
1 0.5 0 0.5
1 t [ms]
10 0 −10 10 0 −10 0
0.5
1 t [ms]
1.5
A calculus of variation of the distributed action results in the describing differential equation for the displacements h (Wauer, 1997). (15) 2 2 A
∂ h ∂ ∂h ε ∂ h − A + Aη ∂t2 ∂z εs−d2 ∂z ∂z∂t
= Af−
∂ d A D ∂z εs−d2
N s ∂h d D− = PU εs−d2 εs−d2 ∂z lP
1 0.5 0 0
0.5
1 t [ms]
1.5
10 0 −10 −20 10 0 −10 −20 0
0.5
1 t [ms]
1.5
Fig. 4. Lumped masses, with hysteresis: Simulations (solid) and measurements (dashed) With this model, the maximum error of the calculated voltage UP is reduced to 9% while the calculated valve opening hN still deviates by up to 25% from the measured values.
4. DISTRIBUTED PARAMETERS WITHOUT HYSTERESIS For the models with distributed parameters, figure 5 shows the simulated arrangement. Tubular Spring
For a mathematical description, all parts of the injector are regarded as one-dimensional homogeneous elements. The piezoelectricity is linearized as before. All parameters of the arrangement are time-invariant functions of the position z: The cross sectional area A(z), the density (z), the piezoelectric charge constant d(z), the electrical permittivity ε(z), the compliance s(z) and the dynamic viscosity η(z).
1 0.5 0
1.5
hN error [%]
P
U error [%]
0
are tied together with a tubular spring. The needle of the injector is pressed against the actuator composite with another spring. As can be seen, some of the individual parts are connected forcefit, but others are connected form-fit only.
Needle
Valve
(16)
Herein, f = f (z) is a distributed force density acting on the continuum. Equation 15 does not only describe the piezoelectric actuator, but with d = 0 also the other pure mechanical parts. If η is chosen within physically reasonable limits, equation 15 is an inhomogeneous hyperbolic partial differential equation of second order with time variant boundary conditions resulting from formfit interconnections. A closed analytical solution for this initial value problem is not available. Therefore, a discretization in space is performed to approximate the partial differential equations by a system of ordinary differential equations. As can be seen in figure 5, some parameters like the cross sectional area may show discontinuities. The arising discretization problems can be solved by using the finite integration (Marklein, 1998): A section-wise averaging ϕ (z) ≈ ϕ (z) = z+ 12 z 1 z z− 12 z ϕ (z) dz with the step size z approximates any function ϕ (z) by a continuous function ϕ (z) which is suitable for discretization. Using the equidistant grid zk = z0 +k ·z and the staggered auxiliary grid zk = z0 + (k + 12 ) · z, the finite integration and subsequent discretization of equations 15 and 16 results in the desired system of differential equations for the discrete displacements h (zk ). (17)
Actuator zP1
Needle Spring zP2
Valve Seat z
Fig. 5. Injector mechanics considered distributed This more detailed perspective locates the actuator between two consecutive metal elements which
As−1 (zk−1 ) Aη(zk ) dh(zk+1 ) d2 h(zk ) A(zk ) − h(zk−1 ) − 2 2 2 dt z dt z
+ −
Aη(zk ) + Aη(zk−1 ) dh(zk ) Aη(zk−1 ) dh(zk−1 ) − 2 z dt z 2 dt
As−1 (zk ) + As−1 (zk−1 ) As−1 (zk ) h(zk+1 ) + h(zk ) z 2 z 2
= A(zk )f (zk ) −
NP Ads−1 (zk ) − Ads−1 (zk−1 ) UP lP z
436
Q=
2 z NP
l2
P
+
NP lP
kP2
A(ε − d2 s−1 )(zk ) UP
(18)
k=kP1 kP2
Ads−1 (zk ) (h(zk+1 ) − h(zk ))
k=kP1
The equations 17 and 18 describe a system of mass-spring elements which are excited by external forces f (zk ) and a position independent actuator voltage UP . As already mentioned, the form-fit connections lead to state dependent discontinuous boundary conditions. Even more, adjacent mechanical parts must be modeled together, while separated parts propagate independently. Hence, with N switching boundary conditions, up to 2N +1 −1 individual differential equation systems with different sizes and dynamics had to be formulated. More feasible is the procedure to first formulate all mathematical models for the individual parts independently. This results in ordinary linear differential equations for all parts which can be integrated analytically for a time step T . If two consecutive parts are in contact after a simulation step, the positions and velocities of the involved discrete masses must be adjusted appropriately, before the next simulation step is performed. In this way, both the mathematical formulation and the simulation process are simplified (Mehlfeldt, n.d.).
1 0.5 0
hN (norm.)
P
U (norm.)
For the simulation the arrangement in figure 5 is discretizised with z = 2mm which results in a mathematical model of order 120 in total. The normalized simulation results of the parameterized equations 17 and 18 are compared to measurements in figure 6.
1 0.5 0 0.5
1 t [ms]
10 0 −10 10 0 −10 0
0.5
1 t [ms]
1.5
5. DISTRIBUTED PARAMETERS WITH HYSTERESIS For the final injector model with distributed parameters and nonlinear piezoelectricity, a review of the previous models helps to find the corresponding mathematical equations. The second model in section 3 is similar to the first model in section 2; basically only the linear piezoelectric equations 3 and 4 are replaced by the nonlinear equations 9 to 11. In the same way, a model with distributed parameters and hysteresis is obtainable by taking the previously described model and substituting analogously the linear piezoelectric relations by the nonlinear ones. This procedure results in the ordinary differential equations A(zk )
+
0
0.5
1 t [ms]
1.5
10 0 −10 −20 10 0 −10 −20
A(sP )−1 (zk−1 ) Aη(zk ) dh(zk+1 ) h(zk−1 ) − 2 2 z dt z
−
Aη(zk−1 ) dh(zk−1 ) A(sP )−1 (zk ) h(zk+1 ) − z 2 z 2 dt
A(sP )−1 (zk ) + A(sP )−1 (zk−1 )
z 2 −
0.5
1 t [ms]
1.5
Fig. 6. Distributed masses, no hysteresis: Simulations (solid) and measurements (dashed) Although the maximum error of the calculated voltage UP still reaches 15%, the calculated valve opening hN now deviates by only less than 14% from the measured values. In other words, the presented distributed approach is mechanically twice as accurate as the lumped one.
h(zk ) = A(zk )f (zk )
(19)
Akσ (sP )−1 (zk−1 ) Akσ (sP )−1 (zk ) )+ ) P (zk P (zk−1 z z
for the displacements h (zk ) together with the group of hysteresis operators (20) NP h(zk+1 )−h(zk ) UP + kE (sP )−1 (zk ) , t, θ lP z
with θ4 = k σ k E (sP )−1 (zk ) for the polarization P and the total electrical actuator charge Q=
0
Aη(zk ) + Aη(zk−1 ) dh(zk ) d2 h(zk ) + 2 2 dt z dt
−
) = Γ−1 P (zk
1 0.5 0
1.5
hN error [%]
P
U error [%]
0
1 0.5 0
The simulations in figure 6 were performed with the program COMPAS which was developed at DaimlerChrysler AG. This versatile simulation tool for piezoelectrically actuated one-dimensional arrangements is configured via plain text files that can be generated automatically by MATLAB scripts.
2 NP z
l2P
kP2
AεP (zk )UP
k=kP1
N z + P lP
kP2
(21)
A(zk )P (zk )
.
k=kP1
Equations 19 to 21 build an implicit system of nonlinear ordinary differential equations. Though a general explicit solution does not exist, the implicit system is always momentary uniquely solvable into an explicit system of nonlinear differential equations (Mehlfeldt, n.d.). Extrapolating that each individual momentary solution is valid for a finite period T , the individual systems of nonlinear differential equations can be integrated numerically over T , e.g. by the forward Euler method. As described in the previous section, the boundary conditions may require additional intermediate state corrections.
437
1 0.5 0
hN (norm.)
P
U (norm.)
For the simulation, the arrangement in figure 5 is discretizised once more with z = 2mm. Again, the simulation is done with COMPAS. The normalized simulation results of the parameterized equations 19 to 21 are compared to measurements in figure 7.
1 0.5 0 0.5
1 t [ms]
10 0 −10 10 0 −10 0
0.5
1 t [ms]
1 0.5 0
1.5
1.5
6. SUMMARY AND CONCLUSION 0
hN error [%]
P
U error [%]
0
1 0.5 0
charging and discharging demonstrates that the form-fit connection between the actuator and the housing is lost for a short time. The following collision stresses the mechanical parts with extremely high pressure peaks. When the injector is opened and closed slowly (graph below), the simulated forces seemingly are much more tolerable for the mechanics.
0.5
1 t [ms]
1.5
1 t [ms]
1.5
10 0 −10 −20 10 0 −10 −20 0
0.5
Fig. 7. Distributed masses, with hysteresis: Simulations (solid) and measurements (dashed) Both the maximum error of the calculated voltage UP is reduced to 9% and the calculated valve opening hN again deviates by only less than 11% from the measured values. In summary, this model is almost twice as accurate as the first one. This model furthermore permits a detailed view inside the mechanical processes within the injector during an injection.
Additional to the existing linear models, this article presents four more realistic mathematical descriptions of the electromechanical behaviour of piezoelectric actuated injectors. Table 4 summarizes the achieved accuracies of the different models. Table 4. Achieved Model Accuracies Electrical Model error (max.) Lumped, no hystersis 11% Lumped, with hystersis 9% Distributed, no hystersis 15% Distributed, with hystersis 9%
Mechanical error (max.) 27% 25% 14% 11%
While the first two models are suited e. g. for model based observer development, the latter ones allow a detailed inspection in the dynamic processes during an injection. Especially the last described model is accurate enough to develop appropriate injector control strategies in simulation.
REFERENCES
Fig. 8. Simulated displacements and forces within the arrangement Therefore, figure 8 presents the two simulations from figure 7 in a different way. The graphs show the mechanical pressure within the needle, the intermediate parts and the actuator (from top to bottom) over time. Positive pressure is dark colored, tension is light colored. The pressures are drawn at the displaced positions, while the displacement is exagerated by a factor of 500. The changes between light and dark color visualize the vibrations within the injector mechanics. The upper graph for the injection with fast
Dahl, P.E. and R. Wilder (1985). Math model of hysteresis in piezo-electric actuators for precision pointing systems. Guidance and Control 57, 61–88. Gottlieb, B., A. Kappel, R. Mock, H. Meixner and B. Fischer (2001). Dosierventil. Patent number DE 199 40 055 C 1. Kasper, R., J. Schr¨oder and A. Wagner (1997). Schnellschaltendes hydraulikventil mit piezoelektrischem stellan¨ trieb. Olhydraulik und Pneumatik 41, 694– 698. Marklein, R. (1998). Ndt related quantitative modeling of coupled piezoelectric and ultrasonic wave phenomena. ECNDT 1998 Conference. Copenhagen, Denmark. Mehlfeldt, D. (n.d.). Modellierung und optimale steuerung piezoelektrisch aktuierter einspritzventile. Ph.D. thesis. Universit¨ at Siegen. Germany. To be published in 2006. Voigt, W. (1966). Lehrbuch der Kristallphysik. Teubner Verlag. Stuttgart. Wauer, J. (1997). Zur modellierung piezoelektrischer wandler mit verteilten parametern. Z. angew. Math. Mech. 77, 365–366.
438
Abstract: In this paper, four mathematical modeling approaches for the electromechanical behaviour of piezoelectrically actuated fuel injectors are discussed. Models with lumped and distributed parameters are investigated together with linear and nonlinear piezoelectric submodels. For all four resulting sets of nonlinear differential equations, numerical integration methods are mentioned. The simulation results of the parameterized models are compared to measurements of a real injection valve. Keywords: Hysteresis, Mathematical models, Numerical simulation, Partial differential equations, Time-varying systems Classification: A (Methods and tools for multi-disciplinary systems development)
1. INTRODUCTION Together with a growing competition in fuel economy aggravated emission restrictions force the automotive engineers to enhance consecutively the thermodynamic processes of internal combustion engines. One of the most recent developments is the gasoline engine with direct injection and stratified combustion. Within this engine concept, the injector is a key component. Figure 1 outlines such a piezoelectrically actuated injector (Gottlieb et al., 2001). The piezoelectric stack actuator with its electrical connections is located in the center of the injector. On the left hand side, it is coupled to a metal element which is in contact to the housing. To the right, the actuator acts against the needle which controls the hydraulic valve at the very right 1
Supported by Conti Temic microelectronic GmbH
Fig. 1. Piezoelectrically actuated injector end. Whenever the actuator is being charged, the induced strain pushes the needle to the right so that the valve is opened. If the actuator is being discharged, a spring pushes the needle back, and the valve closes. Due to the stiff couplings and due to the small inertia of the involved masses, it takes less than 200µs to switch from close to open or back. The common mathematical modeling approach for injectors is an ordinary linear differential equation (Kasper et al., 1997). This approach neither reflects the switching behaviour of the valve nor
433
x˙ i = Ai x + Bi u
the nonlinearity of piezoelectric actuators or the distributed character of the real mechanics. This article introduces and compares four new and improved mathematical modeling approaches. Each mathematical model is parameterized, and the simulation results are compared to measurements. 2. LUMPED PARAMETERS WITHOUT HYSTERESIS At first, the relevant injector mechanics are considered as lumped masses with discrete interactions. Choosing two lumped masses m1 and m2 for both ends of the actuator, together with the consecutive elements and one lumped mass m3 for the valve part of the needle, the injector mechanics become an arrangement as shown in figure 2. m1 k1
FP k2
r1
hN
m2
r3
Fig. 2. Injector mechanics considered as lumped FP is the piezoelectrically induced stress within the actuator. FN is a static force which closes the valve securely when the actuator is discharged. The discrete elements k1 , k2 and k3 represent the stiffness of the housing, the actuator and the needle, respectively. Elements labeled with r stand for velocity proportional dampings. It turned out that the investigated fuel valve has a relatively stiff housing. In this special case, m1 never changes position in the modeled arrangement. Thus, it is omitted in the further examinations. The switching behaviour of the valve requires a distinction of the two cases for the two hydraulic states. Let hN denote the valve clearance. If the valve is closed (hN ≤ 0), the needle is in contact to the valve seat, so that only m2 is free to move. If the valve is opened (hN > 0), the positions of m2 and m3 are variable. Table 1 compares the two resulting simplified arrangements. Table 1. Simplified arrangements Valve state Valve closed hN ≤ 0 (i = 1) Valve opened hN > 0 (i = 2)
Arrangement FP
m2
k2 r2 FP k2 r2
k3 m2
m3 F N k3 r3
Both simplified arrangements can be described separately by sets of ordinary linear differential equations
(1)
i
y = Cx +Du
(2)
with index i ∈ [1, 2] denoting the respective valve state. The input vector u = (IP FN )T is composed by the actuator current IP and the closing force on T the needle FN . The output vector y = (UP hN ) consists of the computed actuator voltage UP and the calculated valve opening hN . As can be seen in table 2, the elements of the state vector x are the actuator charge and the kinetic variables of the lumped masses. Table 2. Elements of the state vector x Element xi1
Physical meaning Electrical actuator charge QP Actuator position hP Actuator velocity h˙
xi2 xi3
P
xi4
m3 F N
k3 r2
i
Valve clearance hN Needle velocity h˙ N
xi5
With these assumptions, the matrices in the equations 1 and 2 become ⎛ 0 0 0 0 0 ⎞ 0 0 1 0 0 r2 ⎜ 1 ∂FP 1 ∂FP − k − k ⎟ 1 A =⎝ − 0 0 ⎠ 2 3 m2 ∂Q m2 ∂h m2 ⎛ 2
A =
0 0
0 0 1 ∂FP
⎜ ⎜ ⎜m2 ⎝
∂QP 0
P
P
1 m2
∂FP ∂hP
1
B
1 0 0 0
0 0
− k2 − k3
⎛1
⎞
0
2
B
−
⎜0 =⎝ 0 0 −
0 0 0 0 1 m3
0 1 r2 m2 0 0
m3
=⎝ 0 0 ⎠ 0 0 0 0
0 0
0 k3
0
⎛
0 0
⎞ ⎟ ⎠
⎛ ⎜ ⎜ C=⎜ ⎝
0 0
0 0 k3
0 0
0 0
0 m2 0 1 k3 r3 − − m3 m3
∂UP ∂QP ∂UP ∂hP 0 0 0
0
⎞ ⎟ ⎟ ⎟ ⎠
⎞T
⎟ ⎟ ⎟ 0⎠ 0
D=0 .
1 0
In this section, a linearized piezoelectricity is considered. Thus, the partial derivatives in the matrices above are to be determined using the linear relations 1 1 s lP d QP − h 2 εs − d2 A NP NP εs − d2 P P d AP 1 ε Q − h FP = NP εs − d2 P εs − d2 lP P
UP =
(3) (4)
between the voltage UP , the charge QP , the force FP and the displacement hP of the actuator (Voigt, 1966). Thereby lP and AP are the length and the cross sectional area of the actuator, and NP is the number of piezoelectric stack layers. The compliance s, the dielectric permittivity ε and the piezoelectric charge constant d are material parameters. Assuming a piece-wise constant input vector u(k) = u(t)|t∈[tk ,tk +T [ for time spans of length T , both linear submodels can be integrated analyti T i i cally. Using Gi = eA T and Hi = 0 eA t Bi dt, equation 1 results in
434
xi (k) = Gi x(k − 1) + Hi u(k − 1)
.
(5)
This equation allows a conclusion from the old state vector at tk−1 to the new vectors at tk = tk−1 + T for each valve state. The appropriate state vector x(k) together wit the output vector y(k) are selected by using the case distinction ⎧ T 1 1 1 ⎪ 0 ⎪ ⎪ x1 (k) x2 (k) x3 (k) 0 ⎪ ⎨ for hN (k) ≤ 0 x(k) =
2 T ⎪ ⎪ x1 (k) x22 (k) x23 (k) x24 (k) x25 (k) ⎪ ⎪ ⎩
elements in the Extended DHM are the operators F Merk and xMerk which are stack memories for past relevant turnaround point coordinates within the (x, F )-plane. Using the operator Γ−1 to denote the numerical solution x = Γ−1 (F, t, θ) of equation 8, the polarization P of the actuator can be given by
(6)
with θ4 = (7)
hN (norm.)
P
U (norm.)
with hN = x24 . The parameterized equations 5 to 7 form a time discrete model for the inspected injector.
0.5
1 t [ms]
10 0 −10 10 0 −10 0
0.5
1 t [ms]
1.5
1 0.5 0 0
0.5
1 t [ms]
0
0.5
1 t [ms]
Table 3. Elements of the state vector x Element xi1
10 0 −10 −20 10 0 −10 −20
3. LUMPED PARAMETERS WITH HYSTERESIS
xi2
Actuator polarization P
xi3
Actuator position hP Actuator velocity h˙ P Valve clearance hN Needle velocity h˙ N T
Using the input vector u = (IP FN ) and output vector y = (UP hN )T again, the matrices for the differential equations in 1 and 2 now are ⎛0 0 0 0 0 0 ⎞ ⎜0 ⎜0 1 A =⎜ ⎜0 ⎝ ⎛0 ⎜0 ⎜0 ⎜ 2 A =⎜ ⎜0 ⎜ ⎝0 ⎛
dx dt
0
+ θ4
describes a friction force F for a displacement x as a nonlinear differential equation. Here, t is the time, θ = (θ1 θ2 θ3 θ4 ) is a parameter vector and sign(·) is the signum function. The very specific
1
B =
∂F
1
m2 ∂P 0 0 0 0
m2
0
0 1 m2
∂F
0 0
⎞
0
P −k −k 2 3 0 k3
⎛
1 ∂P
0 −
∂hP 1 r2 − m2 0 0 −
m3
⎜ ⎟ 2 ⎜∂Q0 P ⎟ ⎜ B = 0 ⎜ 0 ⎠ 0 ⎝ 0 0
∂hP 1 r2 − m2 0 0 0 ∂P
∂hP
0
0
0 0
0
0
⎜∂QP ⎜ 0 ⎝ 0
P −k −k 2 3
∂hP
0
m2 ∂P 0
1 ∂P
0
0
0 1 ∂FP
∂P
0
0 1 ∂FP
0
Based on the Dynamic Hysteresis Model DHM (Dahl and Wilder, 1985), the Extended DHM (Mehlfeldt, n.d.) (8) θ3 Merk Merk −θ4 x−x F −F sign
Physical meaning Electrical actuator charge QP
xi4 xi5 xi6
1.5
The normalized simulation results are compared to measurements in figure 3. This is done for two different injections: In the upper injection, the actuator is charged and discharged in 100µs each, whereas in the lower one, charging and discharging takes 150µs each. On the left hand side, the calculated voltage UP shows errors of up to 11%. On the right hand side, the calculated valve opening hN deviates by up to 27% from the measured values.
θ2
(11)
with the new material properties εP , k σ , k E and sP (Mehlfeldt, n.d.).
1.5
Fig. 3. Lumped masses, no hysteresis: Simulations (solid) and measurements (dashed)
dF = θ1 2− dx
(10)
Since the polarization is given by a differential equation of first order, it is usefull to take P into the state vector x. Thus, x now is composed by the six elements listed in table 3.
1 0.5 0
1.5
hN error [%]
P
U error [%]
0
UP and FP are now given by
UP =
y(k) = C x(k) + D u(k)
1 0.5 0
kσ kE sP .
(9)
lP l 1 1 QP − P P P 2 A εP NP N ε P P kσ AP 1 h F P = AP P P − s lP sP P
for hN (k) > 0
1 0.5 0
1 kE 1 1 Q + h , t, θ NP A P ε P P lP sP P
P = Γ−1
0 0 0 0 0 1 m3
⎞ ⎟ ⎟ ⎟ ⎟ ⎠
⎛∂U
0
0
0
0
0
0
0 0
0 0
0
0
0
0
0
0
0
0
0 k3
1 r3
m3
−
⎞T
⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠
m3
P 0
P ⎜∂Q ⎜∂UP C = ⎜ ∂P ⎜ 0 ⎝ 0 0 0
⎟ ⎟ ⎟ 0⎟ 0⎠ 0
D=0 .
1 0
Contrary to the previous section, the matrices A1 , A2 , B1 and B2 are not constant here. Because the ∂P ∂P partial derivatives ∂h and ∂Q are defined by the P P nonlinear, hard switching differential equation 9, the two sets of differential equations for i=1 and i=2 are also nonlinear and hard switching.
435
Since analytical solutions are not available for this type, the differential equations have to be integrated numerically, e.g. by the forward Euler method. (12) i i i x (k) ≈ T A x(k − 1) x(k − 1) + T B x(k − 1) u(k − 1)
Again, a case distinction selects one single state vector x out of x1 and x2 evaluating hN = x25 . ⎧ T ⎪ x11 (k) x12 (k) x13 (k) x14 (k) 0 0 ⎪ ⎪ ⎪ ⎨ for hN (k) ≤ 0 x(k) =
2 T (13) ⎪ 2 2 2 2 2 ⎪ x (k) x2 (k) x3 (k) x4 (k) x5 (k) x5 (k) ⎪ ⎪ ⎩ 1 for hN (k) > 0
y(k) = C x(k) + D u(k)
(14)
1 0.5 0
hN (norm.)
P
U (norm.)
For figure 4, the time discrete model from the parameterized equations 12 to 14 is simulated and the normalized simulation results are compared to the measurements.
1 0.5 0 0.5
1 t [ms]
10 0 −10 10 0 −10 0
0.5
1 t [ms]
1.5
A calculus of variation of the distributed action results in the describing differential equation for the displacements h (Wauer, 1997). (15) 2 2 A
∂ h ∂ ∂h ε ∂ h − A + Aη ∂t2 ∂z εs−d2 ∂z ∂z∂t
= Af−
∂ d A D ∂z εs−d2
N s ∂h d D− = PU εs−d2 εs−d2 ∂z lP
1 0.5 0 0
0.5
1 t [ms]
1.5
10 0 −10 −20 10 0 −10 −20 0
0.5
1 t [ms]
1.5
Fig. 4. Lumped masses, with hysteresis: Simulations (solid) and measurements (dashed) With this model, the maximum error of the calculated voltage UP is reduced to 9% while the calculated valve opening hN still deviates by up to 25% from the measured values.
4. DISTRIBUTED PARAMETERS WITHOUT HYSTERESIS For the models with distributed parameters, figure 5 shows the simulated arrangement. Tubular Spring
For a mathematical description, all parts of the injector are regarded as one-dimensional homogeneous elements. The piezoelectricity is linearized as before. All parameters of the arrangement are time-invariant functions of the position z: The cross sectional area A(z), the density (z), the piezoelectric charge constant d(z), the electrical permittivity ε(z), the compliance s(z) and the dynamic viscosity η(z).
1 0.5 0
1.5
hN error [%]
P
U error [%]
0
are tied together with a tubular spring. The needle of the injector is pressed against the actuator composite with another spring. As can be seen, some of the individual parts are connected forcefit, but others are connected form-fit only.
Needle
Valve
(16)
Herein, f = f (z) is a distributed force density acting on the continuum. Equation 15 does not only describe the piezoelectric actuator, but with d = 0 also the other pure mechanical parts. If η is chosen within physically reasonable limits, equation 15 is an inhomogeneous hyperbolic partial differential equation of second order with time variant boundary conditions resulting from formfit interconnections. A closed analytical solution for this initial value problem is not available. Therefore, a discretization in space is performed to approximate the partial differential equations by a system of ordinary differential equations. As can be seen in figure 5, some parameters like the cross sectional area may show discontinuities. The arising discretization problems can be solved by using the finite integration (Marklein, 1998): A section-wise averaging ϕ (z) ≈ ϕ (z) = z+ 12 z 1 z z− 12 z ϕ (z) dz with the step size z approximates any function ϕ (z) by a continuous function ϕ (z) which is suitable for discretization. Using the equidistant grid zk = z0 +k ·z and the staggered auxiliary grid zk = z0 + (k + 12 ) · z, the finite integration and subsequent discretization of equations 15 and 16 results in the desired system of differential equations for the discrete displacements h (zk ). (17)
Actuator zP1
Needle Spring zP2
Valve Seat z
Fig. 5. Injector mechanics considered distributed This more detailed perspective locates the actuator between two consecutive metal elements which
As−1 (zk−1 ) Aη(zk ) dh(zk+1 ) d2 h(zk ) A(zk ) − h(zk−1 ) − 2 2 2 dt z dt z
+ −
Aη(zk ) + Aη(zk−1 ) dh(zk ) Aη(zk−1 ) dh(zk−1 ) − 2 z dt z 2 dt
As−1 (zk ) + As−1 (zk−1 ) As−1 (zk ) h(zk+1 ) + h(zk ) z 2 z 2
= A(zk )f (zk ) −
NP Ads−1 (zk ) − Ads−1 (zk−1 ) UP lP z
436
Q=
2 z NP
l2
P
+
NP lP
kP2
A(ε − d2 s−1 )(zk ) UP
(18)
k=kP1 kP2
Ads−1 (zk ) (h(zk+1 ) − h(zk ))
k=kP1
The equations 17 and 18 describe a system of mass-spring elements which are excited by external forces f (zk ) and a position independent actuator voltage UP . As already mentioned, the form-fit connections lead to state dependent discontinuous boundary conditions. Even more, adjacent mechanical parts must be modeled together, while separated parts propagate independently. Hence, with N switching boundary conditions, up to 2N +1 −1 individual differential equation systems with different sizes and dynamics had to be formulated. More feasible is the procedure to first formulate all mathematical models for the individual parts independently. This results in ordinary linear differential equations for all parts which can be integrated analytically for a time step T . If two consecutive parts are in contact after a simulation step, the positions and velocities of the involved discrete masses must be adjusted appropriately, before the next simulation step is performed. In this way, both the mathematical formulation and the simulation process are simplified (Mehlfeldt, n.d.).
1 0.5 0
hN (norm.)
P
U (norm.)
For the simulation the arrangement in figure 5 is discretizised with z = 2mm which results in a mathematical model of order 120 in total. The normalized simulation results of the parameterized equations 17 and 18 are compared to measurements in figure 6.
1 0.5 0 0.5
1 t [ms]
10 0 −10 10 0 −10 0
0.5
1 t [ms]
1.5
5. DISTRIBUTED PARAMETERS WITH HYSTERESIS For the final injector model with distributed parameters and nonlinear piezoelectricity, a review of the previous models helps to find the corresponding mathematical equations. The second model in section 3 is similar to the first model in section 2; basically only the linear piezoelectric equations 3 and 4 are replaced by the nonlinear equations 9 to 11. In the same way, a model with distributed parameters and hysteresis is obtainable by taking the previously described model and substituting analogously the linear piezoelectric relations by the nonlinear ones. This procedure results in the ordinary differential equations A(zk )
+
0
0.5
1 t [ms]
1.5
10 0 −10 −20 10 0 −10 −20
A(sP )−1 (zk−1 ) Aη(zk ) dh(zk+1 ) h(zk−1 ) − 2 2 z dt z
−
Aη(zk−1 ) dh(zk−1 ) A(sP )−1 (zk ) h(zk+1 ) − z 2 z 2 dt
A(sP )−1 (zk ) + A(sP )−1 (zk−1 )
z 2 −
0.5
1 t [ms]
1.5
Fig. 6. Distributed masses, no hysteresis: Simulations (solid) and measurements (dashed) Although the maximum error of the calculated voltage UP still reaches 15%, the calculated valve opening hN now deviates by only less than 14% from the measured values. In other words, the presented distributed approach is mechanically twice as accurate as the lumped one.
h(zk ) = A(zk )f (zk )
(19)
Akσ (sP )−1 (zk−1 ) Akσ (sP )−1 (zk ) )+ ) P (zk P (zk−1 z z
for the displacements h (zk ) together with the group of hysteresis operators (20) NP h(zk+1 )−h(zk ) UP + kE (sP )−1 (zk ) , t, θ lP z
with θ4 = k σ k E (sP )−1 (zk ) for the polarization P and the total electrical actuator charge Q=
0
Aη(zk ) + Aη(zk−1 ) dh(zk ) d2 h(zk ) + 2 2 dt z dt
−
) = Γ−1 P (zk
1 0.5 0
1.5
hN error [%]
P
U error [%]
0
1 0.5 0
The simulations in figure 6 were performed with the program COMPAS which was developed at DaimlerChrysler AG. This versatile simulation tool for piezoelectrically actuated one-dimensional arrangements is configured via plain text files that can be generated automatically by MATLAB scripts.
2 NP z
l2P
kP2
AεP (zk )UP
k=kP1
N z + P lP
kP2
(21)
A(zk )P (zk )
.
k=kP1
Equations 19 to 21 build an implicit system of nonlinear ordinary differential equations. Though a general explicit solution does not exist, the implicit system is always momentary uniquely solvable into an explicit system of nonlinear differential equations (Mehlfeldt, n.d.). Extrapolating that each individual momentary solution is valid for a finite period T , the individual systems of nonlinear differential equations can be integrated numerically over T , e.g. by the forward Euler method. As described in the previous section, the boundary conditions may require additional intermediate state corrections.
437
1 0.5 0
hN (norm.)
P
U (norm.)
For the simulation, the arrangement in figure 5 is discretizised once more with z = 2mm. Again, the simulation is done with COMPAS. The normalized simulation results of the parameterized equations 19 to 21 are compared to measurements in figure 7.
1 0.5 0 0.5
1 t [ms]
10 0 −10 10 0 −10 0
0.5
1 t [ms]
1 0.5 0
1.5
1.5
6. SUMMARY AND CONCLUSION 0
hN error [%]
P
U error [%]
0
1 0.5 0
charging and discharging demonstrates that the form-fit connection between the actuator and the housing is lost for a short time. The following collision stresses the mechanical parts with extremely high pressure peaks. When the injector is opened and closed slowly (graph below), the simulated forces seemingly are much more tolerable for the mechanics.
0.5
1 t [ms]
1.5
1 t [ms]
1.5
10 0 −10 −20 10 0 −10 −20 0
0.5
Fig. 7. Distributed masses, with hysteresis: Simulations (solid) and measurements (dashed) Both the maximum error of the calculated voltage UP is reduced to 9% and the calculated valve opening hN again deviates by only less than 11% from the measured values. In summary, this model is almost twice as accurate as the first one. This model furthermore permits a detailed view inside the mechanical processes within the injector during an injection.
Additional to the existing linear models, this article presents four more realistic mathematical descriptions of the electromechanical behaviour of piezoelectric actuated injectors. Table 4 summarizes the achieved accuracies of the different models. Table 4. Achieved Model Accuracies Electrical Model error (max.) Lumped, no hystersis 11% Lumped, with hystersis 9% Distributed, no hystersis 15% Distributed, with hystersis 9%
Mechanical error (max.) 27% 25% 14% 11%
While the first two models are suited e. g. for model based observer development, the latter ones allow a detailed inspection in the dynamic processes during an injection. Especially the last described model is accurate enough to develop appropriate injector control strategies in simulation.
REFERENCES
Fig. 8. Simulated displacements and forces within the arrangement Therefore, figure 8 presents the two simulations from figure 7 in a different way. The graphs show the mechanical pressure within the needle, the intermediate parts and the actuator (from top to bottom) over time. Positive pressure is dark colored, tension is light colored. The pressures are drawn at the displaced positions, while the displacement is exagerated by a factor of 500. The changes between light and dark color visualize the vibrations within the injector mechanics. The upper graph for the injection with fast
Dahl, P.E. and R. Wilder (1985). Math model of hysteresis in piezo-electric actuators for precision pointing systems. Guidance and Control 57, 61–88. Gottlieb, B., A. Kappel, R. Mock, H. Meixner and B. Fischer (2001). Dosierventil. Patent number DE 199 40 055 C 1. Kasper, R., J. Schr¨oder and A. Wagner (1997). Schnellschaltendes hydraulikventil mit piezoelektrischem stellan¨ trieb. Olhydraulik und Pneumatik 41, 694– 698. Marklein, R. (1998). Ndt related quantitative modeling of coupled piezoelectric and ultrasonic wave phenomena. ECNDT 1998 Conference. Copenhagen, Denmark. Mehlfeldt, D. (n.d.). Modellierung und optimale steuerung piezoelektrisch aktuierter einspritzventile. Ph.D. thesis. Universit¨ at Siegen. Germany. To be published in 2006. Voigt, W. (1966). Lehrbuch der Kristallphysik. Teubner Verlag. Stuttgart. Wauer, J. (1997). Zur modellierung piezoelektrischer wandler mit verteilten parametern. Z. angew. Math. Mech. 77, 365–366.
438