Next Article in Journal
Gas Dynamics of Micro- and Nanofluidic Systems
Next Article in Special Issue
A Generalized Diffusion Equation: Solutions and Anomalous Diffusion
Previous Article in Journal
Numerical Simulation of Carbon Dioxide–Nitrogen Mixture Dissolution in Water-Saturated Porous Media: Considering Cross-Diffusion Effects
Previous Article in Special Issue
CFD Simulation of a Hybrid Solar/Electric Reactor for Hydrogen and Carbon Production from Methane Cracking
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Assessment of a RANS Transition Model with Flapping Foils at Moderate Reynolds Numbers

Department of Industrial Engineering and Mathematical Sciences, Polytechnic Marche University, 60131 Ancona, Italy
*
Author to whom correspondence should be addressed.
Fluids 2023, 8(1), 23; https://doi.org/10.3390/fluids8010023
Submission received: 22 November 2022 / Revised: 29 December 2022 / Accepted: 6 January 2023 / Published: 8 January 2023
(This article belongs to the Special Issue Recent Advances in Fluid Mechanics: Feature Papers, 2022)

Abstract

:
Numerical simulations based on a high-order discontinuous Galerkin solver were performed to investigate two-dimensional flapping foils at moderate Reynolds numbers, moving with different prescribed harmonic motion laws. A Spalart–Allmaras RANS model with and without an algebraic local transition modification was employed for the resolution of multiple kinematic configurations, considering both moderate-frequency large-amplitude flapping and high-frequency small-amplitude pure heaving. The propulsive performance of the airfoils with the two modelling approaches were tested by referring to experimental and (scale-resolving) numerical data available in the literature. The results show an increase in effectiveness in predicting loads when applying the transition model. This is particularly true at low Strouhal numbers when, after laminar separation at the leading edge, vorticity dynamics appears to have a strong effect on the forces exerted on the profile. Specifically, the transition model more accurately predicts the wake topology emerging in the flow field, which is the primary influence on thrust/drag generation.

1. Introduction

Starting from the pioneering work of Knoller [1] and Betz [2], the fluid dynamics aspects of oscillating foils have drawn the attention of researchers such as Burgers and von Kármán [3], Lighthill [4], and many others, whose endeavours have been extensively reported in a recent review by Wu et al. [5].
In particular, several experimental and numerical studies [6,7,8,9,10,11,12] have been carried out at moderate Reynolds numbers, for which transitional effects can be relevant. In this circumstance, the use of the standard, fully turbulent, Reynolds-averaged Navier–Stokes equations (RANS) approach can be erratic, while scale-resolving simulations, namely large eddy simulations (LES) and direct numerical simulations (DNS), are still prohibitively expensive, particularly when a large-scale investigation is required. Hybrid LES-RANS approaches [13,14,15,16] allowed the reduction in the computational cost while retaining large-scale solution accuracy through a mesh size relaxation at solid boundaries; despite this, the required computational resources remain closer to a LES rather than a RANS. A viable alternative consists of using RANS equipped with a transition model [17,18,19], which may require the solution of additional partial differential equations (PDE). Among all the possibilities, the local algebraic approach is well suited for unstructured solvers, because it does not increase the computational effort and, at the same time, significantly improves the prediction capabilities. For instance, when the Reynolds number is of the order of 10 4 or 10 5 , employing a transition model enables capture of phenomena such as formation of a laminar separation bubble over an airfoil and more accurate prediction of the onset of the stall. Despite the availability of results involving transition models applied to airfoils at moderate Reynolds number, even recently [20], the authors are not aware of studies that specifically consider the oscillating motion of airfoils using RANS and a local transition algebraic model. This is the main subject of this work. Note that in the literature it is possible to find flapping foils at moderate Reynolds number simulated solving both the two-dimensional Navier–Stokes (Ashraf et al. [21] at R e = 20,000 and Li et al. [22] at R e = 40,000) and the fully turbulent RANS governing equations (Young and Lay [23] at R e = 10,000–80,000 and Wang et al. [24] at R e = 13,800), so the purpose of this article is also to establish the limits of applicability of the latter approach.
In particular, the analysis proposed here is focused on performance evaluation concerning the prediction of loads, efficiency, and flow structure emerging from the dynamics of two-dimensional oscillating airfoils given by a Spalart–Allmaras (SA) RANS approach with and without the contribution of a local algebraic transition model referred to as SA-BCM. The choice of the SA model relies on its extensive use in aerodynamic applications and its suitability to simulate turbulent flows past an airfoil, all while having a low computational expense. Within the subject of oscillating foils, the same model has been adopted in [10] for comparing numerical results with experimental particle image velocimetry (PIV) data. Therein, the model performance are considered satisfactory; however, in Section 3 we will show how the approach presented here can further increase the agreement between simulations and experimental visualizations for moderate Reynolds numbers. As far as the transition model is concerned, the implementation considered in this work follows the mathematical description sketched in [25], which will be explored in more detail in the dedicated section (Section 2.2). The appeal for algebraic models stems from the absence of any additional transport equation for the intermittency field, which is instead modelled through a suitable empirical correlation with local flow quantities and embedded within the transport equation for an eddy viscosity-like variable. Its effect consists of determining the activation of the eddy viscosity-like production term corresponding to the regions exhibiting turbulent features. This algebraic transition approach proved to increase the accuracy of the solution considerably over a wide range of Reynolds numbers without significant variation in the computational cost [25,26,27].
Both models are implemented on top of a discontinuous Galerkin (DG) solver. This method preserves high-order accuracy independently from the grid characteristics, and it has been exploited here to ensure the resolution independence of the solutions through p-refinement. Moreover, as a side result, it is shown that the transition-modified SA model does not compromise the scheme robustness, despite the well-known stability issues due to the implementation of RANS equations within high-order frameworks.
As far as the airfoil motion is concerned, two oscillatory configurations have been considered: (1) a two-dimensional flapping motion, composed of a vertical translation with a superimposed rotation around a pivotal point and (2) a two-dimensional, high-frequency small-amplitude vertical oscillation with a constant angle of attack. These two prescribed motions allowed testing of the models with kinematic configurations producing rather different effects. Specifically, the flapping motion may be assimilated to the tail movement of oscillatory swimmers, such as tuna and sharks, and it has been thoroughly investigated as a propulsive mechanism, whereas high-frequency heaving has been proposed as a motion strategy aimed at drag reduction and flight stability improvement of air vehicles.
To validate the derived results, enabling the capture of consistent conclusions over the performance of the two models, suitable reference works have been used as a foundation. Concerning the flapping motion, the work of Schouveiler, Hover and Triantafyllou [8] has been taken as a benchmark for the evaluation of hydrodynamic forces and propulsive efficiency. In order to be able to make comparisons with their experimental results, the simulations were performed by changing the motion parameters in such a way to simulate Strouhal numbers in the range 0.1 < S t < 0.45 for two different values of maximum angles of attack, i.e., α m a x = 15 , 30 . A more detailed description of the parameters will be outlined in the dedicated section (Section 3.1). For the second motion condition, namely high-frequency plunging, the work of Visbal [11] and Krais et al. [12] have instead been used as a reference. In [11], a three-dimensional LES analysis was provided for a plunging SD7003 airfoil, considering different values of Reynolds number and motion parameters, together with validation of the experimental results [10]. The comparison between SA and SA-BCM models was carried out using the previously reported configurations with a frequency k = 3.93 at R e = 40,000 and k = 10 at R e = 60,000. A more detailed description will be given in Section 3.2. The work of Krais et al. [12] provides additional data, once again derived from scale-resolving simulations (ILES) employing a nodal DG solver adopting the plunging configuration characterized by k = 3.93 and R e = 40,000 defined in [11] as a test case.
The overall organization of the manuscript is as follows: in Section 2.1 the parametrization of the profile motion is reported. This is followed in Section 2.2 by an outline of the high-order DG method and the implemented turbulence models. In Section 3, the results obtained from simulating both NACA0012 (Section 3.1) and SD7003 (Section 3.2) are provided. Finally, in Section 4, the main findings are summarized and discussed.

2. Simulation Setup

2.1. Kinematic Framework

The airfoil motion is fully defined by the following periodic kinematic functions:
θ ( t ) = α 0 + θ 0 sin ( 2 π f t ) h ( t ) = h 0 sin ( 2 π f t + ψ ) ,
where α 0 is the mean angle of attack; h 0 and θ 0 are the vertical and angular displacement, respectively; f is the frequency of oscillation; and ψ is the phase shift between the heaving and pitching sinusoidal waves. In terms of velocity, the vertical oscillation reads as
u 2 t ( t ) = 2 π f h 0 cos ( 2 π f t + ψ ) = U 0 cos ( 2 π f t + ψ )
where U 0 = 2 π f h 0 is the maximum vertical velocity associated with the oscillating motion. Similarly, the angular pitching velocity is
ω 3 t ( t ) = 2 π f θ 0 cos ( 2 π f t ) .
The effective angle of attack α , according to our sign convention, is related to the other kinematic quantities through
α ( t ) = θ ( t ) + arctan u 2 t ( t ) U .
The maximum value assumed by α during the oscillating cycle is denoted with α m a x . An illustration of the various kinematic variables is given in Figure 1.
The amplitude-based Strouhal number of the oscillation ( S t ) and/or the reduced frequency (k),
S t = U 0 π U , k = π c S t 2 h 0 ,
are the non-dimensional numbers used to characterize the foil motion. In particular, S t represents the ratio between the maximum heave velocity U 0 and the free-stream speed U , except for an additional π in the denominator. Finally, the Reynolds number is, as usual, defined as R e = U c / ν , where c is the foil chord and ν is the kinematic viscosity. Note that c and U are hereafter used as reference scales of length and velocity.
In order to identify the position of the profile during its motion for the different cases explored, the notation adopted in [10,11] is followed. According to that, the instantaneous vertical position of the foil is expressed in terms of a parameter Φ , called the motion phase, which expresses the vertical coordinate as a percentage of the oscillating period. In this way, the vertical positions of interest will be denoted in the following way: Φ = 0 indicates the maximum upward displacement, Φ = 0.25 the neutral position during downstroke, Φ = 0.5 the maximum downward displacement, and Φ = 0.75 the neutral position during upstroke.

2.2. Numerical Framework

The numerical simulations were performed using a DG solver (see [28,29,30,31]). The advantage of using such a class of methods is related to the ability to provide accurate high-order solutions on stretched, curved, and unstructured computational grids. In this work, we consider a two-dimensional incompressible DG code, which is able to deal with a non-inertial reference frame to account for the foil oscillations. Note that the results here reported are not DG specific and the high-order accuracy, up to the seventh order, was only exploited in order to reach the numerical resolution independence. The governing equations are expressed in terms of the absolute velocity u i , a vector in the inertial frame observed by the non-inertial frame attached to the foil (see, for example, [32,33]). Consequently, the inter-frame velocity components are U i = u i t ( t ) + ε i j k ω j t ( t ) ( x k x 0 , k ) , where x 0 , k is the centre of rotation. Note that, for the sake of convenience, in the U i identity we adopted the same u i t ( t ) notation already used in Equation (2); however, here the heaving velocity should be interpreted as defined in the rotating frame and not in the inertial one. The resulting PDE system reads as
u i x i = 0 , u i t + u i u j x j = p x i + x j ν + ν t u i x j ε i j k ω j t u k + U j u i x j , ν ˜ t + ν ˜ u j x j = γ B C C b 1 S ˜ ν ˜ + C b 2 σ ν ˜ x j ν ˜ x j C w 1 f w ν ˜ d 2 + 1 σ x j ν + ν t ν ˜ x j + U j ν ˜ x j ,
where p is the pressure divided by the density, ν t is the kinematic eddy viscosity, d is the minimum distance from the wall, and ε i j k is the Levi-Civita tensor. Note that in Equation (6), taking advantage of the identity U i / x i = 0 and deviating from the process that is usually performed in the literature, the divergence of the convective fluxes does not contain the relative velocity, and all the non-inertial terms are moved to the right-hand side.
Following the standard RANS approach and Boussinesq’s hypothesis, the flow-governing equations are combined with a closure model for the turbulent Reynolds stresses. More precisely, the Spalart–Allmaras one-equation model has been used, with the eddy viscosity-like variable identified by ν ˜ . The eddy viscosity is computed according to the equation ν t = f v 1 ν ˜ , where f v 1 = χ 3 / ( χ 3 + C v 1 3 ) and χ = ν ˜ / ν .
In order to fully define the SA model, the following functions should be also specified:
S ˜ = S + ν ˜ k 2 d 2 , r = ν ˜ S ˜ k 2 d 2 , g = r + c w 2 r 6 r , f w = g 1 + c w 3 6 g 6 + c w 3 6 1 / 6 ,
where k = 0.41 is the Von Kármán constant and S = ( 2 Ω i j Ω i j ) 1 / 2 is the magnitude of the vorticity tensor, Ω i j = 1 / 2 ( u i / x j u j / x i ) . As performed in several SA implementations, the trip term is ignored in Equation (6). Thus, the base model is fully turbulent.
For the sake of compactness, the model constants, such as C b 1 , C b 2 , etc., will not be given in this paper. Interested readers can refer to the original work of Spalart and Allmaras [34]. In the following, the algebraic transition model proposed in [26] and more recently modified in [25,27] will be introduced. The latter version solves an issue regarding the Galilean invariance of the prior version, which is obviously crucial for the type of application handled here.
For the SA standard, fully turbulent model, the intermittency function γ B C is set to 1 in Equation (6), and for the transition model, it is evaluated as
γ B C = 1 e T 1 T 2 , T 1 = max R e θ R e θ c χ 1 R e θ c , 0 , T 2 = max ν t χ 2 ν , 0 ,
where χ 1 , 2 = 0.002 , 50 are calibration constants of the model. Thanks to the observation (see [35]) that the momentum thickness Reynolds number, R e θ , is proportional to the vorticity Reynolds number, R e v = d 2 / ν S , the former is evaluated as R e θ = R e v / 2.193 . Finally, the critical momentum thickness Reynolds number, R e θ c , is computed according to the empirical correlation R e θ c = 803.73 ( T u + 0.6067 ) 1.027 . Note that the value of T u is not assigned to the inflow/far-field boundary conditions, because the SA model does not consider the turbulent kinetic energy. As explained in [26], the term T 1 establishes the onset of transition, whereas the term T 2 propagates high levels of the intermittency function γ B C within the boundary layer.
Moreover, we decided to complement the model with the recent low Reynolds correction of Spalart and Garbaruk [34], which seems very well suited to deal with transitional flows. Following the latter approach, the model constant c w 2 = 0.3 in Equation (7) is changed to the function
c w 2 , L R e = c w 4 + c w 5 ( χ / 40 + 1 ) 2 ,
with c w 4 = 0.21 and c w 5 = 1.5 .
The DG discretization can be obtained by expressing the system of Equation (6) in variational form, decomposing the domain by elements, and defining an appropriate polynomial space of order q within each element. In this work, the convective numerical fluxes that are used to couple the solution at the interior faces of the mesh are computed through the artificial compressibility flux method of [36]. On the other hand, the second form of the Bassi and Rebay scheme (see [28]) has been used for the viscous terms. The DG discretized system can be written in the following compact form
M 0 d W d t = R ( W , t ) ,
where W is the global vector of the degrees of freedom (DOF), M 0 is a modified mass matrix with null entries corresponding to the pressure DOF, and R is the vector of residuals. This system of differential algebraic non-linear equations is integrated in time through a linearly implicit Rosenbrock-type Runge–Kutta scheme named ROS3P by [37]. The advantage of this class of schemes is that the non-linear solution process is completely avoided and the Jacobian matrix R / W has to be computed only once for each time step. The resulting time-discretized equations can be stated as
W n + 1 = W n + j = 1 s m j Y j , M 0 γ i i Δ t R W ( W n , t n ) Y i = R W n + j = 1 i 1 a i j Y j , t n + α i Δ t + + M 0 j = 1 i 1 c i j Δ t Y j + γ i Δ t R t ( W n , t n ) , i = 1 , , s
where the superscript n refers to the marching time step, s = 3 is the number of stages and m j , γ i , γ i i , a i j , α i are the coefficients of the scheme, which can be found in [37]. The last term of Equation (11) considers that the system is non-autonomous. In fact, due to the non-inertial frame, the absolute velocity vector u i changes at the non-slip, wall, and free-stream boundary conditions to follow the airfoil motion. Moreover, in Equation (6), U j and ω j t are not constant in time. At the wall, ν ˜ is set to zero, and the free-stream value depends on the choice of the model. Fully turbulent standard computations are performed with ν ˜ = 3 ν , and for the transitional simulations this value is set to 0.02 ν . These data are consistent with the recommendations found in the current literature (see, for example, [27,38]). The requirement to use smaller ν ˜ boundary data highlights the need of a strategy to prevent negative (spurious and non-physical) values of ν ˜ , which may result in the blow up of the simulations. This stability issue is exacerbated when a high-order scheme is adopted. Our approach is described in [30], where the details of our SA implementation are given.
The solution of the linear systems arising from Equation (11) were handled using a flexible generalized minimal residual (GMRES) solver preconditioned by a highly parallel efficient p-multigrid algorithm, described by [39].

3. Results

In the following section, we present the results obtained for two test cases that include a symmetric (NACA0012) and an asymmetric (SD7003) airfoil. The numerical grids, generated using the open-source software GMSH [40], are unstructured, as illustrated in Figure 2. In both cases, the outer boundary is a circle with a radius equal to 25 c . The boundary layer region is discretized using structured-like quadrangular elements suitably curved to accurately represent the foil’s curvature. To this end, a piece-wise third-order polynomial approximation of the faces of the elements is employed. The height of the first cell, expressed in wall friction units, is at most of the order of y + = 1 for all the computations here reported. The NACA0012 grid was significantly refined behind the airfoil because we paid particular attention (see the following section) to the wake topology. For this reason, the NACA0012 mesh consists of more elements: 7416 instead of 2875. We did not perform a classical grid refinement study, but rather we verified the numerical resolution by raising the polynomial order q of the solution. See Section 3.1 and Section 3.2.2 for examples of the solution convergence analysis. If not otherwise stated, the results here reported are obtained with q = 6 , ensuring a seventh-order space accuracy for the velocity components and the eddy viscosity and sixth-order accuracy for pressure. The time step size was also carefully verified, and in this work it ranges between 10 3 c / U and 10 2 c / U . In regard to the choice of the value of T u in the transition model, we did not find any clear statement in the references available for NACA0012 (e.g., [8]), thus we arbitrarily decided to use a small value such as 0.2 % . For the SD7003 computations, instead, we used T u = 0.1 % , in agreement with the experimental data reported in [10]. However, multiple tests showed that the results are quite insensitive to small changes of T u .
The forces F i acting on the airfoil expressed in the inertial frame and the momentum M 3 are non-dimensionalized using the chord c and the free-stream velocity U , obtaining the following drag, lift, and moment coefficients:
C D = F 1 1 / 2 ρ S U 2 , C L = F 2 1 / 2 ρ S U 2 , C M = M 3 1 / 2 ρ c S U 2 ,
where ρ is the fluid density and S is the reference area, which is the chord c times a unit’s span-wise length.

3.1. NACA0012

The first analysis focused on a sinusoidal pitching/plunging, two-dimensional NACA0012 airfoil at R e = 40,000. The motion laws are given by Equation (1), with α 0 = 0 . The heaving amplitude was fixed at h 0 = 0.75 c . Meanwhile, the axis of rotation was positioned at a distance of 1 / 3 chords from the leading edge. The heaving signal led the pitching harmonic by a quarter period, so ψ = 90 , which is known to be the optimal value in terms of efficiency ([7,8,9]).
Analysis was conducted through the comparison of propulsive performance and flow structure for several kinematic configurations. Specifically, we studied propulsive attributes at α m a x = 15 and α m a x = 30 while ranging S t from 0.1 to 0.45 with an interval of 0.05 . On the other hand, flow fields were inspected for cases at α m a x = 20 and α m a x = 30 , both with S t = 0.45 . The polynomial degree was adapted to the value of the frequency, with a minimum of q = 3 for low S t and a maximum of q = 6 for high S t , because as the Strouhal number increases, the vortical structures become smaller, requiring a greater spatial resolution. Results were then compared with the experiments of Schouveiler et al. [8] for validation, although the experimental data available in the literature appear to be scattered.
The series of surveys, aiming to perform a direct comparison to the models, are focused on the examination of propulsive performance sensitivity to S t , as heaving/pitching amplitude and frequency are the main factors influencing the wake topology [41,42,43,44]. The shape of the wake is in fact directly correlated with the thrust that the airfoil is able to create, as explained in [3]. This means that simulating a broad spectrum of combinations of θ 0 and f constitutes a solid benchmark to test the goodness of the numerical models.
Propulsive performance was assessed in terms of the thrust coefficient C T = C D and propulsive efficiency η p , defined as
η p = C ¯ T C ¯ P .
The terms C ¯ T and C ¯ P are, respectively, the mean thrust and mean power coefficient, computed as
C ¯ T = 1 n c T 0 n c T C T d t , C ¯ P = 1 n c T 0 n c T C L u 2 t U + C M c ω 3 t U d t ,
where u 2 t and ω 3 t are defined in Equations (2) and (3), respectively; n c indicates the number of cycles considered for the average; and T represents the motion period.
Note that the definition of Equation (13) admits negative values of η p when the flapping foil experiences drag, whereas it is positive when thrust is generated.
The results are collected in Figure 3.
At a high angle of attack, both models seem to perform in a quite similar fashion, retracing well the experimental data points at medium/high frequencies, especially if we consider η p . As f decreases, the gap between simulations and experiments starts to widen. Indeed, at S t = 0.1 we have that η p e x p 40.7 % . Instead, the transition and fully turbulent models predict efficiencies of 68.3 % and 63.0 % , respectively. Nonetheless, throughout the frequency spectrum analysed, SA always provides a moderately higher thrust than SA-BCM, which translates into a greater efficiency. This is probably attributed to the more diffusive nature of the fully turbulent model, which postpones separation and contributes positively to thrust generation. At a small angle of attack, the situation changes drastically, as the transition and fully turbulent models largely differ at low frequencies. For example, at S t = 0.1 , SA estimates a mean thrust coefficient of 0.100 , whereas SA-BCM provides a value of C ¯ T = 0.056 , which is very close to the experimental value of 0.061 . In addition to this, the efficiency of SA is characterized by η p = 56.2 % , whereas SA-BCM shows a value of η p = 26.9 % , much closer to the value of 37.1 obtained from the experiments. This seems to imply a better performance of the transition model, as it is also able to replicate the drop of the η p curve at small values of f, unlike the fully turbulent model.
In order to gain insight on this diverse behaviour at a small angle of attack, it is useful to exploit the phase space ( C L , C D ) representation composed of the instantaneous lift and drag coefficient curves, as proposed in [45], coupled with a visualization of the flow fields.
Figure 4 depicts ( C L , C D ) phase maps at the extremes of our frequency range for the case α m a x = 15 . As expected, when S t = 0.1 , the two models trace a completely different orbit. Instead, when S t = 0.45 , the curves are more or less overlapped. At high frequencies, the transition model moves between larger negative and positive values of the lift coefficient than the fully turbulent model and it is capable of revealing ‘richer’ dynamics as the airfoils generate resistance ( C D > 0 ), but the overall paths are coincident. Furthermore, SA-BCM does not hold a perfectly periodic trajectory due to the unsystematic activation/deactivation of γ B C within multiple periods. At small frequencies, the models follow their own trajectory, with SA-BCM featuring an imperfectly symmetric curve with respect to C L = 0 and a much longer portion in drag conditions ( C D > 0 ). The reason why this happens can be sensed by examination of the vorticity fields in Figure 5. It is evident how the transition model captures flow separation. Meanwhile, the classic Spalart–Allmaras model only causes a thickening of the boundary layer. At high frequencies, where just a trailing-edge vortex (TEV) detaches from the profile, this approximation is sufficient for replication of the wake pattern. At low frequencies, separation is anticipated, and there is the creation of a LEV that travels along the suction surface, eventually interacting with the TEV. The lack of this interaction when using SA causes little to no separation and leads to a streamlined wake. This supports the thesis of Guglielmini et al. [46], which states that in order to accurately evaluate performances of oscillating foils, flow separation must be modelled at the leading and trailing edges. Animations showing the physical mechanisms involved in the arise of flow instabilities are provided as Supplementary Material, although we will not focus on those aspects as it is beyond the purpose of this work.
The inability of the fully turbulent model to capture flow detachment is confirmed by Figure 6. At α m a x = 20 and S t = 0.45 , the transition model displays swirling structures forming at the leading edge, as portrayed in Figure 6f, instead in the fully turbulent model the flow remains attached (see Figure 6h). Despite this clear contrast, the fast dynamics appears to uniform the flow field at the trailing edge, yielding a good similarity in the wake topology, which also finds a good agreement with the visual recordings in [8]. This seems to imply that, for α m a x 20 , both models perform well if S t is high enough, but as soon as S t decreases, SA-BCM becomes much more accurate than SA. That said, Figure 6e,g still exhibit some minor differences. In particular, adopting the designation introduced by Williamson and Roshko [41], we have that SA-BCM draws a 2P + 2S wake, with a pair of corotating vortices (denoted by the letter ‘P’) plus a single vortex (denoted by the letter ‘S’) formed during each half cycle, whereas the wake resulting from SA falls into the 2P category, because the single vortex becomes embedded into the tail of the paired eddies. The only configuration where SA also reveals a separated flow upstream the trailing edge is at a high angle of attack and high frequency ( α m a x = 30 , S t = 0.45 ), due to the fact that the turbulence level is so high that even a standard RANS model is suitable. Indeed, SA-BCM and SA both induce a 2P + 2S wake, as shown in Figure 6a,c, which once again manifest a good resemblance to the visualizations obtained experimentally by Schouveiler et al. [8]. We point out that each wake in Figure 6 evolves in a 2S reverse von Kármán vortex street after 10 c for α m a x = 20 and 7 c for α m a x = 30 ; hence, they all generate thrust.
For completeness, it is worth finally providing an example of the convergence study of the solutions here reported. Space and time resolution-independence have been guaranteed through a polynomial refinement of the space discretization q and the choice of the time step size Δ t . It should be noted that by raising the degree of the polynomial approximation, the order of accuracy of the solution increases, although this also corresponds to a significant growth in the number of degrees of freedom. For instance, moving from q = 3 to q = 6 , the total number of variables goes from 296,640 to 830,592. An increase by a factor almost equal to 3 in the dimension of the numerical problem may be observed as a result of such a raise in the polynomial-order approximation. The solutions’ accuracy was verified by examining the behaviour of both the average performance coefficients and the ( C L , C D ) phase diagrams. For the sake of conciseness, in the following we only report the convergence results referring to the particular case with α m a x = 20 and S t = 0.45 . Its stringent requirement for high spatial and temporal resolutions due to the large S t makes it indeed a most challenging scenario. Considering the SA-BCM model, which presents relatively stricter accuracy constraints compared to the standard SA due to its algebraic modification, the average drag coefficient and efficiency have been observed to remain substantially unchanged throughout the performed space–time refinement: moving from q = 3 and Δ t = 10 2 c / U up to q = 6 and Δ t = 2.5 × 10 3 c / U , we have that C ¯ D ranges between 1.016 and 1.017 , whereas η p spans between = 0.548 and 0.551 . In Figure 7, the ( C L , C D ) phase diagrams for different polynomial approximations are reported over one cycle. It is easy to notice that for q 3 the results are almost indistinguishable. When only adopting SA-BCM, very small differences arise; however, as previously mentioned, with this model an imperfectly periodic behaviour of the solution is often observed, and small differences can be seen even between consecutive cycles performed at the same polynomial order, as will be further shown in Section 3.2.2. Note that another detailed space convergence study is also presented for the SD7003 profile in Section 3.2.2.

3.2. SD7003

For the SD7003 airfoil, the turbulence models were investigated by considering a high-frequency plunging at a constant angle of attack. Two different kinematic parametrizations have been explored; in both cases, the pitch angle of the profile was kept constant by switching off its time dependence in Equation (1), placing θ 0 = 0 .

3.2.1. Case 1

In this configuration, the airfoil was set to describe a plunging motion with reduced frequency k = 3.93 , amplitude h 0 = 0.05 c , and a constant angle of attack α 0 = 4 at R e = 40,000. The corresponding Strouhal number can be easily derived by inserting the values of k and h 0 in Equation (5), giving S t = 0.125 .
During the downward motion, the effective angle of attack increases as a consequence of the vertical velocity of the profile. In particular, substituting u 2 t into Equation (4) leads to α > 20 . This high value determines a boundary layer separation at the leading edge resulting in the formation of several LEVs trailing downstream.
Even though an accurate capture of the LEV dynamics, and thus an overall accurate flow description, is thought to be fundamental for precise computation of the forces exerted on the profile, the simulated configuration revealed to be solely dominated by the dynamics of the profile motion rather than the flow structure itself. Such an observation seems to be justified by the small difference in the forces predicted by the two models, even though SA-BCM allowed the derivation of a much more accurate representation of the relevant flow features with respect to the standard SA, especially when compared with the three-dimensional results in [11,12]. To demonstrate this, in Table 1 the thrust coefficients derived from the two models are reported as averages realized over five motion periods together with the aforementioned experimental and numerical reference values. It is possible to note that the models employed in the two-dimensional computations present no significant difference in terms of mean drag coefficient, placing both at a similar distance from the three-dimensional results.
This is further supported by examining the instantaneous behaviour of the hydrodynamic forces given in Figure 8, where the variations occurring in the phase plane ( C L , C D ) are jointly reported with the associated time dynamics over two periods. In these figures, it is possible to note that the curves related to the two models overlap almost entirely, aside from a minor offset in the trajectories drawn in the phase plane. Moreover, in Figure 8a, there is also a fairly good agreement between the C D signal of both models and the scale-resolving data from [11]. This suggests a substantial model independence in terms of instantaneous and mean force computations.
On the other hand, significant differences can be found by examining the details of the vorticity dynamics at the solid surface and the associated vortex development as a result of boundary-layer instabilities. From the instantaneous contours of the span-wise vorticity reported in Figure 9 and Figure 10, it can be noted that while the results of the two models in the aft portion of the profile are fairly congruent, a quite significant deviation is observed at the fore region close to the leading edge.
Considering the transition model computations, it is possible to see that when the profile reaches the top position Φ = 0 , (Figure 9a,b), two vortices develop over the lower surface of the profile due to boundary-layer separation, discernible by the regions of strongly negative vorticity generated by the wall-enforced no-slip condition. On the upper surface, the boundary layer appears to be fully attached on a relatively large zone past the leading edge, but a strongly diffused vorticity structure, a remnant of the previous cycle, travels downstream. Three-dimensional reference data in [11] clearly highlight the span-wise homogeneous character of the laminar boundary-layer instabilities at the onset corresponding to the front part of the airfoil. As a consequence, this region exhibits a local physical behaviour that is purely two-dimensional and is fully captured by the transition model. Conversely, the standard SA model at the same position presents a much smoother vorticity distribution over the profile surface (Figure 10a,b). Even though separation at the lower surface is still visible in the small layer of negative vorticity produced at the wall, the enhanced diffusivity here induces a coalescence of the two vortices into a single vortex. Furthermore, a single core of vorticity is again visible on the upper surface, which seems, however, to remain more closely attached to the solid surface.
As the position Φ = 0.25 is reached, the profile assumes its maximum downward velocity. By examination of the vorticity contour of the SA-BCM model in Figure 9c,d, it is possible to note a separation with the subsequent incipient roll-up of the boundary layer in the region close to the leading edge at the upper surface, resulting in a thinning of the layer in the immediate vicinity. The two formed vortices, upon their interaction with the solid surface, generate two visible regions of positive-sign vorticity arising from the boundary. This region is found to exhibit striking similarities with the field predicted by the three-dimensional simulation (see [11]), once again justified by the two-dimensional character of the instabilities at their onset. On the contrary, by examination of Figure 10c,d referring to the SA model, it is evident that all the aforementioned details of boundary-layer instability are not adequately captured due to the absence of any modifications accounting for the transition. In particular, the characteristic ripples of the boundary layer edge are completely absent, implying that the local flow acceleration toward the solid surfaces resulting from boundary roll-up is completely sheared out, and the flow almost immediately reattaches at the boundary without any flow organization in a coherent structure. Similar observations are also true for the dynamics of the lower surface boundary layer; although the SA-BCM model predicts the presence of three clockwise-rotating structures, the use of a fully turbulent SA closure does not seem to be able to capture any vorticity ejection into the fluid domain. The higher accuracy resulting from the application of a transition model is made evident upon inspection of the three-dimensional results in [11,12]. The early transition state at the upper surface, with the presence of span-wise homogeneous vortices, seems to be exactly captured by the SA-BCM model. In contrast, the scale-resolving results show that the lower vortex trailing downstream exhibits a significant non-homogeneous spanwise character, suggesting a less accurate description given by the two-dimensional simulation due to the absence of three-dimensional mechanisms of vorticity redistribution.
At Φ = 0.5 , the profile reached the maximum downward displacement. Considering again the referenced three-dimensional results [11,12], the presence of two distinct turbulent regions over the upper surface of the profile separated by a laminar, spanwise homogeneous zone can be seen. By examination of Figure 9e, it can be seen that, corresponding to those regions, the SA-BCM model predicts the presence of two large, highly diffused vorticity regions. The vortex close to the leading edge (see Figure 9f), in particular, arises as a result of the diffusion-induced coalescence of the previously formed vortex into a single structure, attributed to the activation of the source term in the model. Behind this single core region, the boundary layer reattaches laminarly—with the intermittency function (not displayed for this phase) locally assuming a value of zero—before separating again. This second separated region, as can be seen by the presence of the positive sign of the vorticity below the rippled boundary layer edge, gives rise to three spanwise-homogeneous distinct vortices that have also been observed in [11,12]. By examining the corresponding field reported by the SA model in Figure 10e,f, the upper surface single vortex is again observable. However, in this case, it emerges as a single structure after boundary-layer separation at the leading edge. Further downstream to such structure, the boundary layer reattaches and then simply increases its thickness along the downstream direction; no significant instability phenomena nor further separation region can be appreciated in the boundary layer by adopting the SA model.
Finally, at Φ = 0.75 the profile reaches the maximum velocity in its upstroke motion. The three-dimensional results [11,12] show the formation of complex vorticity organization over a large portion of the upper surface as a consequence of interactions between the profile and the previously formed vortex structures. In the SA-BCM model, such interactions seem to be captured in a more adequate manner with respect to the SA model, even though not exactly due to the two-dimensional nature of the simulation. By examination of Figure 9g,h, it is possible to observe the single leading-edge vortex being pushed towards the solid surface by the vertical, motion-induced velocity as it becomes advected downstream. As a consequence of this wall–vortex interaction, a region of positive-sign vorticity is formed and subsequently ejected into the fluid domain. The two-dimensional instantaneous vorticity contours derived via implementation of SA-BCM model appear to be in fairly good agreement with the planar vorticity distribution reported in [11,12]. On the contrary, in the SA case, Figure 10g,h, the situation is significantly different: the vorticity distribution is much smoother with respect to the three-dimensional findings. Moreover, the leading-edge vortex appears to be much more heavily diffused and stretched in the streamwise direction with only a small, positive vorticity region confined close to the solid wall. Further downstream, no visible sign of boundary-layer instability can be identified over the upper surface. The overall effect of the upward motion in the SA case is a compression of the boundary layer towards its bounding solid wall without generation of any clearly detectable instability.
To further demonstrate the improvement brought on by the application of a transition model over the common RANS approach, another detail of vortex generation due to leading-edge separation corresponding to the position Φ = 0.35 is reported in Figure 11, and comparison is drawn again from the spanwise vorticity distribution provided in [11,12] at the same position. Once again, striking similarities can be seen between the two-dimensional SA-BCM results and the fields computed in the three-dimensional simulations reported in the previously mentioned works. In more detail, at this stage of the motion the profile is decelerating towards its maximum downward displacement position; three coherent, clockwise-rotating structures are clearly detectable close to the leading edge, separated by regions of positive-valued vorticity. The last of these formed at previous stages of motion become further stretched and are transported away from the solid boundary upon interaction with the upward-moving fluid of the negative vorticity cores. The interactions between these counter-clockwise-rotating regions then trigger instability mechanisms involving the spanwise direction, as clearly highlighted in [11], inducing a transition to a fully turbulent regime characterized by the collapse of the generated structures into smaller scales. Obviously, such mechanisms cannot be captured by two-dimensional simulations; however, the implementation of a transition model enables description with a sufficient degree of accuracy of the onset and the initial evolution of the boundary-layer spanwise instabilities that eventually evolve into the turbulent structures that are resolved in the three-dimensional problems.
Figure 11b provides the instantaneous value of the intermittency factor γ B C modulating the activation of the eddy viscosity-like production. As the parameter reaches unitary value according to the local flow features, it induces a local increase in diffusivity exactly corresponding to the regions that also appear to develop a fully turbulent behaviour based on the three-dimensional results. In more detail, the model activates, triggering a local enhancement of momentum diffusivity corresponding to the region occupied by the three clockwise-rotating structures close to the leading edge and both on the lower and upper surface extending from the trailing edge up to half-profile length. It could be speculated that the localized increase in diffusivity in the fore region determines the coalescence of the three structures therein present, in turn determining the emergence of the single vortex core observed in Figure 9f.
The implementation of the transition SA-BCM model for this case of the plunging airfoil at a moderate Reynolds number is found to provide an overall improved flow structure prediction when compared to the three-dimensional results with respect to the more commonly used Spalart–Allmaras RANS approach. It should be noted that this comes without any variation in terms of the magnitude of neither the mean nor instantaneous force computations.

3.2.2. Case 2

The second explored case for the asymmetric SD7003 profile has been proposed as a mechanism for stall suppression based on high-frequency small-amplitude oscillations. The mean angle of attack was increased to α 0 = 14 , the Reynolds number to R e = 60,000, and the reduced frequency assumed a value of k = 10 . The vertical amplitude of the oscillations has been set to a value equal to 0.5% of the chord length, i.e., h 0 = 0.005 c . Analogously to the previous case, the corresponding Strouhal number may be derived upon substitution of the above parameters into Equation (5), resulting in S t = 0.03 .
For this configuration, the two models have been firstly tested by considering stationary simulations in a static stall condition, fixing the profile at the constant angle of attack reported above. Table 2 collects the values of C D resulting from different polynomial approximations to provide an example, relative to the steady case, of the space convergence study. The presence of a plateau in the drag coefficient value suggests that an independence of the solution with respect to the adopted space discretization is achieved.
Moving to the steady case results, no noticeable difference could be appreciated between the two models, neither in terms of flow structures nor from the viewpoint of mean drag coefficient. Such a situation is depicted in Figure 12: due to the large separated region, the intermittency function of the SA-BCM model is active over a major portion of the domain such that both cases adopt the same eddy diffusivity equation over a critically large region, providing a justification for the essentially overlapping results.
As far as the time-dependent analysis is concerned, the motion prescribed according to the parametrization given above has been observed to exhibit a significant reduction in the average dimension of the separated region with a parallel diminishing of the drag coefficient up to 40 % in three-dimensional, scale-resolving simulations [11].
Figure 13 reports the mean streamwise velocity field of the two-dimensional case for both the SA and the SA-BCM models. The SA-BCM (Figure 13a) shows a significantly larger reduction in the size of the recirculation region behind the profile when compared with the corresponding effect given by the SA model in Figure 13b. Indeed, the local flow reattachment is a result of the interactions between the stall vortex detaching from the leading edge and the moving surface, so an adequate resolution of such structures appears to be fundamental for capture of an effective stall reduction. When adopting a fully turbulent SA model, instabilities are immediately sheared out by the turbulent stresses so that the mean flow strongly accelerates past the leading edge as a result, exhibiting a neat separation.
Conversely to the case previously explored in Section 3.2.1, an adequate capture of stall reduction phenomena seems to be intimately dependent upon an accurate resolution of the flow structures emerging from leading-edge separation, which is in turn better guaranteed by the application of a transition modification rather than a standard, globally diffusive RANS.
By examining the behaviour of the hydrodynamic forces over multiple cycles in the ( C L , C D ) plane (see Figure 14), it is possible to note that whilst the SA curve exhibits a substantial periodicity, the orbit drawn by the transition model moves within a bounded area in the phase plane.
In particular, the non-periodic SA-BCM solution is characterized by a low-frequency variation in the single-period average drag coefficient between two finite values, both of which are reported in Table 3. Furthermore, it is worth noting that the two C ¯ D limits provided by SA-BCM are significantly smaller than the value predicted by SA and is closer to the three-dimensional reference data in [11].

4. Conclusions

In the present work, two-dimensional high-order DG methods have been applied to the problem of oscillating airfoils at moderate Reynolds numbers, adopting both a standard SA-RANS approach and its algebra-based, transition modification SA-BCM as closure models. Two different profiles, namely NACA0012 and SD7003, have been used under different Strouhal regime and motion conditions to explore the relevant differences between the two models in terms of flow structure and force prediction capabilities.
As far as NACA0012 is concerned, for which a composed flapping motion was considered, the transition model seemed to behave adequately in the entire range of kinematic configurations explored. Moreover, at a low Strouhal number, an overall improved agreement with the experimental results has been observed upon the application of SA-BCM. The latter, in particular, proved to be able to capture flow separation more accurately on the surface of the airfoil for low-frequency oscillations. In contrast, the SA succeeded in this task only when separation was very intense, i.e., for high angles of attack and/or fast dynamics. Similar considerations can be made for the SD7003 profile. The high-frequency plunging configurations showed, upon confrontation with the scale-resolving studies [11,12], that the embedding of a transition model into the standard SA equations guarantees a significant increase in the accuracy of the flow field close to the solid boundary. In particular, at k = 10 the reattachment of the mean flow was captured only by the SA-BCM model, whose prediction of the mean drag was significantly closer to the three-dimensional references with respect to the standard SA.
In conclusion, the distinct transitional aspects of the flow induced by oscillating airfoils at moderate Reynolds numbers, featured by laminar separation followed by turbulent reattachment, make them suitable candidates to test transition modification of standard RANS models. The embedding of an intermittency function γ B C locally modulating the magnitude of the production term within the closure eddy viscosity-like equation allows for an exact resolution of the spanwise-homogeneous region at the leading edge during the early transition state. This, in turn, enables more accurate assessment of the dynamics of the coherent structures arising from such instabilities in comparison to the standard SA. For motion configurations characterized by very low Strouhal numbers, where the interactions of the LEV with the TEV are crucial for correct evaluation of the hydrodynamic forces, the application of the SA-BCM provided results more in line with the experimental data. Contrarily, model performances were found to be almost identical when the dynamics were characterized by combinations of high frequency and large amplitude, where an abrupt separation of the boundary layer occurs. In such cases, γ B C reaches the unit value in a considerably large zone, so SA-BCM and SA basically solve an identical set of equations.

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/fluids8010023/s1, Video S1: Time evolution of the vorticity field in Figure 5a. Video S2: Time evolution of the vorticity field in Figure 5b. Video S3: Time evolution of the vorticity field in Figure 5c. Video S4: Time evolution of the vorticity field in Figure 5d. Video S5: Time evolution of the vorticity field in Figure 5e. Video S6: Time evolution of the vorticity field in Figure 5f. Video S7: Time evolution of the vorticity field in Figure 5g. Video S8: Time evolution of the vorticity field in Figure 5h. Video S9: Time evolution of the vorticity field for Case 1, SA-BCM model (Figure 9). Video S10: Time evolution of the vorticity field for Case 1, SA model (Figure 10). Video S11: Time evolution of the intermittency function field for Case 1 (Figure 11b).

Author Contributions

Conceptualization, A.C.; methodology, A.C.; software, A.C., E.C. and L.A.; validation, A.C., E.C. and L.A.; formal analysis, A.C., E.C. and L.A.; investigation, A.C., E.C. and L.A.; resources, A.C.; data curation, A.C., E.C. and L.A.; writing—original draft preparation, A.C., E.C. and L.A.; writing—review and editing, A.C., E.C. and L.A.; visualization, A.C., E.C. and L.A.; supervision, A.C. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Data Availability Statement

The data presented in this study are available on request from the corresponding author.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:
DGDiscontinuous Galerkin
DNSDirect numerical simulation
DOFDegrees of freedom
GMRESGeneralized minimal residual method
ILESImplicit large eddy simulation
LESLarge eddy simulation
LEVLeading-edge vortex
PDEPartial differential equation
PIVParticle image velocimetry
RANSReynolds-averaged Navier–Stokes
SASpalart–Allmaras
SA-BCMSpalart–Allmaras—Bas-Cakmakcioglu-Mura
TEVTrailing-edge vortex

References

  1. Knoller, A. Die gesetzedes luftwiderstandes. Flug Und Mot. 1909, 3, 1–7. [Google Scholar]
  2. Betz, A. Ein beitrag zur erklaerung segelfluges. Z. Flugtech Mot. 1912, 3, 269–272. [Google Scholar]
  3. Von Kármán, T.; Burgers, J. General Aerodynamic Theory-Perfect Fluids. Aerodyn. Theory 1943, 2, 346–349. [Google Scholar]
  4. Lighthill, S.J. Mathematical Biofluiddynamics; Society for Industrial and Applied Mathematics: Philadelphia, PA, USA, 1975. [Google Scholar] [CrossRef]
  5. Wu, X.; Zhang, X.; Tian, X.; Li, X.; Lu, W. A review on fluid dynamics of flapping foils. Ocean. Eng. 2020, 195, 106712. [Google Scholar] [CrossRef]
  6. Triantafyllou, M.; Triantafyllou, G.; Gopalkrishnan, R. Wake mechanics for thrust generation in oscillating foils. Phys. Fluids Fluid Dyn. 1991, 3, 2835–2837. [Google Scholar] [CrossRef]
  7. Read, D.; Hover, F.; Triantafyllou, M. Forces on oscillating foils for propulsion and maneuvering. J. Fluids Struct. 2003, 17, 163–183. [Google Scholar] [CrossRef]
  8. Schouveiler, L.; Hover, F.; Triantafyllou, M. Performance of flapping foil propulsion. J. Fluids Struct. 2005, 20, 949–959. [Google Scholar] [CrossRef]
  9. Anderson, J.M.; Streitlien, K.; Barrett, D.; Triantafyllou, M.S. Oscillating foils of high propulsive efficiency. J. Fluid Mech. 1998, 360, 41–72. [Google Scholar] [CrossRef] [Green Version]
  10. Ol, M.V.; Reeder, M.; Fredberg, D.; McGowan, G.Z.; Gopalarathnam, A.; Edwards, J.R. Computation vs. Experiment for High-Frequency Low-Reynolds Number Airfoil Plunge. Int. J. Micro Air Veh. 2009, 1, 99–119. [Google Scholar] [CrossRef]
  11. Visbal, M. High-fidelity simulation of transitional flows past a plunging airfoil. AIAA J. 2009, 47, 2685–2697. [Google Scholar] [CrossRef]
  12. Krais, N.; Schnücke, G.; Bolemann, T.; Gassner, G.J. Split form ALE discontinuous Galerkin methods with applications to under-resolved turbulent low-Mach number flows. J. Comput. Phys. 2020, 421, 109726. [Google Scholar] [CrossRef]
  13. Fröhlich, J.; Von Terzi, D. Hybrid LES/RANS methods for the simulation of turbulent flows. Prog. Aerosp. Sci. 2008, 44, 349–377. [Google Scholar] [CrossRef]
  14. Chaouat, B. The state of the art of hybrid RANS/LES modeling for the simulation of turbulent flows. Flow Turbul. Combust. 2017, 99, 279–327. [Google Scholar] [CrossRef] [Green Version]
  15. Yang, W.; Fan, Z.; Deng, X.; Zhao, X. RANS/LES Simulation of Low-Frequency Flow Oscillations on a NACA0012 Airfoil Near Stall. In Proceedings of the 2021 International Symposium on Electrical, Electronics and Information Engineering, Seoul, Republic of Korea, 19–21 February 2021; pp. 62–65. [Google Scholar]
  16. Sanchez-Rocha, M.; Kirtas, M.; Menon, S. Zonal hybrid RANS-LES method for static and oscillating airfoils and wings. In Proceedings of the 44th AIAA Aerospace Sciences Meeting and Exhibit, Reno, Nevada, 9–12 January 2006; p. 1256. [Google Scholar]
  17. Menter, F.R.; Langtry, R.B.; Likki, S.R.; Suzen, Y.B.; Huang, P.G.; Völker, S. A Correlation-Based Transition Model Using Local Variables—Part I: Model Formulation. J. Turbomach. 2004, 128, 413–422. [Google Scholar] [CrossRef]
  18. Langtry, R.B.; Menter, F.R.; Likki, S.R.; Suzen, Y.B.; Huang, P.G.; Völker, S. A Correlation-Based Transition Model Using Local Variables—Part II: Test Cases and Industrial Applications. J. Turbomach. 2004, 128, 423–434. [Google Scholar] [CrossRef]
  19. Walters, D.K.; Leylek, J.H. A New Model for Boundary Layer Transition Using a Single-Point RANS Approach. J. Turbomach. 2004, 126, 193–202. [Google Scholar] [CrossRef]
  20. Carreño Ruiz, M.; D’Ambrosio, D. Validation of the γ-Reθ Transition Model for Airfoils Operating in the Very Low Reynolds Number Regime. Flow Turbul. Combust. 2022, 109, 279–308. [Google Scholar] [CrossRef]
  21. Ashraf, M.; Young, J.; Lai, J. Oscillation frequency and amplitude effects on plunging airfoil propulsion and flow periodicity. AIAA J. 2012, 50, 2308–2324. [Google Scholar] [CrossRef]
  22. Li, X.; Liu, Y.; Kou, J.; Zhang, W. Reduced-order thrust modeling for an efficiently flapping airfoil using system identification method. J. Fluids Struct. 2017, 69, 137–153. [Google Scholar] [CrossRef]
  23. Young, J.; Lai, J.C. Oscillation frequency and amplitude effects on the wake of a plunging airfoil. AIAA J. 2004, 42, 2042–2052. [Google Scholar] [CrossRef]
  24. Wang, Z.; Du, L.; Zhao, J.; Sun, X. Structural response and energy extraction of a fully passive flapping foil. J. Fluids Struct. 2017, 72, 96–113. [Google Scholar] [CrossRef]
  25. Cakmakcioglu, S.; Bas, O.; Mura, R.; Kaynak, U. A Revised One-Equation Transitional Model for External Aerodynamics. In Proceedings of the AIAA Aviation 2020 Forum, Virtual Event, 15–19 June 2020. [Google Scholar] [CrossRef]
  26. Cakmakcioglu, S.C.; Bas, O.; Kaynak, U. A correlation-based algebraic transition model. Proc. Inst. Mech. Eng. Part J. Mech. Eng. Sci. 2017, 232, 3915–3929. [Google Scholar] [CrossRef]
  27. Mura, R.; Cakmakcioglu, S.C. A Revised One-Equation Transitional Model for External Aerodynamics—Part I: Theory, Validation and Base Cases. In Proceedings of the AIAA Aviation 2020 Forum, Virtual Event, 15–19 June 2020. [Google Scholar] [CrossRef]
  28. Bassi, F.; Crivellini, A.; Rebay, S.; Savini, M. Discontinuous Galerkin solution of the Reynolds averaged Navier-Stokes and k-ω turbulence model equations. Comput. Fluids 2005, 34, 507–540. [Google Scholar] [CrossRef]
  29. Bassi, F.; Crivellini, A.; Di Pietro, D.; Rebay, S. An implicit high-order discontinuous Galerkin method for steady and unsteady incompressible flows. Comput. Fluids 2007, 36, 1529–1546. [Google Scholar] [CrossRef]
  30. Crivellini, A.; D’Alessandro, V.; Bassi, F. A Spalart-Allmaras turbulence model implementation in a discontinuous Galerkin solver for incompressible flows. J. Comput. Phys. 2013, 241, 388–415. [Google Scholar] [CrossRef]
  31. Crivellini, A.; D’Alessandro, V.; Bassi, F. High-order discontinuous Galerkin solutions of three-dimensional incompressible RANS equations. Comput. Fluids 2013, 81, 122–133. [Google Scholar] [CrossRef]
  32. Gledhill, I.M.A.; Roohani, H.; Forsberg, K.; Eliasson, P.; Skews, B.W.; Nordström, J. Theoretical treatment of fluid flow for accelerating bodies. Theor. Comput. Fluid Dyn. 2016, 30, 449–467. [Google Scholar] [CrossRef] [Green Version]
  33. Bassi, F.; Botti, L.; Colombo, A.; Crivellini, A.; Franchina, N.; Ghidoni, A. Assessment of a high-order accurate Discontinuous Galerkin method for turbomachinery flows. Int. J. Comput. Fluid Dyn. 2016, 30, 307–328. [Google Scholar] [CrossRef]
  34. Spalart, P.; Allmaras, S. A one-equation turbulence model for aerodynamic flows. In Proceedings of the 30th Aerospace Sciences Meeting and Exhibit, Reno, NV, USA, 6–9 January 1992. [Google Scholar] [CrossRef]
  35. Menter, F.R.; Matyushenko, A.; Lechner, R.; Stabnikov, A.; Garbaruk, A. An Algebraic LCTM Model for Laminar–Turbulent Transition Prediction. Flow Turbul. Combust. 2022, 109, 841–869. [Google Scholar] [CrossRef]
  36. Bassi, F.; Crivellini, A.; Di Pietro, D.A.; Rebay, S. An artificial compressibility flux for the discontinuous Galerkin solution of the incompressible Navier-Stokes equations. J. Comput. Phys. 2006, 218, 794–815. [Google Scholar] [CrossRef]
  37. Lang, J.; Verwer, J. ROS3P—An accurate third-order Rosenbrock solver designed for parabolic problems. BIT 2001, 41, 731–738. [Google Scholar] [CrossRef]
  38. Rumsey, C.L. Apparent transition behavior of widely-used turbulence models. Int. J. Heat Fluid Flow 2007, 28, 1460–1471. [Google Scholar] [CrossRef] [Green Version]
  39. Franciolini, M.; Botti, L.; Colombo, A.; Crivellini, A. p-Multigrid matrix-free discontinuous Galerkin solution strategies for the under-resolved simulation of incompressible turbulent flows. Comput. Fluids 2020, 206, 104558. [Google Scholar] [CrossRef]
  40. Geuzaine, C.; Remacle, J.F. Gmsh: A 3-D finite element mesh generator with built-in pre-and post-processing facilities. Int. J. Numer. Methods Eng. 2009, 79, 1309–1331. [Google Scholar] [CrossRef]
  41. Williamson, C.H.; Roshko, A. Vortex formation in the wake of an oscillating cylinder. J. Fluids Struct. 1988, 2, 355–381. [Google Scholar] [CrossRef]
  42. Koochesfahani, M.M. Vortical patterns in the wake of an oscillating airfoil. AIAA J. 1989, 27, 1200–1205. [Google Scholar] [CrossRef]
  43. Ol, M. Vortical structures in high frequency pitch and plunge at low Reynolds number. In Proceedings of the 37th AIAA Fluid Dynamics Conference and Exhibit, Miami, FL, USA, 25–28 June 2007; p. 4233. [Google Scholar]
  44. Andersen, A.; Bohr, T.; Schnipper, T.; Walther, J.H. Wake structure and thrust generation of a flapping foil in two-dimensional flow. J. Fluid Mech. 2017, 812. [Google Scholar] [CrossRef] [Green Version]
  45. Cimarelli, A.; Franciolini, M.; Crivellini, A. On the kinematics and dynamics parameters governing the flow in oscillating foils. J. Fluids Struct. 2021, 101, 103220. [Google Scholar] [CrossRef]
  46. Guglielmini, L.; Blondeaux, P.; Vittori, G. A simple model of propulsive oscillating foils. Ocean. Eng. 2004, 31, 883–899. [Google Scholar] [CrossRef]
Figure 1. Scheme of the NACA0012 airfoil.
Figure 1. Scheme of the NACA0012 airfoil.
Fluids 08 00023 g001
Figure 2. Mesh of NACA0012 (a) and SD7003 (b).
Figure 2. Mesh of NACA0012 (a) and SD7003 (b).
Fluids 08 00023 g002
Figure 3. Propulsive efficiency and mean thrust coefficient as a function of the Strouhal number for α m a x = 15 (a) and α m a x = 30 (b). Experimental data digitized from [8].
Figure 3. Propulsive efficiency and mean thrust coefficient as a function of the Strouhal number for α m a x = 15 (a) and α m a x = 30 (b). Experimental data digitized from [8].
Fluids 08 00023 g003
Figure 4. ( C L , C D ) phase diagram for α m a x = 15 using transition and fully turbulent model with S t = 0.1 (a) and S t = 0.45 (b). The curves cover a total of two time periods; solid circles (•) identify the starting point of the oscillating cycle and empty circles (∘) indicate the point where vorticity fields in Figure 5 are depicted.
Figure 4. ( C L , C D ) phase diagram for α m a x = 15 using transition and fully turbulent model with S t = 0.1 (a) and S t = 0.45 (b). The curves cover a total of two time periods; solid circles (•) identify the starting point of the oscillating cycle and empty circles (∘) indicate the point where vorticity fields in Figure 5 are depicted.
Fluids 08 00023 g004
Figure 5. Contour plot of the instantaneous vorticity field for α m a x = 15 using the transition and fully turbulent models. The snapshots are taken at Φ 0.183 for S t = 0.1 and at Φ 0.385 for S t = 0.45 . (a) S t = 0.1 , SA-BCM; (b) detail I; (c) S t = 0.1 , SA; (d) detail II; (e) S t = 0.45 , SA-BCM; (f) detail III; (g) S t = 0.45 , SA; (h) detail IV.
Figure 5. Contour plot of the instantaneous vorticity field for α m a x = 15 using the transition and fully turbulent models. The snapshots are taken at Φ 0.183 for S t = 0.1 and at Φ 0.385 for S t = 0.45 . (a) S t = 0.1 , SA-BCM; (b) detail I; (c) S t = 0.1 , SA; (d) detail II; (e) S t = 0.45 , SA-BCM; (f) detail III; (g) S t = 0.45 , SA; (h) detail IV.
Fluids 08 00023 g005aFluids 08 00023 g005b
Figure 6. Contour plot of the instantaneous vorticity field using transition and fully turbulent model for S t = 0.45 . All snapshots are taken at Φ = 0.25 . (a) α m a x = 30 , SA-BCM; (b) Detail I; (c) α m a x = 30 , SA; (d) Detail II; (e) α m a x = 20 , SA-BCM; (f) Detail III; (g) α m a x = 20 , SA; (h) Detail IV.
Figure 6. Contour plot of the instantaneous vorticity field using transition and fully turbulent model for S t = 0.45 . All snapshots are taken at Φ = 0.25 . (a) α m a x = 30 , SA-BCM; (b) Detail I; (c) α m a x = 30 , SA; (d) Detail II; (e) α m a x = 20 , SA-BCM; (f) Detail III; (g) α m a x = 20 , SA; (h) Detail IV.
Fluids 08 00023 g006
Figure 7. Single-period ( C L , C D ) phase diagram at increasing polynomial degree q for the case with α m a x = 20 and S t = 0.45 , using transition (a) and fully turbulent model (b). (i) q = 3 , Δ t = 10 2 c / U ; (ii) q = 4 , Δ t = 5 · 10 3 c / U ; (iii) q = 5 , Δ t = 2.5 · 10 3 c / U ; (iv) q = 6 , Δ t = 2.5 · 10 3 c / U . Solid circles (•) identify the starting point of the oscillating cycle, empty circles (∘) indicate the point where vorticity fields in Figure 6e–h are depicted.
Figure 7. Single-period ( C L , C D ) phase diagram at increasing polynomial degree q for the case with α m a x = 20 and S t = 0.45 , using transition (a) and fully turbulent model (b). (i) q = 3 , Δ t = 10 2 c / U ; (ii) q = 4 , Δ t = 5 · 10 3 c / U ; (iii) q = 5 , Δ t = 2.5 · 10 3 c / U ; (iv) q = 6 , Δ t = 2.5 · 10 3 c / U . Solid circles (•) identify the starting point of the oscillating cycle, empty circles (∘) indicate the point where vorticity fields in Figure 6e–h are depicted.
Fluids 08 00023 g007
Figure 8. (a) Time behaviour of drag and lift coefficients for SA, SA-BCM, and scale-resolving 3D model obtained via image digitization of data in [11]. (b) Phase plane representation of drag and lift coefficients for SA and SA-BCM model.
Figure 8. (a) Time behaviour of drag and lift coefficients for SA, SA-BCM, and scale-resolving 3D model obtained via image digitization of data in [11]. (b) Phase plane representation of drag and lift coefficients for SA and SA-BCM model.
Fluids 08 00023 g008
Figure 9. SA-BCM model. Contour of instantaneous spanwise vorticity. From top to bottom, the vertical positions are presented in the sequence: (a,b) Φ = 0 ; (c,d) Φ = 0.25 ; (e,f) Φ = 0.5 ; (g,h) Φ = 0.75 . The right column reports the details of the instability developing at leading edge.
Figure 9. SA-BCM model. Contour of instantaneous spanwise vorticity. From top to bottom, the vertical positions are presented in the sequence: (a,b) Φ = 0 ; (c,d) Φ = 0.25 ; (e,f) Φ = 0.5 ; (g,h) Φ = 0.75 . The right column reports the details of the instability developing at leading edge.
Fluids 08 00023 g009
Figure 10. SA model. Contour of instantaneous spanwise vorticity. From top to bottom, the vertical positions are presented in the sequence: (a,b) Φ = 0 ; (c,d) Φ = 0.25 ; (e,f) Φ = 0.5 ; (g,h) Φ = 0.75 . The right column reports the details of the instability developing at leading edge.
Figure 10. SA model. Contour of instantaneous spanwise vorticity. From top to bottom, the vertical positions are presented in the sequence: (a,b) Φ = 0 ; (c,d) Φ = 0.25 ; (e,f) Φ = 0.5 ; (g,h) Φ = 0.75 . The right column reports the details of the instability developing at leading edge.
Fluids 08 00023 g010
Figure 11. SA-BCM model. Vertical position Φ = 0.35 . Contour of instantaneous spanwise vorticity (a) and intermittency function γ B C (b).
Figure 11. SA-BCM model. Vertical position Φ = 0.35 . Contour of instantaneous spanwise vorticity (a) and intermittency function γ B C (b).
Fluids 08 00023 g011
Figure 12. Stationary case. Contour of streamwise velocity u 1 . SA-BCM model (a) and SA model (b).
Figure 12. Stationary case. Contour of streamwise velocity u 1 . SA-BCM model (a) and SA model (b).
Fluids 08 00023 g012
Figure 13. Time-dependent case. Contour of the mean streamwise velocity u 1 . SA-BCM model (a) and SA model (b).
Figure 13. Time-dependent case. Contour of the mean streamwise velocity u 1 . SA-BCM model (a) and SA model (b).
Fluids 08 00023 g013
Figure 14. Time-dependent case. Instantaneous phase space representation of hydrodynamic forces over multiple motion cycles (a) and relative enlargement (b).
Figure 14. Time-dependent case. Instantaneous phase space representation of hydrodynamic forces over multiple motion cycles (a) and relative enlargement (b).
Fluids 08 00023 g014
Table 1. Comparison of the mean drag coefficient C ¯ D between the different models, together with experimental and numerical reference data.
Table 1. Comparison of the mean drag coefficient C ¯ D between the different models, together with experimental and numerical reference data.
SASA-BCMVisbal [11]Krais et al. [12]
−0.072−0.073−0.083−0.082
Table 2. Steady case. Space convergence shown by means of the trend of the drag coefficient C D at increasing polynomial degree q.
Table 2. Steady case. Space convergence shown by means of the trend of the drag coefficient C D at increasing polynomial degree q.
qSASA-BCM
10.2100.215
20.1960.204
30.1940.204
40.1940.205
50.1940.205
60.1940.205
Table 3. Comparison of C ¯ D between the different models and experiments for the steady (first row) and high-frequency plunge (second row) cases. The values reported for SA-BCM in the plunging configuration are the minimum and maximum value of the single-period average drag coefficient.
Table 3. Comparison of C ¯ D between the different models and experiments for the steady (first row) and high-frequency plunge (second row) cases. The values reported for SA-BCM in the plunging configuration are the minimum and maximum value of the single-period average drag coefficient.
CaseSASA-BCMVisbal [11]
k = 0 , h 0 / c = 0 0.1930.2050.225
k = 10 , h 0 / c = 0.005 0.1580.133–0.1480.133
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Alberti, L.; Carnevali, E.; Crivellini, A. Assessment of a RANS Transition Model with Flapping Foils at Moderate Reynolds Numbers. Fluids 2023, 8, 23. https://doi.org/10.3390/fluids8010023

AMA Style

Alberti L, Carnevali E, Crivellini A. Assessment of a RANS Transition Model with Flapping Foils at Moderate Reynolds Numbers. Fluids. 2023; 8(1):23. https://doi.org/10.3390/fluids8010023

Chicago/Turabian Style

Alberti, Luca, Emanuele Carnevali, and Andrea Crivellini. 2023. "Assessment of a RANS Transition Model with Flapping Foils at Moderate Reynolds Numbers" Fluids 8, no. 1: 23. https://doi.org/10.3390/fluids8010023

APA Style

Alberti, L., Carnevali, E., & Crivellini, A. (2023). Assessment of a RANS Transition Model with Flapping Foils at Moderate Reynolds Numbers. Fluids, 8(1), 23. https://doi.org/10.3390/fluids8010023

Article Metrics

Back to TopTop