International Journal of Rotating Machinery

Volume 2012, Article ID 372364, 11 pages

http://dx.doi.org/10.1155/2012/372364

## A Wake Model for the Prediction of Propeller Performance at Low Advance Ratios

Ocean Engineering Group, Department of Civil, Architectural and Environmental Engineering, University of Texas at Austin, Austin, TX 78712, USA

Received 14 February 2012; Revised 25 May 2012; Accepted 4 June 2012

Academic Editor: Francesco Salvatore

Copyright © 2012 Ye Tian and Spyros A. Kinnas. 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.

#### Abstract

A low order panel method is used to predict the performance of propellers. A wake alignment model based on a pseudounsteady scheme is proposed and implemented. The results from this full wake alignment (FWA) model are correlated with available experimental data, and results from RANS for some propellers at design and low advance ratios. Significant improvements have been found in the predicted integrated forces and pressure distributions.

#### 1. Introduction

For the prediction of propeller performance, the panel method shows great advantage over other numerical tools. On one hand, the panel method recovers the real blade geometry as much as possible, instead of simplifying the propeller blades into lifting surfaces as the Vortex lattice method (VLM) does. On the other hand, the panel method is far more efficient than Reynolds-averaged Navier-Stokes (RANS) simulation, not only in a computational sense, but also on grid generation. Generally, the panel method strikes a good balance on accuracy and computational efficiency. Close to the design condition, the results from panel method are in good agreement with either experimental measurement or RANS simulation.

As the advance ratio becomes smaller, the pitch of the wake sheets becomes smaller. Consequently, the wake shed from one blade could be close to the consequent blade and, thus, strongly affect the pressure distribution on the consequent blade. Therefore, in order to correctly predict the propeller performance at low advance ratios, an accurate wake model is necessary. Intensive investigation has been taken in the past two decades on modeling the wake of a propeller.

Greeley and Kerwin [1] invented a fast wake alignment scheme (PSF-2 type alignment). The wake is aligned on a coarse grid first, and then interpolated onto a fine grid. Pyo and Kinnas [2] presented a fully aligned wake model with biquadratic dipole panels. The high order wake model, although gained success on predicting the rollup of the wake and also the propeller performance near design condition, was more computational intensive than those with low-order panels. Later, Lee and Kinnas [3] developed a 3D wake model based on low-order panel method for marine propellers, which emphasizes more on the modeling of a developed tip vortex cavity. The tip vortex cavity is modeled with a helical vortex tube, of which the complete geometry is determined from requiring the pressure on its surface to be equal to vapor pressure. For the fully wetted cases, the tip vortex core is still necessary. Lee and Kinnas [4] also applied the same wake model to fully unsteady nonaxisymmetric flows. He [5] proposed a fully aligned wake model based on a VLM code, MPUF3A. In this wake model, the tip vortex core is no longer necessary.

Independently, Politis [6] formulated the unsteady motion of propellers with an unsteady wake model in an inertial frame of reference. Propellers under different, unsteady operating conditions (impulsive start, inclined inflow, and heave motion) were simulated. Greco et al. [7] proposed a sophisticated steady wake model, including the fully aligned near wake, helical intermediate wake and far wake disk, for analyzing the steady performance of propellers. This wake model was thoroughly validated with experimental data for propeller near design conditions, but the off-design performance was not addressed.

In this paper, the same idea as that in He’s [5] wake model is improved within the framework of a low-order panel method code, PROPCAV. The panel code without any aligned wake model is first applied to a 3D wing at low angle of attack (AOA) on the purpose of verification. With the new wake model and also the PSF-2 type wake model, the panel code is then applied to a 5-bladed propeller DTNRC 4381, especially at low advance ratios, and a 3-bladed propeller DTMB 4119. In the cases of a 3D wing at low AOA and propeller DTNRC 4381 at different advance ratios, the results from PROPCAV are correlated with RANS. In the cases of propeller, comparisons with available experimental measurements are also carried out.

#### 2. Formulation

The velocity flow field is decomposed into two components: inflow velocity and perturbation velocity, written as follows: where is the total velocity, is the inflow velocity, and is the perturbation velocity due to the presence of obstacles.

##### 2.1. Governing Equations

The methodology applied in this study is based on the inviscid potential flow theory: where is the perturbation potential. The perturbation velocity field, therefore, is governed by the Laplace equation: If the inflow is also irrotational, we can also define the total potential as follows: Appling Green’s second identity, the Laplace equation (3) can be written as the following boundary integrated form: where represents the surface of hydrofoil or propeller blades and represents the surface of trailing edge wake. is the Green’s function, which is defined as in 3D, and is the distance between the two points and . A low-order panel method using constant dipole panels is adopted to solve (5).

##### 2.2. Boundary Conditions

###### 2.2.1. Kinematic Boundary Condition

Applying the kinematic boundary condition to the solid surfaces, we obtain the following:

###### 2.2.2. Pressure Kutta Condition

At trailing edge, is unknown, which is subjected to the pressure Kutta condition: where the + sign denotes the suction side panel, and the – sign denotes the pressure side panel. The physical meaning of (7) is that free vortex sheets cannot bear pressure jump.

For steady state, (7) leads to the modified Morino-Kutta condition [10]: where is the vector from the control point of (−) panel to the (+) panel at the trailing edge.

For unsteady problem in 2D, or even steady problem in 3D, (7) has to be enforced through the unsteady Bernoulli equation: An iterative scheme has to be applied in order to enforce (9), proposed by Kinnas and Hsin [11], namely the Iterative pressure Kutta condition (IPK).

##### 2.3. Evaluation of Local Velocity

The local perturbation velocities in the computational domain are not required for evaluating the pressure on the blade surface, but evoked by the wake model. Theoretically, the local perturbation velocity at any point can be evaluated through the following singular integration: where is the same as Green’s function that is defined in (5). Practically, the numerical evaluation of (10) casts two questions: where and how to evaluate the velocity.

If the field point is located right on the vortex sheet, (10) contains singular integration. For low-order panel method, evaluating the velocity at the control point of a panel can avoid the singularity. Lee and Kinnas [3] implemented this idea in their unsteady wake model. In order to align the wake, the velocities at the nodal points of a panel are interpolated from the velocities at the control points. Repaneling was carried out at each time step.

However, the velocities are finite at the control points only if the panels are regularly and smoothly distributed in space. At the late stage of the unsteady free wake evolution, or in certain intermediate steps to determine the steady wake, the geometry of the vortex sheet could be entangled. The velocity evaluated at the control point of a panel could lead to a huge deviation from the exact value. Numerical tests show that evaluating the velocity at nodal point with certain desingularized vortex kernel can be more robust than the schemes mentioned above.

In 2D, the dipole-induced velocity in (10) is equivalent to the induced velocity due to a vorticity distribution with strength : where is the vector which is perpendicular to ,.

Krasny [12] adopted a desingularized kernel instead of in (11), which led to the *vortex blob method* in the study of the evolution of 2D free vortex sheet. This method, although simple, was successful in many 2D applications. The vortex-induced velocity is limited to a finite value, and the inherent Kelvin-Helmholtz instability is suppressed.

In 3D, the vortex sheet is discretized into constant dipole elements, which are equivalent to vortex loops with constant vorticity strength. Therefore, induced velocity due to the vortex loop can be evaluated using the Biot-Savart law:

Replacing the Biot-Savart kernel , with the Rosenhead-Moore kernel , leads to the 3D extension of the *vortex blob method. *Studies by Ramsey [13], He [5], and Lindsay and Krasny [14] show that this method can give out reasonable results for the 3D evolution of vortex sheets.

##### 2.4. Wake Alignment Scheme

The free vortex sheets shed from the trailing edge from a propeller blade are material surfaces, which have to be aligned with the local flow velocity. Alignment schemes play very important roles in determining the spatial geometry or distribution of the free vortex sheets. For the steady and unsteady cases, the alignment schemes are numerically different.

Engineering community is more interested in the steady propeller performances. As an analogy to the Euler-explicit scheme for the unsteady cases, many researchers align the vortex sheet using the velocity at the upstream nodal points:

Numerical tests show that this scheme is stable and converges fast. However, this scheme has the same deficiency as the Euler-explicit scheme for the unsteady alignment: it continuously enlarges the rotating radius around the center of rollingup, and thus cannot predict the correct location and size of the rolling-up region. Therefore, in order to get better prediction of the wake geometry, alignment scheme in a manner of trapezoidal rule is preferred. However, direct application of the trapezoidal rule turned out to be unstable.

Tian and Kinnas [15] proposed a pseudounsteady alignment approach when simulating the leading edge vortex (LEV) for delta wings with sharp leading edge. In the LEV model, the accuracy of the size of the rollup region is more important than that in the free wake model. In the study of this paper, although no LEV model is applied for predicting the propeller performance, the same idea is adopted for aligning the propeller wake.

As shown in Figure 1, consider a material line defined as follows:

We have which leads to

Discretizing (16) with finite difference scheme for the slopes and with the Euler-Explicit scheme for the unsteady terms, we have where , and . Denoting , and , we have

Notice that (18) is the alignment scheme for vortex sheet under horizontal uniform inflow. In other words, the inflow is parallel to the *x*-axis. The idea behind this approach is that the fixed coordinate has to be in the same direction of the dominant inflow, therefore, (18) behaves in an upwind manner numerically. For example, the slope of the curve at current point is evaluated using current point and the upstream point. However, in the case of the flow around a propeller, the dominant inflow would not be in the axial direction, but along helices with constant pitch. The rotational component of the inflow is important. The flow around the propeller is solved relative to a coordinate system which rotates with the propeller. In other words, the in (1) should be replaced by , where is the angular velocity of the propeller, and denotes the radial vector from the rotation axis to any points on the blades of the propeller. The essence of the alignment scheme described above is to decompose the total velocity vector into two components: component parallel to the inflow and component normal to the inflow; then keep the length of the vortex segment in the inflow direction unchanged. The advance ratio is defined as , where and are the rotational velocity in revolutions per second, and diameter of the propeller respectively.

Consider a point under cylindrical coordinates , with being the grid size parameter in the streamwise direction. If the wake is only aligned with the inflow, the consequent point will be at ; and . Therefore, the effective inflow direction in the Cartesian system on the segment connecting the ()th and the th points is where in Cartesian system is defined as

Now the averaged perturbation velocity can be decomposed into two components: where is the scalar projection of in Cartesian system onto , and is the vector component of normal to .

Defining and we obtaine the alignment scheme for the wake of propeller: where denotes the coordinates of the point at time step. It is easy to show that if the projection of the vector on is , the projection of the vector on is as well. Equation (22) is implemented in the full wake alignment (FWA) model. Although (22) can be used for unsteady problems, in the study of this paper, only steady alignment is addressed. Also, the symmetry or periodicity of the flow field is exploited.

When applying (22) for the steady alignment, the wake shape is recalculated through (22) repeatedly without resolving for the potentials until reaching the maximum number of iterations, or reaching a given convergence criterion. Then, the influence coefficients, due to the wake are reevaluated. With the updated influence coefficients due to the wake, the potentials on the blades and also on the wake are resolved for. Repeating this procedure will lead to converged wake shape and converged forces on the blades.

#### 3. Application to 3D Wing at Low Angle of Attack with Straight Wake

Before applying the FWA model, the basic framework of current panel code is first tested with a 3D wing at low angle of attack (AOA). A tip-swept wing, which is expected to resemble a propeller blade with a large skew angle, is investigated. As shown in Figure 2(a), sweep is superimposed onto an elliptic planform, which has an aspect ratio . All the sections of the wing have NACA 0010 thickness form without camber. A straight wake without alignment is attached at the trailing edge of the wing. In this case, no LEV model is applied.

The results at low AOA (AOA=2DEG) from current BEM model is validated with RANS, which is obtained through commercial software STAR-CCM+.

SST turbulence model is adopted for the RANS simulation, and the wall + is mainly controlled in the range between 30 and 80. For the convenience of grid generation, the trailing edge of the wing is made to be round, which has a radius of 0.5% of the chord length at each section. The RANS case took around 3 hours with 16 Intel Xeon 2.54 GHz CPUs to converge. The BEM code took 3 minutes on a single CPU of the same type to get the results.

Three sections are taken to compare the pressure distribution. As shown in Figures 2(b)–2(d), the pressure distributions on all the sections are in good agreement with both methods.

#### 4. Application to Propeller DTNSR 4381

The panel code is applied to a five-bladed propeller DTNSR 4381, which is investigated by Boswell [8] through experimental measurement. The design advance ratio of this propeller is 0.889.

For the simulations with the panel code, two wake models are adopted: PSF-2 wake model, which is a fast wake alignment scheme developed by Greeley and Kerwin, [1] and a full wake alignment (FWA) scheme, which is implemented in this study. For all the advance ratios, wakes of one revolution are modeled with the FWA scheme. For the PSF-2 model, wakes are trimmed at 2.5 R downstream with an ultimate wake disk. The hub of this propeller is also modeled through BEM panels, but the hub vortex core is omitted in this study.

RANS simulations of the same propeller are carried out with commercial code FLUENT, for and 0.4. K- SST turbulence model is adopted. QUICK scheme is used for spatial discretization, and SIMPLEC scheme is applied for pressure correction. 1.5 million hexahedral cells are used to simulate one sector of the domain with periodic boundary condition. It took 12 to 14 hours on 32 Intel Xeon 2.54 GHz CPUs for the residuals to converge to 1E-6.

In the experimental study by Boswell [8], only integrated forces were measured. Figure 3 compares the integrated forces predicted by PROPCAV and RANS, and the experimental measurements. We can see that that the predicted KT from the panel code with either wake model agree the well with both experimental data and RANS near design . However, as the advance ratio becomes smaller, the KT predicted by the PSF-2 alignment starts deviating from the experimental data and the RANS results. At the meantime, the KT from fully aligned wake model is still in good agreement with the experimental data. The panel code with either wake model tends to slightly overestimate the KQ near the design condition, but the results from PSF-2 alignment deteriorate at low advance ratios. Figure 4 shows the results of convergence analysis of the BEM/FWA model: different spatial discretization on the blade and on the wake leads to almost the same KT and KQ.

Figure 5 shows the aligned wake geometry at two different advance ratios ( and ). Clearly at low advance ratio, the wake sheets are strongly rolledup. Figure 6 plots the wake geometry predicted by PSF-2 alignment at the same advance ratios as in Figure 5. Obviously, PSF-2 alignment did not predict the contraction of the transition wake radius. The contraction rate of the transition wake is a user input parameter in PSF-2 alignment. We simply set it as zero because this parameter differs from case to case. On the other hand, Figure 3 indicates that near the design condition, whether contraction of the wake radius has been predicted (FWA) or not (PSF-2), the predicted total performance of the propeller will not be affected significantly. Figure 7 shows sections of the aligned wake cut at the plane, from different with different spatial discretization on the blade, where propeller frame of reference is defined as: *x* is along the axial direction, *y* is along the radial direction passing the middle of the root of the key blade, and *z* is normal to *x* and *y* directions. The rollup of the tip vortex can be observed. Also, the wake geometry with different spatial discretizations converged to similar results. Notice that in the current wake model, no tip vortex core model is adopted.

More detailed comparisons of the results between BEM and RANS are carried out in two different blade sections at two different advance ratios, as shown in Figure 8. The correlations among different numerical methods are remarkably good. Results from panel method with fully aligned wake model are closer to that of RANS, comparing with the PSF-2 wake model. This is consistent with the results for the total forces. Also for the FWA model, different spatial grids yield almost coincident results. The pressure distribution from and panels is not shown in Figure 8 so that the difference of the pressure from RANS and BEM/FWA can be clarified.

In the case of panels on the blade, the computational time for PROPCAV with PSF-2 wake model was only 1 minute on a single Intel Xeon 2.54 GHz CPU. The FWA model for the same discretization on the blade with 60 streamwise wake panels took 30 minutes to converge on the same type of CPU. Comparing with RANS, it is a big advantage to be able to achieve the results with the same order of accuracy in such a short time. It is also easy to parallelize the FWA model. For the wake model itself, linear parallel speedup could be achieved.

Figure 9 lists the convergence history of the FWA model for the Propeller 4381 at . The averaged moving distance of the points at the last streamwise station on the wake is monitored. Notice that once the wake geometry is updated, the boundary integral equation has to be resolved for the corresponding dipole strength on the wake. However, resolving for the unknowns is not carried out at each step, as mentioned before, inversion of the coefficient matrix is only evoked at selected steps.

It should be pointed out that when the advance ratio is low enough, the leading edge vortex (LEV) could be developed, and affect the pressure distribution and total integrated forces. Greeley [16] proposed a semiempirical LEV model for vortex lattice method. Tian and Kinnas [15] started to predict the LEV effects in a more rational way. At the same time, the boundary layer effects on the pressure distribution on the blade would be noticeable at high loadings. In order to obtain the correct KT and KQ, the boundary layer effects have to be taken into account correctly as well. Sun and Kinnas [17] developed a viscous/inviscid interaction approach by coupling the panel code with a boundary layer solver, XFOIL.

The panel code in this study neither included the LEV model nor coupled with the boundary layer solver. The viscous effects are considered by using an empirical viscous pitch correction [18] and a skin friction coefficient, which corresponds to that of a flat plate at the same Reynolds number based on the chord length. It is expected that incorporating the LEV model and the viscous/inviscid interactive scheme into current code will make current results even better.

#### 5. Application to Propeller DTMB 4119

Jessup [9] carried out experimental measurement of the wake flow of a three-bladed propeller, DTMB 4119. The BEM code with the FWA model is applied to this propeller, in order to see the effects of the desingularized parameter in the Rosenhead-Moore kernel when evaluating the dipole induced velocities. Three different values of are selected (0.03; 0.05; 0.075).

Figure 10 compares the predicted radial location of tip vortex for propeller DTMB 4119 at design advance ratio with the experimental measurement by Jessup [9]. All the results are close to the experimental data. It seems that smaller values of lead to better prediction of the location of tip vortex. Figure 11 shows the KT and 10KQ predicted with different values of . Clearly, the total performance of the propeller is insensitive to the value of , especially near the design condition. Under high-loading conditions, slight deviation among the results from different can be noticed. The reason could be that with small , the rollup of the tip vortex could be overpredicted. Figure 11 indicates that can be reasonable candidates for routine application of the FWA model.

More rational approachs to determine can also be proposed. It is possible to relate with the vorticity strength carried by each vortex segments on the wake. However, further investigation is required to implement this idea.

#### 6. Conclusions and Future Work

In this paper, a low-order panel method is used to predict the performance of propellers. A wake alignment model based on a pseudounsteady scheme is presented. In order to validate the basic frame work of the panel code without the fully aligned wake model, the method is applied first to a 3D wing at low angle of attack. The predicted pressure distribution agrees well with the RANS results. In case of the propeller working at design condition, both PSF-2 and FWA wake models predict reasonable KT and KQ. The results from FWA model are still in good agreement with the experimental data even as the advance ratio becomes smaller. Strong rollup of the wake at low is predicted by the FWA model. Detailed pressure comparisons with RANS show that FWA significantly improves the simulated pressure distribution.

The FWA wake model is also partially validated with the experimental measurement of the wake flow of propeller DTMB 4119 by Jessup [9]. Correct radial location of tip vortex is predicted. Sensitivity study indicates that different values of the desingularized parameter only have limited effects on the predicted total performance of propellers.

On the evaluation of KT and KQ, the viscous effects are simplified with an empirical viscous pitch correction [18] and a skin friction coefficient. In the near future, a viscous/inviscid interaction model will be included into the current model in order to treat the viscous effects in a more rational way. A LEV model will also be added into the current wake model, so that more accurate pressure distribution towards the tip of the propeller blade can be obtained.

Correlation of the FWA model with the latest experimental data can also be very helpful to improve the methods described in this paper. Typically the PIV measurement of the wake flow of the INSEAN E779A propeller by Di Felice et al. [19] reveals many detailed vortical structures. An investigation will be carried out if the geometry of this propeller is made available.

#### Acknowledgments

Support for this research was provided by the U.S. Office of Naval Research (Contract no. N00014-07-1-0616 and N0014-10-1-0931) and Phases V and VI of the “*Consortium on Cavitation Performance of High Speed Propulsors*” with the following current members: American Bureau of Shipping, Daewoo Shipbuilding and Marine Engineering Co. Ltd., Kawasaki Heavy Industry Ltd., Rolls-Royce Marine AB, Rolls-Royce Marine AS, Samsung Heavy Industries Co. Ltd., SSPA AB, Sweden, VA Tech Escher Wyss GmbH, Wärtsilä Propulsion Netherlands B.V., Wärtsilä Propulsion Norway AS, Wärtsilä Lips Defense S.A.S., and Wärtsilä CME Zhenjiang Propeller Co. Ltd.

#### References

- D. S. Greeley and J. E. Kerwin, “Numerical methods for propeller design and analysis in steady flow,”
*Transactions of Society of Naval Architects & Marine Engineers*, vol. 90, pp. 415–453, 1982. View at Google Scholar - S. Pyo and S. A. Kinnas, “Propeller wake sheet roll-up modeling in three dimensions,”
*Journal of Ship Research*, vol. 41, no. 2, pp. 81–92, 1997. View at Google Scholar · View at Scopus - H. Lee and S. A. Kinnas, “Application of a boundary element method in the prediction of unsteady blade sheet and developed tip vortex cavitation on marine propellers,”
*Journal of Ship Research*, vol. 48, no. 1, pp. 15–30, 2004. View at Google Scholar · View at Scopus - H. Lee and S. A. Kinnas, “Unsteady wake alignment for propellers in nonaxisymmetric flows,”
*Journal of Ship Research*, vol. 49, no. 3, pp. 176–190, 2005. View at Google Scholar · View at Scopus - L. He,
*Numerical simulation of unsteady rotor/stator interaction and application to propeller/rudder combination [Ph.D. thesis]*, UT Austin, 2010. - G. K. Politis, “Simulation of unsteady motion of a propeller in a fluid including free wake modeling,”
*Engineering Analysis with Boundary Elements*, vol. 28, no. 6, pp. 633–653, 2004. View at Publisher · View at Google Scholar · View at Scopus - L. Greco, F. Salvatore, and F. Di Felice, “Validation of a Quasi-potential flow model for the analysis of marine propellers wake,” in
*Proceedings of 25th Symposium on Naval Hydrodynamics*, St. John's, Canada, August 2004. - R. J. Boswell,
*Design, Cavitation Performance, and Open-Water Performance of A Series of Research Skewed Propellers*, Naval Ship Research and Development Center. Department of the Navy, 1971. - S. D. Jessup,
*An experimental investigation of viscous aspects of propeller blade flow [Ph.D. thesis]*, The Catholic University of America, Washington, DC, USA, 1989. - J. T. Lee,
*A potential based panel method for the analysis of marine propellers in steady flow [Ph.D. dissertation]*, MIT, 1987. - S. A. Kinnas and C.-Y. Hsin, “Boundary element method for the analysis of the unsteady flow around extreme propeller geometries,”
*AIAA Journal*, vol. 30, no. 3, pp. 688–696, 1992. View at Publisher · View at Google Scholar - R. Krasny, “Computation of vortex sheet roll-up in the Trefftz plane,”
*Journal of Fluid Mechanics*, vol. 184, pp. 123–155, 1987. View at Google Scholar · View at Scopus - W. D. Ramsey,
*Boundary Integral Methods For Lifting Bodies With Vortex Wakes [Ph.D. dissertation]*, MIT, 1996. - K. Lindsay and R. Krasny, “A particle method and adaptive treecode for vortex sheet motion in three-dimensional flow,”
*Journal of Computational Physics*, vol. 172, no. 2, pp. 879–907, 2001. View at Publisher · View at Google Scholar · View at Scopus - Y. Tian and S. A. Kinnas, “Modeling of leading edge vortex and its effects on propeller performance,” in
*Proceedings of the 2nd International Symposium on Marine Propulsors (SMP '11)*, Hamburg, Germany, 2011. - D. S. Greeley,
*Marine propeller blade tip flows [Ph.D. thesis]*, MIT, 1982. - H. Sun and S. A. Kinnas, “Performance prediction of cavitating water-jet propulsors using a viscous/inviscid interactive method,” in
*Proceedings of the SNAME Annual Meeting and Ship Production Symposium (SMTC '08)*, Houston, Texas, USA, 2008. - J. E. Kerwin and C. S. Lee, “Prediction of steady and unsteady marine propeller performance by numerical lifting surface theory,”
*Transactions of the Society of Naval Architects and Marine Engineers*, vol. 86, pp. 218–253, 1987. View at Google Scholar - F. Di Felice, D. Di Florio, M. Felli, and G. P. Romano, “Experimental investigation of the propeller wake at different loading conditions by particle image velocimetry,”
*Journal of Ship Research*, vol. 48, no. 2, pp. 168–190, 2004. View at Google Scholar · View at Scopus