- About this Journal ·
- Abstracting and Indexing ·
- Aims and Scope ·
- Annual Issues ·
- Article Processing Charges ·
- Author Guidelines ·
- Bibliographic Information ·
- Citations to this Journal ·
- Contact Information ·
- Editorial Board ·
- Editorial Workflow ·
- Free eTOC Alerts ·
- Publication Ethics ·
- Recently Accepted Articles ·
- Reviewers Acknowledgment ·
- Submit a Manuscript ·
- Subscription Information ·
- Table of Contents
International Journal of Rotating Machinery
Volume 2013 (2013), Article ID 896148, 11 pages
Optimum Design of Oil Lubricated Thrust Bearing for Hard Disk Drive with High Speed Spindle Motor
1Department of Mechanical Engineering, Tokai University, 4-1-1 Kitakaname, Hiratsuka-shi, Kanagawa-ken 259-1292, Japan
2Department of Mechanical and Manufacturing Engineering, Faculty of Engineering, Universiti Malaysia Sarawak, 94300 Kota Samarahan, Sarawak, Malaysia
Received 29 July 2013; Accepted 15 October 2013
Academic Editor: Masaru Ishizuka
Copyright © 2013 Yuta Sunami et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This paper presents the application of optimization method developed by Hashimoto to design oil lubricated thrust bearings for 2.5 inch form factor hard disk drives (HDD). The designing involves optimization of groove geometry and dimensions. Calculations are carried out to maximize the dynamic stiffness of the thrust bearing spindle motor. Static and dynamic characteristics of the modeled thrust bearing are calculated using the divergence formulation method. Results show that, by using the proposed optimization method, dynamic stiffness values can be well improved with the bearing geometries not being fixed to conventional grooves.
HDD has been used as the main storage multimedia for electronics devices. Currently, HDD widely depends on oil lubricated hydrodynamic bearings. Hydrodynamic bearing is mainly supported by thrust and journal bearings. A schematic view of bearings in a 2.5 inch HDD is shown in Figures 1 and 2.
Hydrodynamic bearing gives better performance characteristics compared to conventional ball bearings as it has high dynamic stiffness with additional much higher damping effects. These damping effect characteristics provide smaller nonrepeatable run-out (NRRO). NRRO is the major contributor to the track misregistration in HDD read-write mechanism. Even though the repeatable run-out (RRO) of oil lubricated hydrodynamic bearing is higher than the ball bearing spindle motor, the RRO can be corrected by a read-write servo. Therefore, to reduce NRRO and increase spindle performance is to improve the bearing performance.
Currently, groove geometries that are being widely used in HDD thrust bearings are mainly a spiral or a herringbone grooved geometries. Several investigators have conducted numerical analysis predictions of these grooves for HDD bearing performance [1–5]. There were also some attempts at improving the dynamic stiffness and damping of HDD spindle by introducing permanent magnetic thrust plates into the bearing spindle structure  or introducing magnetic fluid as lubricants . However, there are very few attempts at finding an optimum design with a novel geometry to replace the conventional herringbone or spiral grooves for HDD. In the attempts at improving HDD performance, Arakawa et al.  proposed a nonuniform spiral groove for journal bearings to expand the critical bearing number for higher revolution speed of HDD spindle. The approach indicated that the novel nonuniform spiral grooves hydrodynamic bearings manage to increase the stability of rotation in high speeds.
This suggests that there is a probability of a further improvement if the groove geometry is not being fixed to any conventional grooves, either spiral or herringbone grooves. However, as far as authors know, there are no attempts at drastically improving the bearing characteristics changing the geometry and dimensions for oil lubricated 2.5 inch HDD spindle motor. Recently, HDD is being demanded to be thinner. If the characteristic of thrust bearing can be drastically improved, it is possible for HDD to be thinner. Therefore, in this paper, by adapting the optimization method initiated by Hashimoto, new optimum groove geometry and dimension to replace conventional grooves and increase the bearing performance of an oil lubricated 2.5 inch HDD had been calculated. The groove geometry and dimension of the thrust bearing for spindle motor are calculated using the hybrid method [9, 10] with the improvement of dynamic stiffness , being set as the objective function.
The optimization results showed that the dynamic stiffness of oil lubricated thrust bearing can be improved by introducing a new optimum geometry into the system. Vibration analyses were also presented to numerically verify the applicability of the improved bearings.
2. Geometrical Optimization and Analysis of Bearing Characteristics
2.1. Modification of Groove Geometry
The geometrical optimization process is carried out using the hybrid method [9, 10], a combination of direct search method and successive quadratic programming (SQP) method. SQP is a powerful method to obtain the optimum solutions for nonlinear optimization problems. However, it sometimes can not find the global optimum solution when the objective function has the multipeaks. Therefore, the direct search method is used first to roughly search for several local optimum solutions; then, SQP is used to find the local optimum solutions when we set the initial values obtained by the direct search method. The global optimum solution is finally obtained by choosing the maximum or minimum value. Furthermore, the confirmation of the optimum solution is needed to clarify the validity of the optimum solution. The validity of the optimum solution is confirmed by comparing the value of the objective function calculated from the optimum design variables with those from randomly selected design variables.
Figure 3 shows the design parameters used in this optimization. The initial prescribed design parameters for thrust bearing in this paper, based on a commercialized 2.5 inch form factor HDD design, are shown in Table 1. Taking the spiral groove geometry as a threshold groove geometry, the groove geometry parameters of angle variation through are initially set to zero.
From there, the groove geometry design parameters are changed arbitrarily using the cubic spline interpolation function . The cubic spline interpolation function is a cubic polynomial equation where each section connected to nodal points is a continuity of a second order derivative of the function.
Figure 4 shows the spiral groove geometry represented as initial values in the feasible region for optimum design of thrust bearing treated in this paper. Referring to the same figure, the initial line representing the spiral groove geometry, is shown with the label initial groove geometry while the new line of geometry after the spline interpolation function is calculated step by step as shown with the label th step. The groove geometry parameter of through is set starting from inner region towards the outer region. The feasible region for the groove geometry is divided equally into four regions in the radial direction as there were no significant effects on the optimization with regions over the divisions of four.
With the most inner point of the initial spiral groove curve fixed at point , the angle variation of groove geometry is then evolved gradually from initial stage to the final stage.
2.2. Analysis of Bearing Characteristics
In this paper, applying the boundary-fitted coordinate system to adjust the geometrical optimization method derives the calculation method of bearing characteristics. Moreover, in the process of analyzing the static and dynamic characteristics of the bearings, the perturbation method is applied to the Reynolds equivalent equation. The Reynolds equivalent equation can be solved using the Newton-Raphson iteration method, and the static component and dynamic component of pressure are obtained. This detailed calculation method is shown in the appendix.
The load-carrying capacity is obtained by the following integration: where indicates the atmospheric pressure.
The minimum oil lubricating film thickness is simultaneously determined from the equilibrium condition between the axial load acting on a thrust bearing and the bearing load-carrying capacity. The minimum oil lubricating film thickness is obtained by solving the following force balance equation:
The spindle motor’s friction torque on bearing surface is given by the following integration: where is the angular velocity of the rotating shaft.
The spring coefficient and damping coefficient can be obtained, respectively, by integrating the real and imaginary parts of the dynamic pressure components, , as follows: where is the squeeze frequency of the shaft, which is angular frequency of harmonic vibration when the bearing is small vibrated in the axial direction under equilibrium condition, in the HDD obtained through the perturbation method the and modified Reynolds equation .
Finally, the bearing stiffness, , is given by the following equation:
In this paper, the oil leak is not considered because the surround of thrust bearing of the spindle motor for 2.5 inch HDD, as shown in Figure 1, is sealed.
3. Formulation of Optimum Design and Vibration Analysis
3.1. Formulation of Optimum Design
In formulating the optimum design problems for groove geometry, the design variable of vector consisting of groove geometry design variables of angle variations through and groove dimension design variables of , , and as shown in Figures 3 and 4 is defined as follows: where , and are the groove depth, the seal ratio, and the ratio of groove width, respectively. The seal ratio means ratio of seal part, and the ratio of groove width means ratio of groove part of the total surface area.
Another design variable considered in this optimization is the number of grooves. However, since design variables of groove numbers are of discrete values, the optimization is preliminarily calculated step by step for the groove number from to . Then, the groove numbers, which gives the highest dynamic stiffness, is then selected as the optimum groove number. We had confirmed that the bearing characteristic, which is especially the dynamic stiffness, is decreased if the groove number is larger than and less than in both cases. In addition, it is difficult to manufacture the grooves on the bearing surface because the width of grooves becomes narrow due to increasing the groove numbers. Therefore, we chose the range of groove numbers from to . This is also to save the computational time for obtaining the optimum solution.
Furthermore, the constraint conditions based on a 2.5 inch HDD thrust bearing are expressed as where the details of the constraint functions are defined as in (8).
In (8), besides the constraint function of groove geometries and dimensions from g1 through , the minimum oil lubricating film thickness obtained by (2) must be larger than the allowable film thickness to avoid bearing contact and seizures shown as . Meanwhile, the final constraint function represents the constraint variable of the damping coefficient being kept positive at all times to avoid self-induced vibrations. These values of prescribed constraint values are shown in Table 2:
The objective function being set to improve the performance characteristics of the spindle in hydrodynamic bearing is the dynamic stiffness because thrust bearing plays an important role in supporting the total axial direction load of the HDD spindle. The dynamic stiffness, , of the thrust bearing can be obtained through the following equation:
The optimization problem for HDD thrust bearing is then formulated as follows:
3.2. Vibration Analysis
Once the optimization is being carried out, the vibration analysis is calculated to numerically verify the improvements. The vibration is conducted assuming that the oil film supporting the bearing is a viscously damped one degree free vibration model. The values of the spring and damping coefficient obtained by the optimum design are applied for the vibration response.
The equation of motion is given by where is the mass of the disk platters, and are the damping and spring coefficients obtained by the above mentioned optimum design, respectively, and is the impulsive external force being applied to the disk platters.
The external impulsive force being applied to the HDD disk platters is assumed as follows: where the natural angular frequency is expressed as ; the damped natural angular frequency is expressed as which is angular frequency of oscillatory wave form with viscous damping; the damping ratio and impulse are, respectively, expressed as and .
The nondimensional expressions for the gain and phase angle are, then, calculated using the following equations:
4. Characteristics of Optimized Grooves
4.1. Groove Geometry and Dimensions of Optimized Bearings
As can be seen from the figures, all geometries of the optimized bearings changed from initial spiral groove geometry to a new geometry, a spiral groove with bends in the outer periphery. In this paper, such geometry is called the modified spiral geometry.
Generally, it is widely known that spiral groove geometry has the ability to elevate the film thickness with its pump in effects as the revolution increases due to pressure generation. This explanation is applicable for the inner periphery of the optimized groove geometry. As the bearing rotates, the oil lubricants flows inwards, starting from the outer bends periphery, and flows along the grooves and collides with the seal, hence generating pressure. From all of the pressure distribution figures, the point where the groove and seal meet has the highest point of pressure generated.
However, the phenomenon is different from the outer periphery for opposite spiral groove geometries. Instead of generating the pressure and increasing the film thickness, the oil lubricant flows away from the outer bends towards the most outer radius of the bearing surface. Hence, as can be seen from the pressure distribution figures, there are no pressure peaks generated in the vicinity of the most outer regions of the bearing surface. This is because the pressure in the outer periphery of the bearing is neutralized to the atmospheric pressure. The pressure equivalent to the atmospheric pressure caused by the outer geometry bends will act as a force to pull down the bearing and lower the film thickness. With the combination of spiral geometry in the inner periphery and the opposite spiral bends in the outer periphery, the film thickness decreased, thus maximizing the dynamic stiffness of the spindle.
Even though all groove geometries basically possess the modified spiral geometry, Figure 5(a) shows that the region of spiral groove of optimized with μm is smaller than those of optimized with μm and μm. The geometry is concentrated in the outer regions only, since the seal radius ratio turned to the maximum constraint value being set in this paper, which is 0.8. The oil pressure distribution on the right hand side in the same figure shows pointed sharp pressure distribution. These peaks are the pressure where the groove part and seal part meet, a result of collisions between oil lubricant flow and seal.
For the groove geometries with allowable film thickness of μm and μm shown in Figures 5(b) and 5(c), the inner part of the bearing possesses a larger portion of spiral groove geometries compared to that optimized with allowable film thickness of μm. Even though the inner division creates comparatively smaller pressure peaks, this time the pressure is well distributed throughout the bearing surface.
If we focus on the inner periphery of the total geometry closely, that is, the end part of inner spiral just before the seal, we can see that another bend with the opposite direction of the outer bends started to appear. These inner bends however are smaller than the outer ones. This can be explained by the fact that in order to maintain the prescribed allowable film thickness, it is essential to increase the pump in effect of the inner spiral. It is insufficient with only the inner spiral geometry and outer bends. Therefore, another bend in the inner periphery, just beside the seal, appears to increase the film thickness of the bearing. The inner bends can be seen clearly for the case of allowable film thickness of μm. The combination of inner bends, outer bends, and inner spiral geometry helps to optimally balance the pressure and improve the dynamic stiffness.
The details of groove geometry and dimension of the optimized bearings from Figure 5 are shown in Table 3. As can be seen from the table, the values of geometry parameters changed from initial to a new groove geometry with slightly different angles of through depending on the allowable film thickness. However, in general all the groove geometry changes from spiral groove geometry to spiral with bends in the outer periphery.
The other design parameters of groove dimensions such as seal radius ratio, groove depth, and groove width ratio also change. The seal radius ratio shows that the lower the allowable film thickness is being set, the more the groove geometry surface will be occupied with the seal area. This eventually made the groove geometry to be concentrated at the outer region of the total bearing surface. On the other hand, for the groove depth and groove width ratio values, results show that the lower the allowable film thickness is being set, the smaller those values turned into.
4.2. Static and Dynamic Characteristics of Optimized Bearings
Figures 6(a), 6(b), and 6(c) show the relations between the rotational speed and the optimized values of dynamic stiffness, , minimum oil lubricating film thickness, , and friction torque, , respectively. The solid lines, alternate long and short dashed lines, and short dashed lines with dots represent the optimized bearings for μm, μm, and μm, respectively. The dots represent the characteristics values for each optimized bearing at that specific rotational speed. The short dashed lines without dots represent the initial spiral bearings.
The objective function for this study is the improvement of the HDD spindle dynamic stiffness. As can be seen from Figure 6(a), with the proposed optimization method, the dynamic stiffness values are improved compared to the values of conventional spiral bearings. The oil lubricated bearing for conventional spiral showed a trend of decreasing with the increase of rotational speed. However, optimized bearings still show an improvement where different rotational speed gives different stiffness characteristics but with all of them being improved and maintained compared to the conventional ones.
Figure 6(b) shows the relation of minimum oil lubricating film thickness for all optimized results with rotational speed. The minimum oil lubricating film thickness of conventional spiral shows an increase with rotational speed. However, compared to the initial spiral groove, the geometry of optimized bearings has a much lower film thickness. This means that the outer bends of optimized bearings have the ability to pull down the bearings towards bearing base. This is what contributes to the increase of the dynamic stiffnesses of the spindle. Moreover, the results obtained are attributable to the groove numbers as shown in Figure 5. Generally, the load carrying capacity is increased with an increase in the groove numbers. Therefore, the oil lubricating film thickness is increased with an increase in the groove numbers because we set constant load (N) in this optimization. On the other hand, the dynamic stiffness is increased with a decrease in the film thickness as shown in Figure 6(a). Therefore, in this optimization, the minimum oil lubricating film thickness comes close to the allowable film thickness to drastically improve the dynamic stiffness changing the groove numbers of , 11, 13 under the allowable film thickness of μm, 3.0 μm, 4.0 μm.
For the friction torque values shown in Figure 6(c), all optimized bearings show slightly larger friction torque values compared to initial spiral groove. However, the increase of friction torque is small if we were to compare the dynamic stiffness improvement obtained using the proposed hybrid method.
To confirm the effectiveness of the optimum design presented here, the values of design variables and objective function for the optimized bearing with μm are compared to those for randomly choosen design variables rotating at 10,000 rpm in Table 4. As shown in the table, the objective function of the optimized bearing shows maximum value compared with other values. Furthermore, we confirmed the values of the optimized bearing in the cases of changing allowable film thickness μm, 4.0 μm to be the largest using the same method.
4.3. Results of Vibration Analysis
Figure 7(a) shows the response of bearings towards impulsive force. The magnitude of the amplitude for impulsive force shows the highest value for conventional spiral groove bearings where a second order system wave response of an underdamped oscillatory response can be seen. On the other hand, the responses towards impulsive force for the rest of the bearings show smaller displacement values where a first order system of an overdamped wave form can be seen. This confirms that the optimized bearings have higher stability than the conventional spiral ones.
Figures 7(b) and 7(c) show the frequency response for the bearings. The optimized bearings have a smoother line of gain and phase angle. Gain shows that optimized 2.0 μm has the highest damping followed by 3.0 μm, 4.0 μm, and conventional spiral, respectively. This proves that the optimized bearings have better system responses.
In this paper, the optimum design method for 2.5 inch HDD was carried out to improve the dynamic stiffness of a thrust bearing spindle. The results can be concluded as follows.(1)Bearing characteristics using the proposed hybrid method showed that, compared with conventional spiral groove bearings, improvements of dynamic stiffnesses can be obtained when a new geometry of modified spiral groove geometry is introduced.(2)To maximize the dynamic stiffness is to set the maximum seal radius ratio, that is, by concentrating the geometry at the most outer periphery of the bearing.(3)Vibration analysis of optimized bearings showed better characteristics than the initial spiral bearings. Therefore, optimized bearings are expected to improve vibration characteristic of HDD.
6. Future Work
In this paper, we confirmed that dynamic stiffness of thrust bearing for 2.5 inch HDD was drastically improved numerically using the proposed hybrid method. However, it is necessary to verify the calculated data by experiments. We have prepared the experiments to verify the calculated data by using original experimental test rigs as shown in Figure 8. Figures 8(a) and 8(b) show the impulse and vibration experimental test rigs, respectively. The response of bearings towards impulse force, which can be applied using see-saw mechanism as shown in Figure 8(a), can be measured. On the other hand, it is very difficult to measure the spring and damping coefficients of spindle motor for 2.5 inch HDD because the displacement of the bearing is very small. We proposed the identification method  of oil film coefficients of bearings for HDD spindle using the response waveform obtained by vibration experiment using the test rig as shown in Figure 8(b). As mentioned above, therefore, we plan to submit the further study including these contents.
When optimizing the groove geometry, it is necessary to perform a sequential analysis of characteristics of a bearing with a groove geometry modified in succession by the finite difference method, but because of the polar coordinate system, direct processing is rather different to realize. Therefore, an analysis of the bearing stiffness is performed first by using a boundary-fitted coordinate system as shown in Figure 9 [12, 13] and transforming a complex groove geometry into a simple fanlike geometry. A boundary transformation function used for the transformation is given as
The following Reynolds equivalent equation can be obtained from the equilibrium between the mass flow rates of oil inflowing into and outflowing from the control volume due to the shaft rotation and the squeezing motion: where , indicate the mass flow rates, across the boundary of = const. and across the boundary of = const. as shown in Figure 10. On the other hand, indicates the mass flow rate due to squeezing motion inside the control volume.
, , and are expressed, respectively, as follows:where, in that case,
Then, in (A.2), subscripts l and 2 and I–IV indicate the domains in the control volume.
Assuming that variations of the bearing clearance are microscopic, the minimum oil lubricating film thickness and pressure can be expressed by the following equation:
In the equations stated above, indicates the amplitude of small variations of the oil lubricating film thickness and and express a static component and a dynamic component, respectively.
The substitution of (A.5) into (A.2) and the negligence of terms of small magnitude of of above second order allow for the introduction of two equations for terms of orders 0 and l as follows:where subscript 0 indicates a static component of the mass flow rate determined from (A.3a), (A.3b), and (A.3c), and subscript similarly indicates a dynamic component.
|:||Groove width (m)|
|:||Land width (m)|
|:||Damping coefficient of oil lubricated film (N·s/m)|
|:||Impulse external force (N)|
|:||Allowable film thickness (m)|
|:||Groove depth (m)|
|:||Minimum oil lubricating film thickness (gap between shaft and lower groove surface) (m)|
|:||Spring coefficients of oil lubricated film (N/m)|
|:||Dynamic stiffness of oil lubricated film (N/m)|
|:||Mass of the rotor (kg)|
|:||Rotational speed of spindle (rpm)|
|:||Oil lubricated film pressure (Pa)|
|:||Atmospheric pressure (Pa)|
|:||Static component of oil lubricated film pressure (Pa)|
|:||Dynamic component of oil lubricated film pressure (Pa/m)|
|:||Coordinate of radius direction (m)|
|:||Bearing outer radius (m)|
|:||Bearing inner radius (m)|
|:||Seal radius (m)|
|:||Friction torque on bearing surface (N·m)|
|:||Load-carrying capacity (N)|
|:||Vector of design variables|
|:||Equipartition space of (m).|
|:||Ratio of groove width;|
|:||Viscosity of lubricant (Pa·s)|
|:||Coordinate of circumferential direction (rad)|
|:||Coordinate of change based on boundary-fitted coordinate system (m)|
|:||Coordinate of change based on boundary-fitted coordinate system (rad)|
|:||Phase angle (deg)|
|:||Groove geometry design parameter; (displacement for th point of shift from the initial groove geometry to a new point) (deg)|
|Natural damped frequency (rad/s)|
|Squeeze frequency of the shaft revolution (rad/s)|
|:||Angular velocity of the shaft (rad/s).|
|max:||Maximum value of state variables|
|min:||Minimum value of state variables.|
- G. H. Jang, S. H. Lee, and H. W. Kim, “Finite element analysis of the coupled journal and thrust bearing in a computer hard disk drive,” Journal of Tribology, vol. 128, no. 2, pp. 335–340, 2006.
- K. Y. Park and G. H. Jang, “Dynamics of a hard disk drive spindle system due to its structural design variables and the design variables of fluid dynamic bearings,” IEEE Transactions on Magnetics, vol. 45, no. 11, pp. 5135–5140, 2009.
- T. Kobayashi and H. Yabe, “Numerical analysis of a coupled porous journal and thrust bearing system,” Journal of Tribology, vol. 127, no. 1, pp. 120–129, 2005.
- J. Zhu and K. Ono, “A comparison study on the performance of four types of oil lubricated hydrodynamic thrust bearings for hard disk spindles,” Journal of Tribology, vol. 121, no. 1, pp. 114–120, 1999.
- T. Asada, H. Saitou, and D. Itou, “Design of hydrodynamic bearing for miniature hard disk drives,” IEEE Transactions on Magnetics, vol. 43, no. 9, pp. 3721–3726, 2007.
- G. H. Jang and C. I. Lee, “Development of an HDD spindle motor with increased stiffness and damping coefficients by utilizing a stationary permanent magnet,” IEEE Transactions on Magnetics, vol. 43, no. 6, pp. 2570–2572, 2007.
- Q. Zhang, S. Chen, S. H. Winoto, and E.-H. Ong, “Design of high-speed magnetic fluid bearing spindle motor,” IEEE Transactions on Magnetics, vol. 37, no. 4, pp. 2647–2650, 2001.
- Y. Arakawa, S. Ikeda, T. Hirayama, T. Matsuoka, and N. Hishida, “Hydrodynamic bearing with non-uniform spiral grooves for high-speed HDD spindle,” in Proceedings of the IEEE 13th International Symposium on Consumer Electronics (ISCE '09), pp. 511–512, Kyoto, Japan, May 2009.
- H. Hashimoto and M. Ochiai, “Optimization of groove geometry for thrust air bearing to maximize bearing stiffness,” Journal of Tribology, vol. 130, no. 3, Article ID 031101, 2008.
- M. D. Ibrahim, T. Namba, M. Ochiai, and H. Hashimoto, “Optimum design of thrust air bearing for hard disk drive spindle motor,” Journal of Advanced Mechanical Design, Systems and Manufacturing, vol. 4, no. 1, pp. 70–81, 2010.
- M. Ochiai, Y. Sunami, and H. Hashimoto, “Experimental study on dynamic characteristics of fluid film bearing for HDD spindle motor,” in Proceedings of the 2nd International Conference on Design Engineering and Science, pp. 217–222, 2010.
- H. Hashimoto and M. Ochiai, “Theoretical analysis and optimum design of high speed gas film thrust bearings (static and dynamic characteristic analysis with experimental verifications),” Journal of Advanced Mechanical Design, Systems, and Manufacturing JSME, vol. 1, no. 1, pp. 102–112, 2007.
- N. Kawabata, “Study on generalization of calculations of lubricant flow using boundary fitted coordinate system (part 1, basic equations of DF method and case of incompressible fluids),” Transactions of the Japan Society of Mechanical Engineers C, vol. 53, no. 494, pp. 2155–2160, 1987.