Predicting rotor dynamic coefficient of liquid annular seal
Funded by Zhejiang Province college innovation program(Grant No: T200918)
Rotor dynamic coefficient of liquid annular seal is important for analyzing its influence on the multistage pump. This coefficient usually includes stiffness and damping coefficient. Referring to the method calculating coefficient of bearing, the seal coefficient can be predicted through CFD analysis.
Predicting stiffness coefficient
The force exerted on rotor surface in a certain eccentric position can be written as Fx , Fy. Making rotor in the seal a tiny displacement Δx in x-axis direction, force differences are ΔFxx ,ΔFyx. Similarly, giving rotor a displacement Δy in y-axis direction, force differences are ΔFxy ,ΔFyy. Thus, principal and cross stiffness can be written as:
Kxx=ΔFxx /Δx , Kyx=ΔFyx /Δx.
Kyy=ΔFyy /Δy , Kxy=ΔFxy /Δy.
For simulation, the seal clearance is 0.5mm, rotor diameter in the seal is 30mm with 50% eccentricity, seal length is 40mm. Simulation parameters are: pressure boundary in both inlet and outlet. The inlet pressure is 0.4Mpa, and outlet pressure is 0.101325 Mpa. The rotating speed is 150 rad/s. Liquid in the seal is water and temperature is 25 Celsius degree. The flow in seal clearance is considered to be laminar flow.
Simulation result shows that it is necessary to monitor the pressure exerted on rotor surface. Only when surface pressure doesn't change significantly in iterating steps, the calculated force is reasonable. Fig. 1~ Fig. 3 show the monitor of residual, force in x-axis and y-axis direction in each iterating step. Comparing with velocity and energy residual, it is difficult to converge continuity residual (seen in Fig. 1). Force monitors in x-axis and y-axis direction show that force doesn't change significantly when continuity residual reaches 3E-6. This residual value is far beyond the recommended one. Thus, monitering surface pressure is completely necessary.
Kxx=ΔFxx /Δx , Kyx=ΔFyx /Δx.
Kyy=ΔFyy /Δy , Kxy=ΔFxy /Δy.
For simulation, the seal clearance is 0.5mm, rotor diameter in the seal is 30mm with 50% eccentricity, seal length is 40mm. Simulation parameters are: pressure boundary in both inlet and outlet. The inlet pressure is 0.4Mpa, and outlet pressure is 0.101325 Mpa. The rotating speed is 150 rad/s. Liquid in the seal is water and temperature is 25 Celsius degree. The flow in seal clearance is considered to be laminar flow.
Simulation result shows that it is necessary to monitor the pressure exerted on rotor surface. Only when surface pressure doesn't change significantly in iterating steps, the calculated force is reasonable. Fig. 1~ Fig. 3 show the monitor of residual, force in x-axis and y-axis direction in each iterating step. Comparing with velocity and energy residual, it is difficult to converge continuity residual (seen in Fig. 1). Force monitors in x-axis and y-axis direction show that force doesn't change significantly when continuity residual reaches 3E-6. This residual value is far beyond the recommended one. Thus, monitering surface pressure is completely necessary.
In order to keep the prediction of stiffness coefficient accurate, the
tiny displacement of rotor should be small enough so that the force difference is propotional to displacement. Several displacements (8E-4mm, 1E-3mm, 0.002mm, 0.005mm, 0.006mm, 0.008mm) are selected. Forces in these dispacement are calculated. Linear fit of these forces, the stiffness coefficient can be figured out, seen in Fig. 4 ~ Fig. 7. The results are: Kxx=0.277 MN/m, Kxy=0.032 MN/m, Kyy=0.305 MN/m, Kyx= -0.045 MN/m. |
Predicting damping coefficient
Damping coefficient of annular seal is predicted by harmonic excitational method. In contrast to experiment which exerts harmonic force on rotor system and measures vibration response, the simulation firstly makes rotor whirl in a circle orbit. Then, the damping coefficient can be figured out through calculating force exerted on rotor surface. Whirling orbit can be written as:
x = Rsin(ωt +φ) , y = y0 + Rcos(ωt +φ) (1)
In this orbit, the rotor performs a circle whirl around position (0, y0). The whirling angular velocity is ω. At the begining (t = 0), the rotor locates at (Rsinφ, y0 + Rcosφ), where R is the whirl radius, φ can be setted according to practical beginning location. For example, if the beginning location is (-R, y0), then φ should be 3π/2. In this section, the parameter φ is kept to be 3π/2. The calculating force can be written as:
Fx = A + Bsin(ωt - θx) , Fy = C + Dcos(ωt - θy) (2)
where A, B, C, D, θx, θy can be getten from simulation result.
The motion equation of rotor can be written as:
[Fx, Fy] = [Fx0, Fy0] + [kxx, kxy; kyx, kyy][Δx, Δy] + [cxx, cxy; cyx, cyy][Δx', Δy'] (3)
A = Fx0, C = Fy0 (4)
where Fx0, Fy0 is the force exerted on rotor when it doesn't whirl. Equation (4) can be used to examine the accuracy of whirling simulation.
Combined equation (1), (2), (4) with (3), there is
[Bsin(ωt - θx), Dcos(ωt - θy)] = [kxx, kxy; kyx, kyy][Rsin(ωt +φ), Rcos(ωt +φ)] +[cxx, cxy; cyx, cyy][Rωcos(ωt +φ), -Rωsin(ωt +φ)] (5)
where kxx, kxy, kyx, kyy have already been caclulated above.
In equation (5), the sine and cosine coefficient in equation left should equal with that in right, then
cxx = (Bcosθx - kxyR)/Rω, cxy = (kxxR - Bsinθx)/Rω, cyx = (Dsinθy - kyyR)/Rω, cyy = (Dcosθy + kyxR)/Rω (6)
The following figures show a simulating case with 50% eccentricity.
x = Rsin(ωt +φ) , y = y0 + Rcos(ωt +φ) (1)
In this orbit, the rotor performs a circle whirl around position (0, y0). The whirling angular velocity is ω. At the begining (t = 0), the rotor locates at (Rsinφ, y0 + Rcosφ), where R is the whirl radius, φ can be setted according to practical beginning location. For example, if the beginning location is (-R, y0), then φ should be 3π/2. In this section, the parameter φ is kept to be 3π/2. The calculating force can be written as:
Fx = A + Bsin(ωt - θx) , Fy = C + Dcos(ωt - θy) (2)
where A, B, C, D, θx, θy can be getten from simulation result.
The motion equation of rotor can be written as:
[Fx, Fy] = [Fx0, Fy0] + [kxx, kxy; kyx, kyy][Δx, Δy] + [cxx, cxy; cyx, cyy][Δx', Δy'] (3)
A = Fx0, C = Fy0 (4)
where Fx0, Fy0 is the force exerted on rotor when it doesn't whirl. Equation (4) can be used to examine the accuracy of whirling simulation.
Combined equation (1), (2), (4) with (3), there is
[Bsin(ωt - θx), Dcos(ωt - θy)] = [kxx, kxy; kyx, kyy][Rsin(ωt +φ), Rcos(ωt +φ)] +[cxx, cxy; cyx, cyy][Rωcos(ωt +φ), -Rωsin(ωt +φ)] (5)
where kxx, kxy, kyx, kyy have already been caclulated above.
In equation (5), the sine and cosine coefficient in equation left should equal with that in right, then
cxx = (Bcosθx - kxyR)/Rω, cxy = (kxxR - Bsinθx)/Rω, cyx = (Dsinθy - kyyR)/Rω, cyy = (Dcosθy + kyxR)/Rω (6)
The following figures show a simulating case with 50% eccentricity.
From Fig. 9 and Fig. 10, force in x-axis and y-axis direction can be written as:
Fx = A + Bsin(ωt - θx) = -5 + 0.58sin(500t - 3.355), Fy = C + Dcos(ωt - θy) = 69.84 + 0.59cos(500t - 3.38) (7)
Considering equation (4), Fx0, Fy0 in this case is -7N and 69.6N. After whirling simulation, it is gotten that A is -5N and C is 69.84N. The equation (4) is well met. Difference between A and Fx0 may come from error of calculation since force in x-axis direction is far smaller than that in y-axis direction. According to result in equation (7) and equation (6), the damping coefficient can be figured out as:cxx = (Bcosθx - kxyR)/Rω = 1383.5 N*s/m, cxy = (kxxR - Bsinθx)/Rω = 469.14 N*s/m, cyx = (Dsinθy - kyyR)/Rω = -523 N*s/m, cyy = (Dcosθy + kyxR)/Rω = 1382.4 N*s/m.
where B=0.58,R=0.0000008, ω = 500, θx = 3.355, D = 0.59, θy = 3.38, Kxx=0.277 MN/m, Kxy=0.032 MN/m, Kyy=0.305 MN/m, Kyx= -0.045 MN/m.
Summary
A method about predicting rotor dynamic coefficient of annular seal is introduced. Harmonic excitation method is used to calculate damping coefficient. Under circle whirling, the calculating force change in a harmonic style, which fits theoretical prediction quite well. This method is conducted by CFD, and can be used to replace bulk flow model in predicting seal coefficient. Back to list >>