Some effects of domain size and boundary conditions on the accuracy of airfoil simulations

This paper investigates a specific case of one of the most popular fluid dynamic simulations, the incompressible flow around an airfoil (NACA 0012 here) at a high Reynolds number (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$6 \times 10^6$$\end{document}6×106). OpenFOAM software was used to study the effect of domain size and four common choices of boundary conditions on airfoil lift, drag, surface friction, and pressure. We also examine the relation between boundary conditions and the velocity, pressure, and vorticity distributions throughout the domain. In addition to the common boundary conditions, we implement the “point vortex” boundary condition that was introduced many years ago but is now rarely used. We also applied the point vortex condition for the outlet pressure instead of using the traditional Neumann condition. With the airfoil generating significant lift at incidence angles of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$5^\circ , 10^\circ$$\end{document}5∘,10∘, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$14^\circ$$\end{document}14∘, we confirm a previous finding that the boundary conditions combine with domain size to produce an induced (pressure) drag. The change in the pressure drag with domain size is significant for the commonly-used boundary conditions but is much smaller for the point vortex alternative. The point vortex boundary condition increases the execution time, but this is more than offset by the reduction in domain size needed to achieve a specified accuracy in the lift and drag. This study also estimates the error in total drag and lift due to domain size and shows it can be almost eliminated using the point vortex boundary condition. We also used the impulse form of the momentum equations to study the relation between drag and lift and spurious vorticity, which is generated as a result of using non-exact boundary conditions. These equations reveal that the spurious vorticity throughout the domain is associated with cancelling circulation around the domain boundaries.


Introduction
Predicting the flow around an airfoil plays an important role in multiple engineering applications, ranging from the design of wind turbine blades, propellers, and wings to turbomachinery and hydrodynamic engineering.Despite being a two-dimensional (2D) problem, airfoil simulation has a number of important challenges.
First, accurate results require using a suitable turbulence model.There are different approaches for solving turbulent flows, e.g., Direct Numerical Simulation (DNS), structured mesh gave more accurate drag predictions, but an unstructured mesh gave the lift to a comparable accuracy.Also, they reported that a structured mesh provides better convergence and higher resolution compared to an unstructured mesh.Regarding the domain shape, Lu et al. [6] reported that since the O-mesh cannot accurately capture the wake region, the lift and drag errors are higher for this mesh type compared to the C-mesh and H-mesh.Using a C-mesh with boundaries less than 30c from the airfoil, where c is the airfoil chord, is a common practice [7][8][9][10][11], but doing a domain independence study is essential to make sure that the boundaries do not affect the results.In this study, we show that using these domain sizes can result in significant errors in the aerodynamic forces.
Another important consideration is to select appropriate BCs.In 2D flows, the velocity disturbances decay asymptotically as the inverse of the distance from the body while they decay as the square of the inverse distance in three-dimensional (3D) flows.Consequently, determining the domain size and the outer or far-field BCs plays a crucial role in achieving accuracy in a 2D simulation when the domain size is comparable to 30c.An added complication is that the correct outer BCs for pressure and velocity depend on the solution, and there are no simple physical principles, such as the no-slip condition at the inner or body boundary, that can be used to set them for airfoil flow.
An important velocity disturbance for a lifting airfoil is the circulation given by the Kutta-Joukowski (KJ) theorem.To approximate this disturbance, Thomas and Salas [12] applied the "point vortex" (PV) BC, which comprises simple equations for the velocity perturbations along the domain boundaries due to a vortex placed at the aerodynamic center -the quarter chord position -with a strength given by the KJ theorem.Despite its simplicity, the PVBC is rarely used in modern airfoil simulations.Destarac's [13] Euler solutions of the flow around a NACA 0012 airfoil showed that reducing the domain size while keeping the same BCs increased the pressure drag.Applying the PVBC reduced the magnitude of the increase, which is proportional to the lift squared (as in 3D) and inversely proportional to the domain size.
The present study extends Destarac's [13] important insight by showing that the viscous drag varies much less with the domain size.Thus, the error in the pressure drag can be used to estimate the domain size needed to achieve a desired accuracy for the total drag.We demonstrate that the PVBC is easily implemented in a modern CFD software, OpenFOAM, leading to a significantly smaller domain size for a given accuracy than the common BCs.We also consider the associated issue of vorticity generation within the domain and along its boundaries by the choice of BCs.Moreover, an impulse formulation for the lift and drag is used to determine the significance of the spurious vorticity.
In the first part of Section 2, we survey the common BCs used for airfoil simulations and select four for further study.Our implementation of the PVBC is described and the choice of a baseline domain size explained.Section 2 is completed by explaining the choice of computational domain and domain size.Section 3 describes the numerical solutions and a demonstration of grid convergence.Section 4 contains the results and their analysis using the impulse form of the momentum equations.Section 5 lists the conclusions.

Boundary conditions
We employed four BCs that are commonly used for airfoil simulations: BC-1, BC-2, BC-3, and BC-4, as described in Table 1, applied to the velocity vector, U , pressure, p, turbulent viscosity or eddy viscosity, ν t , and viscosity-like variable, ν at the inlet (I), top (T), bottom (B), and outlet (O) boundaries of our square domains (Fig. 1).The reason for using rectangular domains rather than C-grids is given at the end of this subsection and Section 2.2.The values for the free-stream velocity, U ∞ , ν t , and ν are explained in Section 3.
BC-1 is based on the recommendation of Versteeg and Malalasekera [14].It comprises the Dirichlet boundary condition for U at I, T, and B and the Neumann bound- ary condition at O. For the pressure, BC-1 uses the Neumann boundary condition at I, T, and B and a fixed value at O.
One of the OpenFOAM tutorials [15] suggested BC-2 for an airfoil simulation.The "freestream" boundary condition in OpenFOAM is a mixture of the fixed value and the zero gradient boundary conditions.OpenFOAM chooses between these two conditions depending on the sign of the mass flux across the boundary.
BC-3 uses "slip" boundary condition [16] at T and B. The slip boundary condition for vector quantities enforces a zero-gradient condition for the tangential component and a zero value for the normal component.For scalars, it is the same as the "zero-Gradient" boundary condition in OpenFOAM.
BC-4, which uses the "symmetry" boundary condition [17] at T and B, can be regarded as analogous to BC-3.It imposes a parallel flow at boundaries as well.
Table 2 lists the three additional BCs we developed, partly through trial and error.The PVBC, which uses the transitional periodic boundary condition [18] at T and B, is easily derived from the Biot-Savart (BS) law.If we decompose the velocity to its unperturbed and induced components, � U = (U ∞ + u, v) , then u and v at a distance r from the point vortex are given by Table 1 Typical boundary conditions used in the airfoil simulations.The CV was the same for all cases BC-1   where θ is the angle between the induced velocity and the y-axis, and Ŵ is the circulation, which is defined as where is the surface whose boundary is C, d σ is the vector surface element, d l is the unit tangent vector to C, and is the vorticity vector, � = ∇ × � U.The drag and lift coefficients are defined as respectively, where D and L are the drag and lift per span (the z-direction), and ρ is the density of the fluid.Also, one of the most important equations in aerodynamics applications is the KJ theorem, which is Substituting the equation for C l into Eq.( 4) gives Using Eqs.(1) and ( 5), the PVBC for the velocity vector is Fig. 1 The computational domain, boundary labels, and coordinate system.The flow enters the inlet (I) on the left and exits the outlet (O).The top (T) and bottom (B) boundaries are also identified.The airfoil at the origin is shown for a positive α .Note that this domain has 1 unit length in the z-direction that is not shown where � x = (x, y) is the position vector of an arbitrary point on the boundary measured from the aerodynamic center.Neither Thomas and Salas [12] nor Destarac [13] clearly described their implementation of the PVBC for the pressure.Numerical experiments led us to the following for the outlet.Using the PVBC with Bernoulli's equation, the pressure is given by Ignoring the second order terms gives Equation ( 8) is applied at the outlet.As far as we know, this is the first use of the PVBC for the outlet pressure.
One reason for the unpopularity of the PVBC may be that it is iterative; the value of C l is required in Eqs. ( 6) and (8).We show below, however, that it is easy to implement the PVBC as defined in Table 2, in OpenFOAM.
Suppose we change perspective and consider the airfoil as one of an infinite cascade of lifting bodies separated vertically by distance 2A [5].In that case, the natural computational domain is rectangular, and the T and B BCs become periodic.This boundary condition is PeBC in Table 2.The lift and drag will be affected by this change [5], but periodic BCs are exact in the sense that the no-slip condition at the airfoil surface is exact, and so PeBC is a useful test case here.
To implement the PV equations for PeBC, we apply periodicity at the T and B boundaries and the following velocity components, which are modified by a straightforward extension of the BS law analysis to sum over all PVs in the cascade (6) • Table 2 The additional boundary conditions used in the present work.The CV was the same as in Table 1 Boundaries and the pressure at the outlet These equations are used in PVPeBC, the last case listed in Table 2.
The PVBC and PVPeBC were applied using the following procedure.First, we included the "forceCoeffs" function and set the appropriate values for "liftDir", "dragDir", "CofR", "pitchAxis", "magUInf ", "lRef " and "Aref " in the "controlDict" file.Adding this function generated a file called "coefficient" inside the "postprocessing" folder.The coefficient file stored C d and C l calculated in each time step.At every 100 iterations, C l was used to determine the boundary values of u, v, and p for the next 100 iterations.Changing the determination of C l from every 100 steps to every 1, 200, and 1000 steps was found not to alter the converged results or to have a significant effect on the execution time.

Computational domain and domain size
A rectangular domain is required for the periodic BCs, so we used this shape for all simulations.A square domain of sides 2A, where A is a multiple of c, has the advantage of specifying domain size by a single parameter, which is useful in estimating the error in C d .A rectangular domain also simplifies the application of the momentum equations in Section 4.3.
Many studies of airfoils use domains of dimensions less than 30 far from an airfoil [7][8][9][10][11].Therefore, A = 30 was chosen as the base case to determine the appropriate number of cells and begin assessing the various BCs.
For reasons that are explained in Section 4.3, most results were obtained with structured meshes.Results obtained with unstructured meshes will be highlighted as such.

Numerical method
We used OpenFOAM software [16] for all simulations.The simpleFoam algorithm for incompressible, steady, turbulent flows, was employed together with the Spalart-Allmaras one-equation RANS turbulence model.The aerodynamic center of the airfoil was at the origin of the coordinate system.The simulations were done for a range of angles of attack α = 5 • , 10 • , and 14 • to examine how variations in the lift influence the accuracy of the boundary conditions.These angles were specifically chosen to avoid flow separation.c = 1 m; the freestream velocity, U ∞ , was 51.48 ms −1 ; the kin- ematic viscosity, ν , was 8.58 × 10 −6 m 2 s −1 .The values for the turbulence viscosity, ν t , and ν are within the ranges recommended by Spalart and Rumsey [19], Spalart [20], and Menter [21].With these values, Re was 6 × 10 6 , which causes turbulent flow over most of the airfoil and so avoids having to simulate transition.( 9) Because it is a common value in modern practice, our initial simulations used A = 30 combined with BC-1.Without attempting to find the "best" A, Golmirzaee and Wood [5] found A = 500 with BC-1 and the Spalart-Allmaras turbulence model gave results very close to the experimental data from Ladson [22] for the NACA 0012.The computed C d and C l were 0.01220 and 1.07687, respectively, for α = 10 • and their finest grid.Table 3 also presents the experimental data [22] for C d and C l .The need to use a large A to achieve this agreement was part of the motivation for the present study.We aim to investigate the effect of the domain size and boundary conditions on the force coefficients.
A grid convergence study is necessary to quantify the numerical uncertainty arising from the choice of the number of cells.Roache [23] proposed a method for quantifying the uncertainty of the grid convergence using Richardson extrapolation.The resulting estimated values for C d and C l at zero grid spacing are 0.012937 and 1.07273 with the numerical uncertainty bounds of 0.03% and 6 × 10 −6 % , respectively.Using Roache's [23] grid convergence index, it was established that the data in Table 4 were in the asymptotic range.Figure 2 also shows the monotonic behavior of C d and C l with the number of cells.We used the fine grid and same turbulence model and only changed the BCs, as explained in Section 2.1.
Figure 3 shows that the grids were refined in regions with large spatial gradients, that is, around the airfoil and in the wake region.To ensure precise resolution of velocity gradients near the airfoil using the Spalart-Allmaras model, Eça et al. [24] suggested employing a value of y + = yu τ /ν < 1 , where y, ν , and u τ are the height of the first cell, kinematic viscosity, and friction velocity, respectively.Additionally, Eça et al. [24] recommended a value of y + ≃ 0.1 for the k − ω SST model.The values of min y + and max y + for different grid resolutions, namely coarsest, coarse, intermediate, and fine, are presented in Table 4.All further results in this study utilized the fine grid, which has a max y + of approximately 0.1, for all simulations.As presented in Table 4, the growth ratio represents the expansion of adjacent cells perpendicular to the airfoil for the initial  18, 36, 72, and 144 layers from the airfoil for the coarsest, coarse, medium, and fine grids, respectively.Furthermore, Table 4 shows the first grid height, which is twice the distance from the center of the cell next to the airfoil surface to the airfoil.This height was constant for all adjacent cells.

Results
This section discusses the impact of different boundary conditions and domain sizes on the execution time, velocity and pressure profiles, and lift and drag coefficients.Section 4.1 first compares the time, C d , and C l for four BCs that are usually used for an airfoil simulation.Subsequently, BC-3 is selected among these four conditions to investigate the effect of domain size on C d , C l , and the simulation time.Addition- ally, the results of PVBC were compared with those obtained from BC-3 at different domain sizes.We also derive an approximation for C d for A → ∞ , denoted as C d,∞ , along with an approximation for C l for three values of α.
In Section 4.2, we study the effect of different boundary conditions on the velocity, pressure, and spurious vorticities inside domains and along boundaries.Section 4.3 derives a relation between the generated vorticities, lift, and drag using the impulse equation.We indicate that PVBC allows for a reduction in domain size without losing accuracy, in contrast to BC-3.

Effect of domain size and boundary conditions on the lift, drag, and execution time
C d , C l , and execution time for BC-1, BC-2, BC-3, and BC-4 are listed in Table 5 for the baseline case of α = 10 • and A = 30 .All simulations were run on the University of Calgary's high-performance cluster, and no attempt was made to optimize the parallelization.Each simulation was run separately on one compute node with forty cores.Table 5 indicates that C l and C d for BC-2 were significantly different to the other BCs.This table also shows that BC-1 and BC-4 gave very similar results to BC-3.This was true for all A and so no further results for these two BCs are shown in the interests of brevity.Since the forces and execution time are very similar for BC-1, BC-3, and BC-4, whereas the longer running BC-2 is less accurate, we will take BC-3 as representative of current practice in applying BCs.Table 6 compares the execution time between BC-3 and PVBC for different values of A. Using PVBC increases the execution time compared to BC-3, but we conclude that the reduction in A allowed by the PVBC more than offsets the increase in execution time.
Table 7 compares C d , the pressure drag, C dp , the frictional drag, C df , and C l for BC-3 and PVBC for varying A and α .Over the large range of A in Table 7, C df for BC-3 varies much less than C dp ; a similar trend was found for BC-1, BC-2, and BC-4.BC-2 results are presented in Table 8 as an example.C df for the PVBC also has a weak dependence on A, but it takes a reduction to A = 3 for the PVBC pressure drag to rise significantly.
Table 9 shows C d,∞ and C l,∞ , the values of C d and C l as A → ∞ .They were calcu- lated using the least square extrapolation for BC-3 and PVBC.The values of C d,∞ and C l,∞ are in good agreement for the two BCs implying that either will give accurate results for sufficiently large A. C d,∞ and C l,∞ that are used in all subsequent equations refer to the values for BC-3, as listed in Table 9.
The lift and drag data for BC-3 and PVBC in Table 7 were processed in the following ways.Figure 4 plots the variations of C l with 1/A, and Fig. 5 depicts the relation- ship between C d = C d − C d,∞ and C 2 l /A .First, we note that the relative error in the lift changes less with A than does the drag error.Also, it is important that the trend of C l with A for BC-3 is opposite that for the PVBC.Since the trend of C d is opposite that for C l for BC-3, the lift-to-drag ratio, which often determines aerodynamic efficiency,  has a higher error than either the lift or drag.In the Introduction, it was mentioned that Destarac [13] found a similar behaviour in the pressure drag so the main new result is the near-constancy of C df .Destarac [13] also noted that C dp scales with 1/A  and his least squares fit to the data suggests C dp ∼ C 2 l /A .Using Desterac's scaling for BC-3 and the three values of α in Table 7, leads to with an error of less than 1% for A ≥ 5 .Thus, achieving better than a 2% error in C d at α = 10 • using BC-3, requires A = 91 , a value that is rarely used in modern practice.On the other hand, the PVBC requires only A = 5 for this level of accuracy.It is possible that Eq. ( 11) is independent of the airfoil shape, which would make it extremely useful.
Figure 6 shows the relationship between C l = C l − C l,∞ and C l /A .Using squares to derive the best fit for PVBC, we obtained with an error of less than 0.2% .Figure 6 shows that it is not possible to derive a similar data fit for BC-3 that uses only C l and A. The error in C l is much larger than for PVBC and is of opposite sign.

Velocity, pressure, and vorticity over the domain
The effects of the different BCs are apparent in the contour plots of the velocities, Fig. 7, and pressure, Fig. 8.We did not consider PeBC and PVPeBC in the previous section on airfoil results, but they are included here to clarify the role of the BCs.Note that all parts of the contour plots show the complete domain.The main differences in the velocities and pressure occur near the boundaries; applying the PVBC Eq. ( 8) at the outlet makes u(y) and v(y) equal to the corresponding inlet u(y) and −v(y) , respectively, even though the normal velocity gradient was set to zero at the outlet.This Neumann condition was used in all BCs, but was just shown to have little effect as part of the PVBC and the same is true for PVPeBC.It does have a big effect, however, for BC-3 where the velocities at the corners, u(A, A) and u(A, −A) have changed sign from those for PVBC.We return to the significance of the corner velocities in discussing the impulse analysis in Section 4.3.
The PVBC velocity at the inlet, Eq. ( 6), ensures that the inlet pressure matches the outlet pressure.For PeBC, the spurious U velocities near the outlet of Fig. 7 part (c) can be attributed to using the zero pressure condition instead of Eq. ( 10) at the outlet.Furthermore, the zero pressure condition leads to an abrupt pressure drop near the outlet in Fig. 8       as discussed in the next section.Interestingly, the BC that minimizes the spurious vorticity is PVPeBC, showing that the "exact" PeBC is not sufficient on its own to minimize the spurious vorticity but becomes effective when combined with Eq. ( 9) at the inlet and Eq. ( 10) at the outlet.Figure 10 plots the maximum spurious or "boundary-condition" vorticity, BC,max , determined at the inlet; (A 2 � BC,max c)/U ∞ is approximately constant for each α and independent of A. This constancy and the scaling of BC,max on C l means that the effect on the corner velocities mentioned above does not disappear as A increases, so Figs.7 to 9 show results that are qualitatively the same for all A.
It is emphasized that the levels of the vorticity contours in Fig. 9 were chosen to be very small to exclude the vorticity in the boundary layers and wake.To put the levels into context, assume that the maximum vorticity, max , in the airfoil boundary layer occurs at the surface.It is then easy to show that � max c/U ∞ = 2C f ,max Re , where C f ,max is the maximum friction coefficient.From the simulation results, C f ,max = 0.0315 for α = 10 • , from which it follows that � max c/U ∞ = 378, 000 .The magnitude of the boundary layer vorticity shed into the wake, , is U ∞ /δ te , where δ te is the upper sur- face boundary layer thickness at the trailing edge: �c/U ∞ ≈ 20 for α = 10 • .The small values are important, however, because they act over large areas -this is the reason for testing the A 2 scaling in Fig. 10 -and influence the momentum balances to be presented below.

Impulse equations for the lift and drag
Equations (3.55) and (3.56) of Noca [25] give the impulse equation for the force exerted by an incompressible flow on a body contained within an arbitrary CV.We need only the steady state version.He derived the equations by removing the pressure from the conventional Reynolds transport theorem for the momentum equation by using the ( In an ideal airfoil simulation, the vorticity is non-zero only in the boundary layers and wake, so we have and The first of the integral constraints of Eqs.(1.5) and (1.6) of Liu et al. [26] which are combined as Eq.(9.1.20) of Wu et al. [27] removes from Eq. ( 17), leaving To deal with the quadratic or nonlinear terms involving u 2 , v 2 , and uv, note that, the PV approximation gives u(−A, y) = u(A, y) , v(−A, y) = −v(A, y) , u(x, A) = −u(x, −A) and v(x, A) = v(x, A) , which reduces some of the quadratic terms.For example, the integral of u 2 along T is cancelled by that along B. uv is always antisymmetric about y = 0 for I and O, and about x = 0 for T and B, so its integral along each face is zero.Also, from conservation of mass, the sum of the U ∞ u and U ∞ v terms in Eq. ( 16) becomes zero.For a sufficiently large A, the PV approximation requires the quadratic terms to scale as C 2 l /A .This scaling suggests a simple origin for the error in drag given by Eq. ( 11), but we were unable to turn this observation into a prediction for the coefficient (0.0205) in Eq. (11).Neglecting the quadratic terms for sufficiently large A gives where the asymptotic form is Eq.(1.8) of Liu et al. [26].See also Eq. (9.1.15) of Wu et al. [27].Without the quadratic terms, Eq. ( 19) reduces to Eq. ( 4) for the lift, where We note, for future purposes, that the PV approximation requires all integrals on the right of Eq. ( 21) to be equal.( 16) O (U ∞ + u)�dy = 0 and The presence of spurious vorticity for BC-3 alters the lift and drag equations from Eqs. ( 14) and (15).At I, u = 0 and v = 0 .At T and B, v = 0 , so Eq. ( 14) becomes and, for lift, Again, the perturbations, u and v, must decay asymptotically at least as rapidly as 1/A to leave Now Green's theorem, which gives the relationship between the area integral of vorticity and circulation, cannot depend on the BCs, so Ŵ in Eq. ( 24) must equal the difference in circulation of the true bound vortical flow and the circulation due to the spurious vorticity outside the boundary layers and wake.To preserve the KJ equation, the latter circulation then must equal the sum of the two integrals which can only occur if Invariance also follows from the streamfunction ( ψ ) form of Eq. (25).Define (U ∞ + u) = ∂ψ/∂y and v = −∂ψ/∂x , and arbitrarily set ψ = 0 along B. It follows that because BC-3 also makes T a streamline where ψ(−A, A) = 2U ∞ A as determined at the inlet.Thus, I 1 is invariant.Part of the significance of invariance comes from Eq. (7.4.3) of Batchelor [28] which shows I 1 gives the gain in total head of the fluid along T due to the spurious vorticity.From Fig. 10, this gain scales as 1/A and so is CV-dependent and is clearly unphysical.Equation ( 23) can be rewritten as where Ŵ L is the lift-generating circulation due to the bound vorticity.Furthermore, the integral in the drag equation, Eq. ( 20), is altered by the boundary condition vorticity.Define I 2 as where the streamfunction form shows that I 2 is also invariant, but only for the flow ahead of the airfoil. ( Tables 10 and 11 show the values of the terms in Eqs.( 23) and ( 22), when α = 10 • and A = 30 .D/ρ and L/ρ from the impulse equation differ from the OpenFOAM direct determination of airfoil forces by less than 1% .Table 11 shows the terms in Eq. ( 22) eval- uated at x = 5 , 10, 20, and 30.As expected, the difference in the quadratic terms gets smaller with increasing x and is absorbed by the x-dependence of the integral in Column 5.That integral is, therefore, not invariant downstream of the airfoil and is added to I 2 at the inlet to determine D/ρ .From Eq. ( 20), the dominant term in D/ρ gives in combination with the Neumann outlet condition ∂v/∂x = 0 .The same condition applied to the vorticity equation gives Because the drag and the spurious vorticity are both positive, so are the integrals.Thus, u(A, −A) ≥ 0 , whereas the BS law requires u(A, −A) ≤ 0 , so u(A, −A) = 0 , and u(A, A) = −u(A, −A) , which requires = 0 .This argument explains the substantial dif- ference in U near the corner points in Fig. 7 between part (a) for BC-3 with non-zero and part (b) for PVBC with much closer to zero.The behaviour of the corner velocities is related to the cancellation of circulation along T and B for BC-3 which is implied by Fig. 7(a) and is quantified in the next paragraph.
Columns 2 and 5 in Table 10 suggest the invariance of I 1 and this was found to be the case at all x to within 1% .Again, the quadratic terms are very small, and L is dominated by I 1 and the Ŵ term, which is approximately one-quarter the magnitude of L/ρ .The rea- son for this is apparent from Table 12 which shows the contributions to Ŵ from the four CV faces for BC-3 and PVBC.As explained above, the contributions are nearly equal for the PVBC, but BC-3 has nearly cancelled the contributions from the inlet, top, and   bottom boundaries, leading to the approximately one-quarter contribution to L/ρ from the circulation in the last column of Table 12.The cancellation for T and B is associated with the interaction of BC-3 with the Neumann outflow condition on v.It is now clear that this cancellation is directly associated with the generation of BC , but we have not found any direct causal connection that would lead to a prediction of the magnitude of BC .
One further aspect of the generation of BC is that the inequality of the contributions to Ŵ shown in Table 12 implies that the irrotational flow has gained strain as well as vorticity.For a rectangular domain, the contributions to Ŵ from the top and bottom boundaries are from ∂u/∂y and those to the inlet and outlet from ∂v/∂x .Since must be constant along a streamline in the irrotational flow and ∂u/∂y = 0 at the inlet, and ∂v/∂x = 0 at the outlet, the magnitudes of rotation and strain must be comparable.The dynamics of such a flow may make an interesting study, but one that is obviated if the PVBC is substituted for BC-3.
So far we have assumed the spurious vorticity arises only from the BCs.It is possible, however, for it to be generated numerically for any BC, and we found this to be the case in our initial investigations using a structured mesh near the airfoil and an unstructured mesh for the outer flow.The differences in spurious vorticity are shown in Fig. 11 for typical simulations at α = 10 • plotted using contour levels much larger than in Fig. 9. Parts (a) and (b) of Fig. 11 show the structured and unstructured grids, respectively.As can be seen in part (c), the structured mesh generates little spurious vorticity, whereas part (d) reveals the presence of randomly generated vorticity in the unstructured mesh.Any vorticity that penetrates the structured grid closer to the airfoil is then convected downstream without change in level, presumably because the numerics of the structured  mesh prevent vorticity from being augmented or destroyed.The structured grid reduces the generation of spurious vorticity, enabling a more precise examination of different boundary conditions.For this reason, all other results in this study were obtained using structured meshes.

Airfoil pressure and skin friction
All the results presented so far relate to global quantities like lift and drag, whereas it is important to document the BC effects on the generation of forces around the airfoil.

Summary and conclusion
The imposition of boundary conditions (BCs) at the outer edge of the computational domain for a lifting airfoil simulation can never be exact as they depend on the solution.This is in stark contrast to the inner BCs at the airfoil surface where no-slip and impermeability provide exact conditions on the velocity and pressure.This paper considered high Reynolds number ( 6 × 10 6 ) flow around a NACA 0012 airfoil at three angles of The PVBC provides lift and drag values that are significantly more accurate at any domain size A than BC-3, which was found to be typical of the four common BCs that we tested.This improvement in accuracy extends to the local skin friction and pressure on the airfoil surface.As expected, the choice of domain size and BCs was found to be more critical in relative terms for the drag than the lift.The main outcomes of this study start with the confirmation of the finding by Destarac [13] that the BCs interact with the domain size to generate an (induced) pressure drag that scales with 1/A and the square of the lift coef-Our second outcome was that the skin friction drag is much less sensitive to the value of A. We quantified the domain size required for any desired level of accuracy in the simulation and suggest that with BC-3, A must be at least 90 at α = 10 • to achieve a 2% accuracy in the drag.This value of A is much larger than that used in many studies.
An important manifestation of inexactness was the generation of spurious vorticity within the computational domain.The magnitude of this vorticity was small, but it acts over a large area and influences the momentum balances for a control volume coinciding with the computational domain.All BCs generate spurious vorticity in the supposedly irrotational flow outside the airfoil boundary layers and wake.For the point votex condition, this vorticity was confined to the immediate vicinity of the top and bottom boundaries.The only relevant exact outer BC is the periodic one for a cascade of airfoils, rather than a single airfoil, and it was shown that this BC had to be supplemented with Using the impulse form of the momentum equations, the spurious vorticity throughout the domain was shown to be associated with the cancellation of circulation around the domain boundaries, but we could not formulate an expression for the error in lift that it caused.
The final main outcome of this study was the identification of the interaction of the BCs, the outlet Neumann condition on the velocity with the BCs on the top and bottom boundary.This interaction was associated with the cancellation of circulation along the top and bottom boundaries when BC-3 is used and the generation of spurious vorticity.It appears that one of the main benefits of the PVBC is that it mitigates the effect of the Neumann condition at the outlet velocity, which is almost universally used in airfoil simulations.
As future work, investigating the effect of different turbulence models on boundary conditions and domain sizes would be invaluable.Additionally, it would be interesting to understand whether different airfoil geometries have the same dependence of drag error on the lift and domain size.

Fig. 2
Fig. 2 Grid convergence study for a NACA 0012 airfoil at α = 10 • , A = 30 , and BC-1; a C d vs number of cells, and b C l vs number of cells

Fig. 3
Fig. 3 Representation of the intermediate grid at a the whole domain, b the leading edge, and c the trailing edge of a NACA 0012 airfoil.The domain size is A = 30 parts (a) and (c) in contrast to parts (b) and (d) of the same figure.

Fig. 7
Fig. 7 Comparison of the normalized velocity components for different boundary conditions for A = 30 .The parts show the entire domain with the airfoil at the center

Fig. 8
Fig. 8 Comparison of the normalized pressure for different boundary conditions for A = 30

Fig. 11 a
Fig. 11 a Structured and b unstructured grids for a NACA 0012 airfoil at α = 10 • with A = 30 and BC-3.Note that the grids at the boundary layer are structured for both parts.Normalized vorticity for the c structured and d unstructured grids

Figures 12 and 13
compare the pressure coefficient, C p , and the skin friction coefficient, C f , for three cases at α = 10 • : (i) A = 500 and BC-3, (ii) A = 5 and PVBC, and (iii) A = 5 and BC-3.These two figures indicate that using BC-3 and reducing the domain size from A = 500 to A = 5 results in a significant change for C p and C f .In contrast, PVBC with A = 5 achieves values of C p and C f that are similar to those for the larger domain size of A = 500 with BC-3.

Fig. 12
Fig. 12 Comparison of C p for airfoil simulations with A = 500 and BC-3, A = 5 and PVBC, and A = 5 and BC-3.a C p along the upper surface, b Zoom of C p along the upper surface, c C p along the lower surface, and d Zoom of C p along the lower surface.For parts b and d, every 5 points are shown

Fig. 13
Fig. 13 Comparison of C f for airfoil simulations with A = 500 and BC-3, A = 5 and PVBC, and A = 5 and BC-3.a C f along the upper surface, b Zoom of C f along the upper surface, c C f along the lower surface, and d Zoom of C f along the lower surface.For parts b and d, every 5 points are shown

Table 4 C
d and C l for a NACA 0012 airfoil at α = 10 • having BC-1, A = 30 and different numbers of cells, N. GR is the growth ratio

Table 5 C
d and C l for a NACA 0012 airfoil at α = 10 • and A = 30 for different boundary conditions

Table 6
Comparison of the execution time of a NACA 0012 airfoil with BC-3 and PVBC for different domain sizes

Table 8
Comparison of force components for a NACA 0012 airfoil at α = 10 • with BC-2 and different domain sizes

Table 11
Comparison of D/ρ between the impulse equation and OpenFOAM results for a NACA 0012 airfoil at α = 10 • having BC-3 and A = 30 .All quantities have units of m 3 s −1 .I 2 = 3.33163 m 3 s −1 at the inlet and the direct determination of drag gives D = 17.14205 m 3 s −1

Table 12
Comparison of the circulation between BC-3 and PVBC for a NACA 0012 airfoil at α = 10 • having A = 30