CFD Prediction of the Airflow at a Large-Scale Wind Farm above a Steep, Three-Dimensional Escarpment

Show more

1. Introduction

The author has developed the numerical wind diagnosis technique named RIAM-COMPACT [1] - [13] which is based on LES turbulence modeling. An exclusive license for RIAM-COMPACT has been granted by Kyushu TLO Company, Limited to RIAM-COMPACT Co., Ltd. (http://www.riam-compact.com/), a venture corporation which originated at Kyushu University. Efforts have been made to promote RIAM-COMPACT, mainly in the wind power industry (e.g., private wind power providers, local governments, and wind turbine manufacturers) in Japan.

On another front, open-source CFD software packages are more widely used than in the past. One of the most widely used software packages is OpenFOAM (Open Field Operation And Manipulation) [14] . OpenFOAM is an open-source CFD toolbox which has been released and distributed under the GNU GPL (General Public License) [15] by the OpenFOAM Foundation, a non-profit organization. OpenFOAM solves governing equations which are discretized with the finite volume method on a three-dimensional unstructured collocated grid that consists of arbitrary polygons. Applications of this software extend to a very wide variety of fields such as incompressible flows, compressible flows, multiphase flows, direct numerical simulations, chemical reactions, combustion, buoyant flows, conjugate heat transfer, Lagrangian particle simulations, molecular dynamics, direct simulation Monte Carlo studies, magneto-hydro-dynamics, and stress analyses. Furthermore, this open-source software is equipped with a wide variety of models that are comparable to those used in commercially-available software, including turbulence models ranging from RANS to LES to DES. Compared to commercially available software, the ability of this software to produce converging and numerically stable simulations is poor due to its meager grid generation capabilities. Another issue is that, compared to software developed by Japan’s national projects, the large-scale parallel computing capabilities of this open-source software is poor (i.e., Intra-node parallelization is not supported by the standard version of the software; the coordinate list format is used for the storage format of coefficient matrices; dynamic libraries are frequently used). Nonetheless, since this open-source software can readily be used to obtain CFD results with a moderate knowledge of C++ and the ability to run numerical computations, the business community, mainly small and medium-sized businesses, has started adopting this software in recent years.

In the present paper, simulations are performed for the airflow at a large-scale wind farm constructed above a steep escarpment using two software packages: RIAM-COMPACT (Turbulence model: Standard Smagorinsky LES) and OpenFOAM v.2.1.0 (Turbulence model selected for the present study: SST k-ω RANS). Subsequently, the results obtained from these simulations are compared.

2. Description of the Wind Farm

In the present study, airflow at Duogu Wind Farm, a large-scale wind farm in China is investigated (Figure 1). Located on this wind farm are 33 wind turbines (1.5 MW each; hub-height: 65 m; rotor diameter: 82.6 m) manufactured by Ming Yang Wind Power Group Limited, a Chinese manufacturer. China Ming Yang Wind Power Group Limited is the ninth largest wind turbine manufacturer in the world and held 3.7% of the global market share in fiscal year 2013.

3. Overview of the Software Packages and Numerical Simulation Set-Up

Numerical wind simulations are performed for airflow at Duogu Wind Farm, which is located above a steep escarpment (Figure 1), using RIAM-COMPACT (Turbulence model: Standard Smagorinsky LES) and OpenFOAM (Turbulence

(a) (b)

Figure 1. Large-scale wind farm investigated in the present study (Duogu Wind Farm). (a) Overall view; (b) Enlarged view. A group of 33 wind turbines are located above a steep escarpment. Arrows show the wind turbines under investigation in the present study.

Figure 2. Computational grid used for the present study.

model selected for the present study: SST k-ω RANS). Figure 2 illustrates the computational grid used for the study.

Next, the numerical simulation techniques and simulation set-ups that are used in the two software packages for the present study will be described. The present study uses the RIAM-COMPACT natural terrain version software package, which is based on a collocated grid in a general curvilinear coordinate system, to numerically predict local wind flow over complex terrain with high accuracy while avoiding numerical instability. In this collocated grid, the velocity components and pressure are defined at the grid cell centers, and variables that result from multiplying the contravariant velocity components by the Jacobian are defined at the cell faces. For the numerical simulation method, FDM (Finite-Difference Method) is adopted, and an LES (Large-Eddy Simulation) model is used for the turbulence model. In the LES model, a spatial filter is applied to the flow field to separate eddies of various scales into GS (Grid-Scale) components, which are larger than the computational grid cells, and SGS (Sub-Grid Scale) components, which are smaller than the computational grid cells. Large-scale eddies, i.e., the GS components of turbulence eddies, are directly numerically simulated without the use of a physically simplified model. In contrast, dissipation of energy, which is the main effect of small-scale eddies, i.e., the SGS components, is modeled according to a physics-based analysis of the SGS stress. For the governing equations of the flow, a filtered continuity equation for incompressible fluid and a filtered Navier-Stokes equation are used as follows.

$\frac{\partial {\stackrel{\xaf}{u}}_{i}}{\partial {x}_{i}}=0$

$\frac{\partial {\stackrel{\xaf}{u}}_{i}}{\partial t}+{\stackrel{\xaf}{u}}_{j}\frac{\partial {\stackrel{\xaf}{u}}_{i}}{\partial {x}_{j}}=-\frac{\partial \stackrel{\xaf}{p}}{\partial {x}_{i}}+\frac{1}{Re}\frac{{\partial}^{2}{\stackrel{\xaf}{u}}_{i}}{\partial {x}_{j}\partial {x}_{j}}-\frac{\partial {\tau}_{ij}}{\partial {x}_{j}}$

${\tau}_{ij}\approx \stackrel{\xaf}{{{u}^{\prime}}_{i}{{u}^{\prime}}_{j}}\approx \frac{1}{3}\stackrel{\xaf}{{{u}^{\prime}}_{k}{{u}^{\prime}}_{k}}{\delta}_{ij}-2{\nu}_{SGS}{\stackrel{\xaf}{S}}_{ij}$

${\nu}_{SGS}={\left({C}_{s}{f}_{s}\Delta \right)}^{2}\left|\stackrel{\xaf}{S}\right|$

$\left|\stackrel{\xaf}{S}\right|={\left(2{\stackrel{\xaf}{S}}_{ij}{\stackrel{\xaf}{S}}_{ij}\right)}^{1/2}$

${\stackrel{\xaf}{S}}_{ij}=\frac{1}{2}\left(\frac{\partial {\stackrel{\xaf}{u}}_{i}}{\partial {x}_{j}}+\frac{\partial {\stackrel{\xaf}{u}}_{j}}{\partial {x}_{i}}\right)$

${f}_{s}=1-\mathrm{exp}\left(-{z}^{+}/25\right)$

$\Delta ={\left({h}_{x}{h}_{y}{h}_{z}\right)}^{1/3}$

Because high wind conditions with mean wind speeds of 6 m/s or higher are considered in the present study, the effect of vertical thermal stratification, which is generally present in the atmosphere, is neglected (neutral atmosphere). The effect of the surface roughness is not taken into consideration, and the terrain surface, including the ground surface, is treated as a smooth surface. For the computational algorithm, a method similar to a FS (Fractional Step) method [16] is used, and a time marching method based on the Euler explicit method is adopted. The Poisson’s equation for pressure is solved by the SOR (Successive Over-Relaxation) method. For discretization of all the spatial terms except for the convective term in the Navier-Stokes equation, a second-order central difference scheme is applied. For the convective term, a third-order upwind difference scheme is used. An interpolation technique based on four-point differencing and four-point interpolation by Kajishima [17] is used for the fourth-order central differencing that appears in the discretized form of the convective term. For the weighting of the numerical diffusion term in the convective term discretized by third-order upwind differencing, α = 3 is commonly applied in the Kawamura-Kuwahara scheme [18] . However, α = 0.5 is used in the present study to minimize the influence of numerical diffusion. For the LES subgrid-scale modeling, the standard Smagorinsky model is adopted with a model coefficient of 0.1 in conjunction with a wall-damping function [19] .

In contrast, the SST (Shear Stress Transport) k-ω model [20] [21] [22] , which is a form of RANS model, is used for OpenFOAM in the present study. This model is a two-equation model in which transport by turbulence shear stress is taken into account and the eddy viscosity is determined by solving the transport equations for turbulence kinetic energy, k, and the specific dissipation rate, ω. In order to achieve high prediction accuracy and numerical stability, a k-ε model is applied for airflow away from the boundary layer, and

$\{\begin{array}{l}\frac{\partial k}{\partial t}+{U}_{j}\frac{\partial k}{\partial {x}_{j}}={P}_{k}-{\beta}^{*}k\omega +\frac{\partial}{\partial {x}_{j}}\left[\left(\nu +{\sigma}_{k}{\nu}_{T}\right)\frac{\partial k}{\partial {x}_{j}}\right]\\ \frac{\partial \omega}{\partial t}+{U}_{j}\frac{\partial \omega}{\partial {x}_{j}}=\alpha {S}^{2}-\beta {\omega}^{2}+\frac{\partial}{\partial {x}_{j}}\left[\left(\nu +{\sigma}_{\omega}{\nu}_{T}\right)\frac{\partial k}{\partial {x}_{j}}\right]+2\left(1-{F}_{1}\right){\alpha}_{\omega 2}\frac{1}{\omega}\frac{\partial k}{\partial {x}_{i}}\frac{\partial \omega}{\partial {x}_{i}}\end{array}$

${P}_{k}=\mathrm{min}\left({\tau}_{ij}\frac{\partial {U}_{i}}{\partial {x}_{j}},10{\beta}^{*}k\omega \right)$

$C{D}_{k\omega}=\mathrm{max}\left[2\rho {\sigma}_{\omega 2}\frac{1}{\omega}\frac{\partial k}{\partial {x}_{i}}\frac{\partial \omega}{\partial {x}_{i}},{10}^{-10}\right]$

${F}_{1}=\mathrm{tanh}\left[{\left(\mathrm{min}\left(\mathrm{max}\left(\frac{\sqrt{k}}{{\beta}^{*}\omega y},\frac{500\nu}{{y}^{2}\omega}\right),\frac{4{\alpha}_{\omega 2}k}{C{D}_{k\omega}{y}^{2}}\right)\right)}^{4}\right]$

${F}_{2}=\mathrm{tanh}\left[{\left(\mathrm{max}\left(\frac{2\sqrt{k}}{{\beta}^{*}\omega y},\frac{500\nu}{{y}^{2}\omega}\right)\right)}^{2}\right]$

$S=\left[\frac{1}{2}\left({\partial}_{j}{u}_{i}+{\partial}_{i}{u}_{j}\right)\right],\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\nu}_{T}=\frac{{\alpha}_{1}k}{\mathrm{max}\left({\alpha}_{1}\omega ,{F}_{2}S\right)}$

${\alpha}_{1}=\frac{5}{9},\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\alpha}_{2}=0.44,\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\beta}_{1}=\frac{3}{40}$

${\beta}_{2}=0.0828,\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\beta}^{*}=\frac{9}{100},\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\sigma}_{k1}=0.85$

${\sigma}_{k2}=1,\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\sigma}_{k1}=0.5,\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\sigma}_{k2}=0.856$

a k-ω model is applied for airflow in and in the vicinity of the boundary layer. The specific equations that are used for the models are shown as follows. As is the case for many existing commercial software packages, the SIMPLE (Semi-Implicit Method for Pressure Linked Equations) method, which is a simultaneous relaxation method, is adopted for the pressure-velocity coupling algorithm in OpenFOAM. The present study uses pisoFoam, which is one of the incompressible flow solvers available for OpenFOAM. The pisoFoam solver is based on the PISO (Pressure Implicit with Splitting of Operators) algorithm [6] , which has better convergence properties than the SIMPLE method. For the linear solver, the present study uses the GAMG (Generalized Geometric-Algebraic Multi-Grid) method. The second-order backward difference scheme is used for temporal discretization. For spatial discretization, the first-order upwind difference scheme is blended with the second-order central difference scheme (using the OpenFOAM Gauss limited Linear scheme with a value of 0.05) so that the numerical solutions become TV (Total Variation) stable. The SST k-ω model implemented in OpenFOAM is a high-Reynolds number RANS model. High prediction accuracy can be ensured by setting the location of the first layer of grid points, ${z}_{\mathrm{min}}^{+}$ , at ${z}_{\mathrm{min}}^{+}$ > 30 and additionally using the wall-function. For the values of k and ω over walls, a standard wall function for each of the variables is adopted.

Figure 3 illustrates the characteristic scales used in the present study: h is the difference between the minimum and maximum terrain elevations within the computational domain, U_{in} is the wind velocity (m/s) at the inflow boundary at height z_{max}. In the present study, the Reynolds number, U_{in}h/ν, is set to 10^{4} for the simulations in both software packages, where ν is the kinematic viscosity of

Figure 3. Characteristic scales used in the present study.

air (m^{2}/s). In addition, the same boundary conditions are set for the simulations in both software packages. At the inflow boundary, a vertical wind velocity profile is given with a power law exponent of 5. At the side and upper boundaries, slip conditions are applied, and convective outflow conditions are applied at the outflow boundary. On the terrain surface, a viscous boundary condition (no-slip boundary condition) is imposed. The time step is set to Δt = 2 × 10^{−3} h/U_{in} for RIAM-COMPACT and to Δt = 10^{−4} h/U_{in} for OpenFOAM.

4. Comparison of the Simulation Results from the Two CFD Software Packages

Figures 4-6 show comparisons of the simulation results from RIAM-COMPACT (Turbulence model: Standard Smagorinsky LES) and those from OpenFOAM (Turbulence model: SST k-ω RANS). Figure 4 illustrates the distribution of streamwise (x) wind velocity (non-dimensional) at the inflow boundary and 65 m above the terrain surface (wind turbine hub-height). Figure 5 illustrates the distribution of streamwise (x) wind velocity (non-dimensional) on a vertical cross-section which includes the location of the wind turbine under investigation (No.12). Figure 6 illustrates the vertical profile of streamwise (x) wind velocity (non-dimensional) at the location of the wind turbine under investigation (No.12). Here, these figures show the time-averaged values. The comparisons reveal that the trends of the results from the two simulations are generally in good agreement. In addition, at the investigated wind turbine location, the vertical profiles of the streamwise wind velocity from both of the simulations showed local increases compared to the inflow streamwise wind velocity (terrain effects on wind).

Next, the influence of the spatial grid resolution on the vertical profile of streamwise (x) wind velocity (non-dimensional) at the location of the wind turbine under investigation (No.12) was considered. For this purpose, the RIAM-COMPACT natural terrain version software package based on a large-eddy simulation (LES) was done (see Table 1). Through the comparison between Figure 7 and Figure 8, in order to predict the vertical profile of wind velocity at the location of the wind turbine location precisely, the minimum horizontal grid spacing should be at least 12 m.

Figure 4. Distribution of streamwise (x) wind velocity (non-dimensional) at the inflow boundary and 65 m above the terrain surface (wind turbine hub-height). (a) Time-averaged values over t* = 100 - 200 in non- dimensional time for RIAM-COMPACT; (b) Time-averaged values over t* = 50 - 100 in non-dimensional time for Open FOAM.

Figure 5. Distribution of streamwise (x) wind velocity (non-dimensional) on a vertical cross-section which includes the location of the wind turbine under investigation (No.12), Time-averaged values from the same time periods as those of Figure 4.

Figure 6. Vertical profile of streamwise (x) wind velocity (non-dimensional) at the location of the wind turbine under investigation (Figure 4 and Figure 5), Time-averaged values from the same time period as those of Figure 4 and Figure 5.

Table 1. Parameters used in the RIAM-COMPACT (Turbulence model: Standard Smagorinsky LES).

Recently, it has been reported that the utilization rates of wind turbines on wind farms situated on complex terrain fall short of expectations; that is, reports of damage and breakage of the exteriors and interiors of wind turbines as well as wind turbines with notably low power output have surfaced. Terrain-induced turbulence is considered as the major cause of these issues from the author’s research results [1] . The source of terrain-induced turbulence is small variations in the topographical relief in the vicinity of wind turbines at which turbulence is mechanically generated. Finally, an example of wind risk (terrain-induced turbulence) diagnostics is presented in the case of the wind turbine No.7 as shown in Figure 1. The dimensions of the computational domain in this case are 10.0

Figure 7. Distribution of streamwise (x) wind velocity (non-dimensional) on a vertical cross-section which includes the location of the wind turbine under investigation (No.12), Instantaneous values t* = 100 in non-dimensional time.

Figure 8. Vertical profile of streamwise (x) wind velocity (non-dimensional) at the location of the wind turbine under investigation (No.12).

(x) × 2.0 (y) × 4.7 (z) km, and the number of grid points is 501 × 101 × 101 points (approximately five million points). The minimum horizontal and vertical grid widths are 7.5 m and 2.75 m, respectively. Other simulation settings are the same as those used for the above simulation.

As shown in Figure 9(a), the distribution of the streamwise (x) wind velocity (non-dimensional) on a vertical cross-section which includes the location of the wind turbine under investigation suggests that the No.7 wind turbine in the figure are subject to significant influence from separated flow (terrain-induced turbulence) which is generated locally and periodically upwind of the wind turbine. In addition, the vertical profiles of the streamwise (x) wind velocity (non-dimensional) in Figure 9(b) do not follow the so-called wind profile power law; a large velocity deficit can be seen between the hub center and the lower end of the swept area. The power curve of a wind turbine is calculated using the value of the inflow velocity at the hub-center and without taking into account the presence of the wind turbine. Furthermore, the calculation assumes that the values of the velocity shear present in the vertical wind profile are equal to those in a wind profile power law with exponents of approximately 5 to 7. Therefore, when the values of the vertical wind shear deviate significantly from those predicted by the wind profile power law, it can be anticipated that the actual output of the wind turbine falls significantly below the output predicted by the numerical model used to create the power curve. It is likely that very large velocity shear will be an important research topic in the future in connection with the issues of the tower vibration and the fatigue strength of wind turbines.

5. Conclusion

In the present study, numerical wind simulations were performed for a large-scale wind farm (33 wind turbines) located above a steep escarpment in China.

(a)(b)

Figure 9. Distribution of streamwise (x) wind velocity (non-dimensional) on a vertical cross-section which includes the location of the wind turbine under investigation (No.7), Instantaneous values t* = 100 in non-dimensional time. (a) Overall view; (b) Enlarged view

For the simulations, the RIAM-COMPACT software package (Turbulence model: Standard Smagorinsky LES) and the Open FOAM software package (Turbulence model: SST k-ω RANS, selected from various turbulence model options available in this software package) were used. Comparisons of the simulation results revealed that their tendencies were in good agreement. Of particular note, the vertical profiles of the streamwise wind velocity from both of the simulations showed local increases (terrain effects on wind) at the investigated wind turbine location. Next, the influence of the spatial grid resolution on the vertical profile of streamwise (x) wind velocity (non-dimensional) at the location of the wind turbine under investigation was considered. As a result, in order to predict the vertical profile of wind velocity at the location of the wind turbine location precisely, the minimum horizontal grid spacing should be at least 12 m. Finally, an example of wind risk (terrain-induced turbulence) diagnostics was presented. The vertical profiles of the streamwise (x) wind velocity (non-dimensional) do not follow the so-called wind profile power law; a large velocity deficit can be seen between the hub center and the lower end of the swept area.

Acknowledgements

This work was supported by JSPS KAKENHI Grant Number 17H02053. The author expresses appreciation to them.

References

[1] Uchida, T. (2017) High-Resolution LES of Terrain-Induced Turbulence around Wind Turbine Generators by Using Turbulent Inflow Boundary Conditions. Open Journal of Fluid Dynamics, 7, 511-524.

https://doi.org/10.4236/ojfd.2017.74035

[2] Uchida, T. (2017) Large-Eddy Simulation and Wind Tunnel Experiment of Airflow over Bolund Hill. Open Journal of Fluid Dynamics, in Press.

[3] Uchida, T. (2017) High-Resolution Micro-Siting Technique for Large Scale Wind Farm Outside of Japan Using LES Turbulence Model. Energy and Power Engineering, 9, 802-813.

https://doi.org/10.4236/epe.2017.912050

[4] Uchida, T. and Ohya, Y. (2008) Micro-Siting Technique for Wind Turbine Generators by Using Large-Eddy Simulation. Journal of Wind Engineering & Industrial Aerodynamics, 96, 2121-2138.

https://doi.org/10.1016/j.jweia.2008.02.047

[5] Watanabe, K., Ohya, Y., Uchida, T. and Nagai, T. (2017) Numerical Prediction and Field Verification Test of Wind-Power Generation Potential in Nearshore Area Using a Moored Floating Platform. Journal of Flow Control, Measurement & Visualization, 5, 21-35.

https://doi.org/10.4236/jfcmv.2017.52002

[6] Watanabe, F. and Uchida, T. (2015) Micro-Siting of Wind Turbine in Complex Terrain: Simplified Fatigue Life Prediction of Main Bearing in Direct Drive Wind Turbines. Wind Engineering, 39, 349-368.

[7] Uchida, T., Ohya, Y. and Sugitani, K. (2011) Comparisons between the Wake of a Wind Turbine Generator Operated at Optimal Tip Speed Ratio and the Wake of a Stationary Disk. Modelling and Simulation in Engineering, 2011, Article ID 749421.

http://doi.org/10.1155/2011/749421

[8] Uchida, T., Maruyama, T. and Ohya, Y. (2011) New Evaluation Technique for WTG Design Wind Speed using a CFD-Model-Based Unsteady Flow Simulation with Wind Direction Changes. Modelling and Simulation in Engineering, 2011, Article ID 941870.

http://doi.org/10.1155/2011/941870

[9] Uchida, T. and Ohya, Y. (2011) Latest Developments in Numerical Wind Synopsis Prediction Using the RIAM-COMPACT CFD Model-Design Wind Speed Evaluation and Wind Risk (Terrain-Induced Turbulence) Diagnostics in Japan. Energies, 4, 458-474.

https://doi.org/10.3390/en4030458

[10] Uchida, T. and Ohya, Y. (2008) Verification of the Prediction Accuracy of Annual Energy Output at Noma Wind Park by the Non-Stationary and Non-Linear Wind Synopsis Simulator, RIAM-COMPACT. Journal of Fluid Science and Technology, 3, 344-358.

https://doi.org/10.1299/jfst.3.344

[11] Uchida, T. and Ohya, Y. (2006) Application of LES Technique to Diagnosis of Wind Farm by Using High Resolution Elevation Data. JSME International Journal [Environmental Flows], Series B, 49, 567-575.

[12] Uchida, T. and Ohya, Y. (2003) Large-Eddy Simulation of Turbulent Airflow over Complex Terrain. Journal of Wind Engineering & Industrial Aerodynamics, 91, 219-229.

https://doi.org/10.1016/S0167-6105(02)00347-1

[13] Uchida, T. and Ohya, Y. (1999) Numerical Simulation of Atmospheric Flow over Complex Terrain. Journal of Wind Engi-neering & Industrial Aerodynamics, 81, 283-293.

https://doi.org/10.1016/S0167-6105(99)00024-0

[15] http://www.gnu.org/licenses/gpl-3.0.en.html

[16] Kim, J. and Moin, P. (1985) Application of a Fractional-Step Method to Incompressible Navier-Stokes Equations. Journal of Computational Physics, 59, 308-323.

https://doi.org/10.1016/0021-9991(85)90148-2

[17] Kajishima, T. (1994) Upstream-Shifted Interpolation Method for Numerical Simulation of Incompressible Flows. Bulletin of the JSME, 60-578, 3319-3326. (In Japanese)

[18] Kawamura, T., Takami, H. and Kuwahara, K. (1986) Computation of High Reynolds Number Flow around a Circular Cylinder with Surface Roughness. Fluid Dynamics Research, 1, 145-162.

https://doi.org/10.1016/0169-5983(86)90014-6

[19] Smagorinsky, J. (1963) General Circulation Experiments with the Primitive Equations, Part 1, Basic Experiments. Monthly Weather Review, 91, 99-164.

https://doi.org/10.1175/1520-0493(1963)091<0099:GCEWTP>2.3.CO;2

[20] Menter, F. and Esch, T. (2001) Elements of Industrial Heat Transfer Prediction. 16th Brazilian Congress of Mechanical Engineering, November 2001.

[21] Menter, F. (1994) Two-Equation Eddy-Viscosity Turbulence Models for Engineering Applications. AIAA Journal, 32, 1598-1605.

https://doi.org/10.2514/3.12149

[22] Issa, R.I. (1985) Solution of the Implicitly Discretized Fluid Flow Equations by Operator-Splitting. Journal of Computational Physics, 62, 40-65.

https://doi.org/10.1016/0021-9991(86)90099-9