{'fraction_non_alphanumeric': 0.0233806055960138, 'fraction_numerical': 0.010732081257186662, 'mean_word_length': 4.4267886855241265, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A fast and accurate method to predict 2D and 3D aerodynamic boundary layer flows

A quasi-simultaneous interaction method is applied to predict 2D and 3D aerodynamic flows. This method is suitable for offshore wind turbine design software as it is a very accurate and computationally reasonably cheap method. This study shows the results for a NACA 0012 airfoil. The two applied solvers converge to the experimental values when the grid is refined. We also show that in separation the eigenvalues remain positive thus avoiding the Goldstein singularity at separation. In 3D we show a flow over a dent in which separation occurs. A rotating flat plat is used to show the applicability of the method for rotating flows. The shown capabilities of the method indicate that the quasi-simultaneous interaction method is suitable for design methods for offshore wind turbine blades.


Introduction
In this paper we will show the feasibility of a fast and accurate method for the prediction of aerodynamic forces on airfoils and wings. The application of this method will -after further development -be in an offshore wind turbine rotor aerodynamics simulation code. The use of this method requires little user expertise and computational effort, but will be able to compute in detail the unsteady aerodynamic characteristics of rotor blades. The simulation of separated flows will be possible.
The rotor aerodynamics simulation code under development is a combination of a panel method flow solver for the unsteady, incompressible, inviscid external flow (outer region) and an integral boundary layer solver for the unsteady viscous flow near the blade surface (inner region), see Figure 1. The interaction between these two flow regions in separated flows will be accounted for by a so-called viscous-inviscid interaction scheme.
The viscous-inviscid interaction method is a quasi-simultaneous method. The concept of it was developed by Veldman [1,2]. In this method, the boundary layer is solved together with an approximation of the external flow: the interaction law. This method is applicable to both 2D and 3D flow simulations.
In this paper we focus on the applicability of the interaction scheme.

Viscous-inviscid interaction
The viscous-inviscid interaction between the inviscid and viscous flow domains can be performed in several ways. Among the many methods that exist, two of them form the idea of the quasi- simultanous interaction method: the direct and simultaneous method. They are discussed briefly here.
In the direct method the velocity resulting from the external flow calculation (u e ) is prescribed to the boundary layer. The boundary layer calculation gives the displacement thickness (δ * ) to the external flow calculation. The main drawback of this method is that it fails to calculate separated flows due to the so-called Goldstein singularity [2].
Simultaneously solving the external and boundary layer flow solves the problem of the direct method at separation, but it is rather complex in terms of software and it lacks flexibility in flow modelling. The simultaneous method is applied in the XFOIL code of Drela [3].
Avoiding the complicated software of the simultaneous method can be achieved by solving the boundary layer together with an approximation of the external flow and subsequently solving the complete external flow (see Figure 2). This is the idea of the quasi-simultaneous method. In this method, an interaction law is needed as an approximation. See Veldman [2] for an overview of the method.
We will show that the addition of the interaction law to the set of 2D boundary layer equations prevents the singularity at separation to occur by keeping the eigenvalues of the system of equations positive.     The interaction law is defined in such a way that this has no influence on the converged solution 1 : where I are the interaction law equations and B the boundary layer equations. The superscript n indicates the number of the iteration. If convergence is achieved, δ * (n) = δ * (n−1) causing I to drop out of Equation 1.
The number of non-zero terms in I determines the approximation of the external flow [4]. If I = E, then the matrix is full and no approximation of the external flow is used. By setting more terms to zero, the approximation gets more crude and the convergence of (1) becomes slower. The interaction law is based on thin-airfoil theory, see also Coenen [4].

Two-dimensional interaction law
In 2D applications, the interaction law with only non-zero terms on the diagonal of I turns out to be the most robust [4]. The interaction law formulation becomes: where u e i is the boundary layer edge velocity at point i and u ext i the external velocity at point i. The constant c = 4 π∆x is the result of the thin-airfoil approximation of the external flow and forms the elements on the diagonal of I.

Three-dimensional interaction law
The formulation of the three-dimensional interaction law is based on thin-airfoil theory as well. The difference is that only the effect of the flow direction itself is taken into account in the approximation (we can do that because it is an approximation). The coefficient of the other direction is set to zero.
From 2D results, the interaction law, which only takes into account the effect of the external flow of the corresponding point, turns out to be a very good choice [4]. In the three-dimensional formulation this idea is used as well. The 3D interaction law formulation becomes: ve−v qe dz and q e the total value of the boundary layer velocity.

Effect of the interaction law on the set of equations
The coefficients of the interaction law depend only on the local grid size. The positive coefficient ensures that the line of I in Figure 3 has a positive slope. The slope of I is required to have an intersection with curve B. We can even pose that the slope of I needs to be larger than the  As a result of this, the system of equations composed of I and B in 2D has always positive eigenvalues. We will show this in the results. A detailed description of the requirements for the interaction law coefficient can be found in [5].
In 3D simulations, the sign of the eigenvalue is determined by the sign of the crossflow direction, independent of the interaction law equation. The limit of this coefficient to ∆y → ∞ reduces to the 2D coefficient. From this we can also expect that the 3D interaction law coefficients have to obey the same restrictions as the 2D interaction law coefficient.

Models
In the numerical method, we will test the quasi-simultaneous method by executing simulations of flows over a NACA0012 airfoil, 2D and 3D dented plates. For the airfoil simulations, the flow is assumed to be turbulent except for a small region around the stagnation point. For the simulations of a flow over a dented plate (2D and 3D), we assume the flow to be fully turbulent.

2D
The laminar model consists of the unsteady von Kármán integral equation and Thwaites' model is used for the closure of the shape factor H, where H is defined as H = δ * /θ. In the turbulent part we use the unsteady variants of the von Kármán equation and Head's Entrainment equation. For the determination of the turbulent shape factor H 1 we use the relation of Houwink and Veldman [6]. The skin-friction coefficient is given by Green [7]. The set of equations is given in more detail in [5].
The boundary layer is modelled in the wake and the computational domain is extended for 0.5 chord length.

3D
For the 3D simulations we use two steady integral boundary layer equations with the steady 3D version of Head's entrainment equation (see [8]). The same closure models for H 1 and c f as in 2D are used here.

Numerical Method
The quasi-simultaneous interaction scheme is very suitable for a modular composition of the numerical program. The external flow is calculated in a separate routine than the boundary layer flow calculation, see Figure 4 and 5. The resulting program is flexible when it comes to switching between different sets of equations or models. This is an advantage over XFOIL, because it is very likely that aerodynamic models for offshore wind turbine rotor simulations will change in the (near) future. The flows chart of the 2D and 3D numerical programs are shown in Figures 4 and 5.
The boundary layer and interaction law equations are solved with two solvers. Solver a solves the partial differential boundary layer equations and algebraic interaction law equation simultaneously. For solver b, the interaction law equation is substituted into the partial differential equations describing the boundary layer. This results in a set of equations of partial differential equations only. This is solved using a flux-splitting method which is also used by Van der Wees and Van Muijden [9].
Both the 2D and 3D program use a finite difference upwind scheme. In the 3D program, the direction of integration in crossflow direction is determined by the sign of the local value of the cross flow velocity (solver a) or the sign of the eigenvalues (solver b). In the 2D program it is always an upwind scheme which is possible because the eigenvalues are positive. The potential   flow is solved by applying a panel method. The method is equal to the method used by Coenen [4] in which also a quasi-simultaneous interaction method is applied.

Results
In this section we present results for the application of both 2D and 3D simulations.

2D
Flows over a NACA0012 airfoil and a dented plate have been simulated using solvers a and b.
In Table 1 the input parameters for simulations with the airfoil are given.    [4]. The mismatch with the experimental data for higher angles of attack is caused by the discretization error. From the figures it follows that for the finest grid the results are very close to the experimental data. Simulations with solver b with the grids of 241 and 481 points at angles of attack larger than 12.0 • , use solver a in a few points around the trailing edge. The simulations with the finest grids with solver a at high angle of attack (α > 12.0) take about 1000 steady iterations and one time step to converge. The computational time is about one minute on a modern pc. A flow over a dented plate is performed to show the behavior of the eigenvalues of the system of equations. Figure 7 shows the eigenvalues of this simulation. The figures show that both eigenvalues remain positive in a flow in which separation occurs. The effect of the dent is clearly visible. We can also conclude that both solvers yield the same result.

3D results
In this section we show the results for simulations past a 3D dented plate. The depth of the dent is chosen such that separation occurs in the dent. The number of points in streamwise direction is 80 and in cross flow direction the number of points is 60. The Reynolds number is 11.5 · 10 6 . Figure 8 shows the streamwise component of the displacement thickness δ * s at 25% and 50% width of the plate. The dent is in the middle and the 50%-line is the symmetry line. In the figure on the right differences between the two solvers are clearly visible. It is expected that this is due to a too coarse mesh in that area. In the 2D simulations it was also observed that fine meshes in separated regions are of great importance. Figure 9 shows the skin-friction coefficient at 25% and 50% width of the plate. Figure 9b shows that separation occurs inside the dent.

Rotation
In order to test the applicability for rotating flows a simulation is performed in which a flat plate rotates around (0, 0), see Figure 10. The model of the undisturbed external flow takes the rotation into account as well. All other input parameters are kept equal to those of the simulations without rotation. Results of simulations presented here are obtained with solver a. Solver b also gives converged results, but these are not presented here. Figure 11 shows the streamwise and crossflow boundary layer displacement thickness (δ * s and δ * n resp.) at y = 0.50 for a flow with Ω = 0.0, Ω = 0.1 and Ω = 0.5. From the figure it can be learned that rotation decreases the value of the streamwise boundary layer displacement thickness δ * s and increases the crossflow boundary layer displacement thickness δ * n . This is the result that we expected, because there is an increase in crossflow velocity compared to the non-rotating flat plate.
A smaller value of the displacement thickness will delay separation. For wind turbines this is a favorable property.

Conclusion
The quasi-simultanuous interaction scheme is suitable for the application in a rotor aerodynamics code. It is a fast scheme that is robust and can capture the aerodynamic features which are of interest for rotor applications very well.
In 2D the Goldstein singularity is not present at separation due to the well-chosen value of the interaction law equation. The 3D simulations do not show a singularity at separation either.
The method is suited to switch easily between different flow models due to its high modularity.