Authors: Shigong Fan, Yixi Cai, Yunxi Shi, Zongzhen Zhang
Categories: Article, WOA, deconvolution, impulse function, piston slap force reconstruction, regulation method, variational mode decomposition
Source: Sensors (Basel, Switzerland)
Doi: 10.3390/s24123833
Authors: Shigong Fan, Yixi Cai, Yunxi Shi, Zongzhen Zhang
In a diesel engine, piston slap commonly occurs concurrently with fuel combustion and serves as the main source of excitation. Although combustion pressure can be measured using sensors, determining the slap force is difficult without conducting tests. In this study, we propose a method to identify the slap force of the piston to solve this difficult problem. The traditional VMD algorithm easily receives noise interference, which affects the value of parameter combination [k, α] and thus affects the extraction accuracy of the algorithm. First, we obtain the transfer function between the incentive and vibration response through percussion tests. Secondly, a variational modal decomposition method based on whale algorithm optimization is used to separate the slap response from the surface acceleration of the block. Finally, we calculated the slap force using the deconvolution method. Deconvolution is a typical inverse problem of mathematics, often prone to ill-conditioning, and the singular value decomposition and regularization method is used to overcome this flaw and improve accuracy. The proposed method provides an important means to evaluate the angular distribution of the slap force, identify the shock positions on the piston liner, and determine the peak value of the waveform which helps us analyze the vibration characteristics of the piston and optimize the structural design of the engine.
The diesel engine offers a comfortable driving experience, particularly for trucks and agricultural machinery, because of its high torque at low speeds. Manufacturers consider the acoustic optimization of diesel engines as the most critical task, in addition to reducing fuel consumption and exhaust gas emissions. The major sources of noise in diesel engines include combustion noise, mechanical noise, and aerodynamic noise. One of the mechanisms of the generated noise in a reciprocating piston engine is the motion of the piston because its clearance is perpendicular to its normal running direction. The gas and inertial forces, along with the connecting rod angle and necessary installation clearance, inevitably result in the transverse motion of the piston. As the piston and the crank mechanism reciprocate through intake, compression, expansion, and exhaust strokes, a connecting rod reaction force (Fa) acts on the piston with alternating direction. This reaction force contains a lateral component (Fs), which causes the piston to move across the clearance and impact the cylinder wall. When the piston contacts the cylinder wall, it generates impacts that excite structure-borne noise within the engine structure.
The generation of slap noise depends on various factors, including the piston–liner gap, type of lubricant used, piston pin offsetting [1], and piston geometry [2]. The analysis of lateral and rotary motion of the piston within the gap between the cylinder liner and piston has been conducted. A model that can predict the forces and response of the engine block resulting from slap has been discussed. Parameters, such as mass, spring, and damping constant, have been predicted using a vibrational mobility model [3]. Several methods have been employed to estimate the impact forces. The finite element method (FEM) has been used to analyze the impact force [4,5]. Furthermore, a two-dimensional lumped parameter model has been proposed to calculate the impact force [6]. Meanwhile, the numerical simulation results are compared with measured vibration velocity level at engine block surfaces to verify the proposed model. As shown in Figure 1, the piston moves in three directions, namely vertical, transverse, and rotating directions. The transverse and rotating directions are called the secondary movement of the piston, which has a great impact on the impact of the piston.
In engineering practice, obtaining the actual slap force through experimental testing is a challenging task. Thus, most methods rely on numerical simulations to estimate the slap, which is extremely valuable in recovering the piston slap force from external structural vibration measurements on an engine block. A good method for identifying the cylinder pressure is applying an inverse filter of a specific structure to the vibration signal in an IC engine [7,8,9,10], which is based on angular sampling and cyclic signal processing, shedding light on the condition monitoring of internal combustion engines.
In this study, we conducted three-part tasks to recover the piston slap. In Section 2, the vibration signal produced by the slap force is separated from the total response using a WOA-based VMD method. Then, the impulse response function is acquired between the slap and the structure vibration of the cylinder block through hammer testing. Finally, a deconvolution method was used to reconstruct the slap force.
Variational mode decomposition is a time-frequency analysis method for processing nonlinear and nonstationary signals, which overcomes the disadvantage of frequency aliasing of empirical mode decomposition [11].
In order to obtain the most favorable degradation results, we use the whale optimization algorithm (WOA) [12] to optimize parameter combinations and reduce the significant impact of parameter selection on the results of variational mode decomposition (VMD). The assumption of this method is that each intrinsic mode function (IMF) has a finite bandwidth and different center frequencies. The conversion process solves the problem of variation in order to minimize the total estimated bandwidth of IMF [13]. Each IMF is demodulated to its corresponding baseband, and each IMF and its corresponding center frequency is derived.
Suppose each component has a limited bandwidth centered around a certain madness that is constantly updated during the separation process. Therefore, the original signal ft can be divided into modal components ukt with a central frequency ωkt, where k represents the standard quota [13,14,15]. The modal degradation algorithm of fractional order constructs and solves the variation multiple function. Equations (1) and (2) can be used to describe the relationship between the impact forces on the inner wall and the surface vibration of the object. (1)at=ast+act (2)ast=ht∗fst
(3)δt+jπt∗ukt where δt is an impact function.
(4) δt+jπt∗ukte−jωkt
(5)minuk,ωk∑k=1K∂tδt+jπt∗ukte−jωkt22s.t.∑k=1Kuk=f where uk=u1,…,uk is the IMF components, and ωk=ω1,…,ωk is the central frequency of each component.
(6) Luk,ωk,λ=α∑k∂tδ(t)+jπt∗uk(t)e−jωkt22+f(t)−∑kuk(t)22+λt,f(t)−∑kuk(t)
The alternating direction method of the multiplier is used for iterative sub optimization to handle sub signals and center frequencies. Subsequently, Equation (6) is transformed into the frequency domain using the Parseval/Plancherel Fourier Isometric Transform:(7)H=h1h2h3⋯hk⋯hmh1h2⋯hk−1⋯hm−1h1⋯hk−2⋯hm−2⋱⋮⋮⋮0h1⋯hm−k0hmhm−1hm⋮⋮⋱hm−k+1hm−k+2⋯hm
The first variable of positive frequency is eliminated, and two optimization problems can be solved as (8)u^kn+1(ω)=f^(ω)−∑i≠ku^i(ω)+λ^(ω)21+2α(ω−ωk)2
The same solving process is used, and the center frequency solving problem is transformed into a frequency domain. (9)ωkn+1=argminωk∫0∞(ω−ωk)2u^k(ω)2dω
Deriving the updated formula for center frequency from Equation (10):(10)ωkn+1=∫0∞ωu^k(ω)2dω∫0∞u^k(ω)2dω
The iterative process continues until the convergence criterion is met, which is defined as the error between the modal responses obtained in two consecutive (11)∑k=1Ku^kn+1−u^kn22u^kn22<ε where ε is the selected limit of convergence criteria.
A new heuristic optimization algorithm called the whale optimization algorithm (WOA) was proposed because Mirjalili et al. were inspired by the social behavior of humpback whales [16]. Since its launch, WOA has gained tremendous appeal in the academic community and has been widely applied in various fields. WOA can be divided into three different stages in mathematics.
Prey first identify the prey. After successfully identifying the prey’s position, use the following formula to mathematically represent (12)D=C.X∗t−Xt (13)Xt+1=X∗t−A.D
X represents the position vector of the whale, X∗ represents the optimal position of the captured target, t represents the current iteration, and A and C represent the coefficient vectors. The iterative update of whale position to obtain the optimal solution is represented by the following (14)A=2a−a (15)C=2⋅r
Bubble net attack “Bubble attack” is a special catching method unique to humpback whales. The spiral updating method of the tight encirclement mechanism and position can be used to analyze the bubble web of cetaceans. Assuming that these two behaviors have an equal probability of 50% for whale position updates, the digital model is as (16)X=t+1=X∗t−A.D if 0≤p<0.5D′⋅ebl⋅cos2πl+X∗t if 1≥p≥0.5 (17)D=C.Xrandomt−Xt (18)Xt+1=Xrandomt−A⋅D
The steps of the suggested technique are depicted in Figure 2.
The discrete time series was denoted as
Impact force fsnn=0,1,⋯,N−1.
Acceleration asnn=0,1,⋯,N−1.
Response function hnn=0,1,⋯,N−1.
The transfer function can be calculated by Equation (19). (19)H(k)=A(k)F(k) (20)A(k)=FFT(as(n)) (21)h(n)=IFFT(H(k)) where FFT and IFFT are the fast Fourier Transform and the inverse FFT.
Research on load identification is very valuable in engineering. Load identification is a typical inverse problem with the nonuniqueness of the solution. In the two previous chapters, the response of the piston slap and the transfer function between the input and output were determined. In this section, we can obtain the slap force by using the deconvolution method of Equation (22).
The vibration discrete time series, response function series, and slap force series were denoted as as=as1, as2, ⋯, asll∈Z h=h1, h2, ⋯, hmm∈Z fs=fs1, fs2, ⋯, fskk∈Z
The discrete convolution equation was obtained using Equation (22).
In Equation (15), the slap force series fs(n) is an unknown series needed to be estimated. According to the convolution algorithm, we can change Equation (22) into a matrix multiplication as Equation (23).
(22)asn=fsn∗hn,n=1,2,⋯,l
In Equation (15), the slap force series fsn is an unknown series needed to be estimated. According to the convolution algorithm, we can change Equation (22) into a matrix multiplication as Equation (23). (23)as=fs⋅H where H is an N by M matrix. No inverse of the matrix H exists. However, a pseudo-inverse matrix denoted as H can be obtained to address this problem. Thus, the slap force can be acquired using the following equation. (24)fs=as·H+
This way, the convolution operation is converted into solving the pseudo-inverse matrix.
The pseudo-inverse matrix in Equation (24), often called the discrete ill-posed problem, was solved. The matrix H may be ill-conditioned or singular and may yield a nonunique solution that can be judged by the condition number of the matrix H (the ratio of the largest singular value to the smallest value). Tikhonov regularization is the most commonly used method for regularizing ill-posed problems [17,18]. The standard approach to solving an underdetermined system of linear equations, such as Equation (23), is known as linear least squares. It seeks to minimize the residual regularization. (25)fs·H−as2 where ||·|| is the Euclidean norm. Tikhonov regularization determines a useful approximation of f^s to provide preference to a particular solution with desirable properties by using a penalized least-squares problem of the form. (26)fs·H−as2+βfs2
The scalar β>0 is the regularization parameter. The key to this approach is making an appropriate choice of regulation parameters β to find a compromise between the data fitting effect and minimum norm solution. The generalized cross-validation (GCV) method is a very popular and successful error-free method for choosing the regularization parameter [19,20]. It does not require knowledge of noise variance. The SVD of H is presented by Equation (10). (27)H=VΛWT where V=V1,V2,⋯,Vk∈Rk×k and W=W1,W2,⋯,Wl∈Rl×l are the orthogonal matrices that satisfy the important relations WWT=I and VVT=I (WT is the Hermite transpose), and Λ is a (always rectangular) diagonal matrix whose diagonal entries σj are the singular values. Equation (27) is substituted into Equation (24), which is rewritten as (28)fs=asWΛ+VT=as∑i=1MP1σiWiViT (29)Λ=σ100⋯00⋯0σ20⋯00⋯00σ3⋯00⋯⋮⋮⋮⋱⋮⋮⋮000⋯σk0⋯000⋮0k×l∈Rk×l
We can obtain the regulation solution using the GCV method. (30)fs=asWΛ+VT=as∑i=1MPσiσi2+βWiViT
The GCV function was denoted as (31)GCVβ=1/kI−B(β)P21/kTr(I−B(β))2 where Bβ=ψψTψ+βIψT and Tr is the trace of the matrix. We selected the best regularization parameter when the GCVβ reached the maximum value.
In our study, the measurements were carried out on an in-line four-stroke four-cylinder diesel engine. The parameters are shown in Table 1. Minimizing the effect of additive noise is crucial. The best way to minimize the effect is to carefully choose the location of the vibration transducer. The cylinder block, which mainly concerns the piston slap, exhibits a satisfactory signal-to-noise ratio. Transducers were placed on the thrust side, and transversal acceleration was considered because it represents the most significant direction for studying the piston slap phenomenon [21,22].
One accelerometer sensor was placed on the thrust side surface of the first cylinder block, which is at the centerline of the cylinder and near the top dead center. Figure 3 shows the location of the accelerometer. The first cylinder is equipped with a pressure sensor in the preheating plug hole and an encoder installed at the end of the crankshaft to obtain the angular position of the crankshaft.
The phenomenon of piston slap occurs when the piston collides with the internal wall of the cylinder. Burning and piston slap occur simultaneously in the process of gas combustion in a cylinder. At this time, the valves corresponding to the cylinder being tested are closed. Thus, its vibration is neglected. Moreover, the gears and the engine accessory system are far from the position. The vibration corresponding to it is also neglected. Therefore, the highest SNR signals obtained at this location consist of only two signals produced by combustion and slap.
All signals were recorded with a sampling frequency of 32 kHz. The discrete time series is denoted as an and pn, n=1,2,3…,N, representing the vibration and the combustion pressure discrete signal series, respectively. N is the number of total sample points. Figure 4 shows a recorder in the case of 1800 rpm and 160 Nm. The signals perform strong transient characteristics. One cycle of the signals, which was synchronous to the rotation speed of the crankshaft, was truncated in the angular window [−20, 60] of the crankshaft, as shown in Figure 4.
A truncated signal a(n) which occurs within the range of −20 °CA–80 °CA is decomposed into four IMF components, and the corresponding frequency characteristics of it are shown in the same figure. The central frequencies of IMFs are 1000 Hz, 3500 Hz, 4500 Hz, and 8000 Hz. The spectrum of combustion pressure and velocity of pressure rise (VPR) are shown in Figure 5 and Figure 6. Obviously, the energy of them is mainly distributed in the range of 0–1000 Hz.
Two articles [22,23] have studied the relationship between combustion pressure and the vibration of cylinder blocks (Figure 7 and Figure 8).
Pressure has most of its energy concentrated at a low frequency, below a few hundred hertz, and its vibration is found to have little energy there due to the high rigidity of the engine block.
The signals of block vibration and cylinder pressure were analyzed using a CWT [23], and the results show that the maximum heat release rate has a strong correlation with the magnitude of the vibrations. A specific bandwidth, vibration signals of 0.3~1.5 kHz, was affected by variation in the heat release rate. The vibrations excited by combustion lasted over 50 °CA.
Based on the analysis of these data, we can draw the following
(1) The energy distribution of combustion excitation is in the range of 0–2000 Hz, and the first component, IMF1, is produced by combustion pressure and piston slap.
(2) The spectrum of the other components, IMF2, IMF3, and IMF4, is higher than the excitation of combustion, so they are produced only by the piston slap.
Due to the fact that the energy of IMF1 is a small part of the total energy of the vibration, this component is ignored, and the other components are selected to construct the response of the piston slap, as shown in Figure 9. Using this method, the characteristics of slap force response are analyzed in detail under different loads; for example, Figure 8 shows the slap force response on four different kinds loads. Firstly, at the same speed and the same torque, the large injection advance angle will cause an increase in the velocity of pressure rise and lead to the increase in slap force response, shown in Figure 8a,b. Secondly, at the same speed and the same injection advance angle, the vibration of slap force increases with increasing torque, shown in Figure 8a,c. If we compare Figure 8b,c, we can reach the conclusion that a lower torque with a large injection advance angle will cause a larger vibration of piston slap. Thirdly, at the same torque and the same injection advance angle, the vibration of slap force increases with increasing speed. From the above analysis, it can be seen that the main factors affecting the vibration of slap are torque, speed, and velocity of pressure rate.
The method for reconstructing the slap force strongly depends on the nature of the input–output relationship between the slap force and its resulting vibration on the cylinder block. For simplicity, a linear relationship is assumed in the vibration transfer path in the engine block. A hammer test is the most convenient and effective way of obtaining the response function of the block.
One of the most important things in the test is selecting the impact position on the cylinder. Piston slaps always occur in the power stroke and on the thrust side of the cylinder bore because the maximum pressure occurs in this process with the rapid rise in pressure (Figure 10 and Figure 11).
According to the force analysis of the piston in the cylinder liner using the formula Fs=Fgtanθ (ignoring the friction between the piston and cylinder bore), the maximum slap force must occur in this period. Figure 12, Figure 13, Figure 14, Figure 15 and Figure 16 display the vibration response distribution from −20 °CA to 60 °CA. Thus, three positions were selected on the thrust side. The corresponding angles of crank are 15 °CA, 35 °CA, and 55 °CA. The distance between the hammer impact location and crank axis can be determined using formula Rcosθ+l2−Rsinθ2, as shown in Figure 8. The position of the accelerometer is shown in Figure 4. The vibration frequency response is up to 10,000 Hz. The useful frequency range of the transducer and the hammer must reach the limit. A test was conducted using a Siemens LMS Test-lab system. The selected accelerometer and hammer are Kistler products. All signals were sampled at 32 kHz.
Figure 10 shows the frequency response function between the input and output. The vibration energy is concentrated above 3000 Hz. The impulse response function was calculated using Equation (6) and is depicted in Figure 11. The vibration decays rapidly to zero, indicating that a large internal damping occurs in the cylinder body.
Using the deconvolution method, the slap forces were obtained at six different engine load (1) 1200 rpm/10 Nm and 120 Nm, (2) 1800 rpm/10 Nm and 160 Nm, and (3) 2400 rpm/10 Nm and 160 Nm, as shown in Figure 12, Figure 13, Figure 14, Figure 15 and Figure 16. All of these engine load conditions have the same injection advance angle of 11 °CA.
The calculation results reveal that, under the same load, the slap force increases with the rotational speed. Similarly, at the same speed, the slap force increases with the engine torque, because the greater torque requires a higher gas pressure. The relationship between the cylinder pressure and slap force can be expressed by the formula Fs=Fgtanθ, so Fs increases with Fg (here, ignoring the friction between the piston and cylinder bore).
In all the obtained results, the piston slaps occur within the range of 0 °CA–30 °CA after the top dead center of the power stroke. The slap force waveforms consistently exhibit two peaks, indicated by the dotted circles in Figure 12, which are attributed to the fact that, within this range, the combustion of the combustible gases is the most intense. The piston slap did not happen at the moment of maximum combustion pressure but at the maximum velocity of pressure rise (VPR). Meanwhile, the greater injection advance angle is always accompanied by greater slap force. The delay period exerts a great influence in the diesel engine combustion phenomenon. The longer the delay, the more rapidly and higher the temperature rises because more fuel will be present in the cylinder before the rate of burning comes under control, which causes rougher running and stronger diesel knock.
A novel method is presented for the estimation of piston slap force in this article. This method can optimize the Κ,α parameter combination of the traditional VMD method using WOA, which improves the intelligence of the algorithm and simplifies the selection process. The vibration signal on the block surface can be decomposed into multiple IMF components by the WOA-VMD method. The slap force response can be obtained using the correlation coefficient between IMF1 and cylinder pressure.
The transfer function between the incentive and vibration response points is obtained using the percussion test. Subsequently, the slap force was calculated using the deconvolution method. In this process, converting the convolution operation into solving the pseudo-inverse transfer matrix is the key step. Sometimes, the transfer matrix is ill-conditioned. To address this issue, Tikhonov regularization is commonly employed to overcome this limitation and improve the accuracy of the calculations.
This method provides a valuable means of evaluating the angular distribution of the slap force, identifying the shock positions on the piston liner, and determining the peak value of the waveform. All of this work can help us to study the vibration characteristics of pistons and enhance the structural design of engines. This method can be extended to various engine configurations and types, including gasoline and hybrid engines, to validate its versatility and robustness. The findings can be incorporated into predictive maintenance frameworks to preemptively address potential issues in engine components, thereby reducing downtime and maintenance costs.