Larabi Abderrahim

Ravelet Florent

Azzam Tarik

Oualli Hamid

Menfoukh Laiche

Bakir Farid
email: farid.bakir@ensam.eu

Numerical Investigation of Parietal Pressure Distribution on NACA0012 Wing Controlled by Micro-cylindrical Rod Arranged in Tandem

Keywords: Laminar separated flow, SST transition model, Passive flow control, Bodies' wake interaction, Laminar separation bubble, Drag reduction, Lift enhancement

The aim of this study is to investigate the influence of disturbed freestream flow by a small cylinder on the laminar separated boundary layer over NACA0012 wing operating at a Reynolds number of Re c = 4.45 × 10 5 . A detailed parametric investigations for the rod are performed using numerical simulations coupled with transition sensitive closure model (γ -Re θ,t ) seeking for the optimal passive control parameters. Firstly, the use of such steady RANS model has been successfully accurate in capturing the separation induced transition on the baseline wing suction surface. Secondly, the rod location was scaled according to the formation length of vortices behind the micro-cylinder for which the aerodynamic loads are very sensitive. The effects of three rod diameter ratios (d/c = 0.67%, 1.33% and 2%) on the laminar separation bubble and aerodynamic performances were examined. It was observed that the qualitative analysis of the flow structures revealed the mechanisms of the control device for the aerofoil performance improvements in which the rod wake exerted considerable effects on LSB size, pressure coefficient and flow streamlines. Particularly, it contributes to eliminate the boundary layer separation with pronounced decrease of 75% by energizing the shear layer over a significant extent, resulting in a mean drag dropping of 73% at 12 • incidence, and a lift enhancement of about 23% at 15 • .

Introduction

The laminar separation bubble is one of major issues met in several industrial applications in particular, those operating at low Reynolds number such us, UAV (Unmanned Air Vehicle), Micro Air Vehicle, race cars and even wind turbine. The boundary layer transition mechanism from laminar to turbulent state of such winged devices is deeply influenced by the aerodynamic characteristics of the wing in terms of boundary layer thickness, transition and separation locations, and wake thickness. In particular, the dynamic behaviour underlying behind these low Reynolds flows showed occurrence of slowly recirculating fluid characterized by a stationary rotating swirl known as a laminar separation bubble identified to be responsible of several negative effects on the overall aerodynamic performances, such as lift decreasing, drag increasing, aircraft stability reducing, vibration and even disturbing noise. This separation bubble is basically developed when the laminar boundary layer cannot overcome the local flow deceleration just coming after its acceleration. Hence, when it bypasses the leading edge of the wing at an incidence; the boundary layer peels off from the wall surface under the influence of an adverse pressure gradient. However, at higher Reynolds numbers this laminar boundary layer grows drastically to a turbulent state due to the presence of large inertial forces, and is therefore better able to overcome this adverse pressure gradient. After that, development of Kelvin-Helmholtz instabilities in this separated shear layer initiates turbulence which yields three-dimensional flow motions. This turbulent mixing is a complex process resulting in the separated streamline reattachment, through interactions between turbulent motions and the wing that incites the reattachment of the separated shear layer to the aerofoil surface, leaving an enclosed region of dead recirculating air. Therefore, understanding the physics involved in this phenomenon and possibility to control the resulting bubble are determinant for efficient engineering design. A review on laminar separation bubbles is presented by [START_REF] Tani | Low-Speed Flows Involving Bubble Separations[END_REF] where effects of Reynolds number on the laminar separation bubble occurrence on an aerofoil are discussed. The bubble is observed to appear for a Reynolds number ranging in between [6×10 4 -6×10 6 ]. [START_REF] Russell | Length and Bursting of Separation Bubbles: A Physical Interpretation of Pressure Distribution[END_REF] introduced a theoretical model characterizing the laminar separation bubble in which he advocated a way of depicting the different critical points describing the LSB throughout the on-surface pressure measurements. The analysis clearly shows that pres-sure distribution on the suction side remains fairly constant when the laminar boundary layer separates, after that the boundary layer reattaches to the surface when shear layer transition is triggered as revealed by drastically pressure increase. [START_REF] Yang | Large-Eddy Simulation of Boundary-Layer Separation and Transition at a Change of Surface Curvature[END_REF] used numerical simulations to compare between averaged over time and instantaneous flow-field near the aerofoil surface close to the separation region obtained by RANS-LES turbulence model. Their main results show typical Kelvin-Helmholtz instabilities in the shear layer is twisting into vortices that force the separated flow to stick towards the wall surface leading to a bubble shape in the time averaged outcome. In experimental studies, [START_REF] Zhang | Hot-Film Measurements of Boundary Layer Transition, Separartion and Reattachment on a Low-Pressure Turbine Airfoil at Low Reynolds Numbers[END_REF] conducted a hot wire investigation on the suction side of a wing at low Reynolds number which showed that a transient boundary layer can reattach to form a bubble and then continues downstream evolving into fully turbulent state. [START_REF] Hu | An experimental Study of the Laminar Flow Separation on a Low-Reynolds-Number Airfoil[END_REF] conducted time averaged and instantaneous PIV experiments to provide information about transient behaviour of the laminar separated flow on low speed NASA GA (W-1) aerofoil at a Reynolds number of 7 × 10 4 based on the chord length. It is found that the bubble appears only when the incidence is in the range [8 • -11.8 • ]. Above this range at around 12 • the pressure gradient becomes to be so significant to induce sudden bubble burst causing the aerofoil complete stalling. The authors reported also that LSB moves upstream towards the leading edge as the incidence increases in such a way its overall length remains unchanged. The laminar portion of this bubble, however, slightly extends in size, and the turbulent portion shrinks with the increasing angle of attack. By the same way, experimental results presented by [START_REF] Hain | Dynamics of Laminar Separation Bubbles at Low-Reynolds-Number Aerofoils[END_REF] in form of PIV measurements at low Reynolds number flow around SD7003 aerofoil considering laminar separation bubble formation, pointed out the existence of a shear layer rolling up into a vortex sheet in the separated boundary layer originating from the presence of Kelvin-Helmholtz instabilities just after triggering of the boundary layer transition. At the meantime, laminar separated flow control over aerofoils at moderate Reynolds numbers has gathered widespread research interest over decades, aiming to enhance their performances. A literature review on flow control shows existence of several techniques to manoeuvre a boundary layer, either in purpose of flow separation postponement so that the flow remains attached to the aerofoil surface, or to achieve lift enhancement, drag decreasing and noise damping. In respect to flow control techniques, these are categorized into two major classes: active and passive techniques. [START_REF] Gad-El-Hak | Flow Control: Passive, Active, and Reactive Flow Management[END_REF] added a reactive class to the above and showed up more specific ways of control depending upon the sought main goal. Many works have attempted to control the unfavourable phenomena that destroy the aerofoils aerodynamic efficiency. Suction and blowing are the most popular adopted active flow control. Although their benefits are obvious in terms of lift enhancement, drag reduction and stall delay as shown in [START_REF] Fatahian | merical Study the Effect of Suction at a Compressible and High Reynolds Number Flow to Control the Flow Separation over NACA2415 Airfoil[END_REF] numerical investigations, the high cost, geometrical complexity and even energy consuming are most of their inherent weakness. Hence, naturally, a good flow control technique needs to be simple to implement, cheap and keeping the overall geometry unchanged. For these reasons, the so-called passive control techniques are the most widely used for simplicity, affordness and do not require any external energy supplying. Several numerical and experimental studies have shown up the effectiveness of such simple devices to control laminar separated flows around different aerofoils. [START_REF] Shan | Numerical Study of Passive and Active Flow Separation Control over a NACA0012 Airfoil[END_REF] focused their numerical investigation on both passive and active on-surface vortex generators used to control the flow separation over NACA0012 aerofoil under conditions of Re c = 10 5 and 6 • incidence angle. The results showed a reduction of the separation zone over a significant extent up to 80% when the passive vortex generator is added. On the other hand, the active vortex generator made the separation zone almost to vanish in the time averaged flow results. However, this control technique did not impart drag reduction neither lift enhancement. [START_REF] Hafien | Passive Separation Control on a Symmetric Airfoil via Elastic-Layer[END_REF] carried out numerical simulations to consider the effect of self-adapting flexible flaps on flow separation over NACA0012 aerofoil at moderate Reynolds number. The results showed a large modification of flow structure when elastic flap deforms under fluid-structure interaction forcing the wake to reduce in size and to lower intensity for shedding vortices. This resulted in a maximum lift enhancement of 69.49% for variable length flaps. Another passive control method was proposed by [START_REF] Vuddagiri | Flow Analysis of Airfoil having Different Cavities on its Suction Surface[END_REF] aiming to investigate the aerodynamic characteristics of low Reynolds number flow over NACA0018 aerofoil. They have conducted numerical investigations using k -ω SST turbulence formulation over a wing model having a cavity on its suction side at various locations where a trapped vortex is acting as passive flow control technique. The geometry of the cavity and its aspect ratio is optimized seeking for higher performances. Elliptical cavity is more promising in enhancing the wing performance, whereas, circular one gave better lift to drag ratio when it is placed near the TE only. Leading edge protuberances is another technique adopted recently for aerofoil laminar flow control. The effectiveness of such vortex generator on the overall performance of S823 wind turbine model was highlighted from the numerical and experimental outcomes conducted by [START_REF] Supreeth | Experimental and Numerical Investigation of the Influence of Leading Edge Tubercles on S823 Airfoil Behavior[END_REF] at 2 × 10 5 Reynolds number. They postulated the main part of tubercles in smoothening and mitigating the stall effects to higher incidences which widen the scope of engineering operations for this modified wing. Moreover, the technique offered an additional increase of 19.38% in the overall aerodynamic efficiency. The control effect of permeable leading edge of NACA008 wing at moderate Reynolds number of 3.45 × 10 5 was illustrated recently by [START_REF] Aldheeb | Aerodynamics of Symmetric Permeable Airfoils and Wings: CFD Simulation[END_REF]. Despite its structure complexity, this method showed a strong influence on the aerodynamic coefficients where an improve of 65% in lift to drag ratio can be achieved. The studies discussed above are mostly based on classical control techniques that introduce geometrical changes to the body itself. Lately, another passive flow control method has gathered considerable attention in many engineering applications. It consists of an off-surface control element located far away from a body interacting with its boundary layer. This concept was first inspired each of these authors ( [START_REF] Strykowski | On the Formation and Suppression of Vortex Shedding at Low Reynolds Numbers[END_REF]; [START_REF] Igarashi | Enhancement of Heat Transfer and Reduction of Drag of a Square Prism in an Air Stream[END_REF]; [START_REF] Igarashi | Drag Reduction of a Square Prism by Flow Control using a Small Rod[END_REF]; [START_REF] Igarashi | Enhancement of Heat Transfer and Reduction of Drag of a Square Prism Arranged Diamond-Shape in an Air Stream[END_REF]) who have suggested a control method of the separated shear layer from bluff bodies by placing a cylindrical rod in tandem arrangement in the upstream crossflow. The idea consisted in forcing the free shear layer to reattach on the body surface. This technique was applied to a cylinder, flat plate, square prism, diamond arranged square prism in order to investigate the flow response to this control method. The outcome recorded a pressure drag reduction of 63%, 30%, 75% and 34% respectively to each bluff body compared to the baseline bodies. Similar observations have been recently reported by [START_REF] Gupta | Suppression of Vortex Shedding in Flow around a Square Cylinder using Cylinder[END_REF] in their DNS investigations to supress the vortex shedding of a square cylinder by small control rod. A remarkable reduction of almost 15% in the mean drag coefficient was achieved. As for laminar separated flow over slender bodies, [START_REF] Zhou | A Numerical Study of a Circular Cylinder in the Wake of an Airfoil[END_REF] examined numerically the wake body interaction between a cylinder placed in the near wake of NACA4412 aerofoil for a Reynolds number of 200 based on the cylindrical body diameter. The effect of the profile incidence, lateral and longitudinal distances between the two bodies in tandem arrangement is focused on the unsteady aerodynamic forces and carefully examined the vortex patterns of the cylinder along with its shedding frequencies. The main outcomes indicate a significant impact on the aerodynamic loadings, vortex patterns and shedding frequencies acting on the cylinder with an upstream streamlined body. In the same manner, [START_REF] Zhang | Aerodynamic Loading on a Cylinder behind an Airfoil[END_REF] experimentally investigated aerodynamic performances of a cylinder placed downstream a NACA4412 aerofoil at chord based low Reynolds number ranged between [1.47 × 10 4 -1.4 × 10 5 ]. They pointed out an abrupt decrease in both lift and drag forces as the side distance is increased. Furthermore, an exhaustive review on effects of wakes interaction and shedding frequencies between the added bluff body and the cambered aerofoil arranged in tandem in the cross flow was given by [START_REF] Bajalan | Experimental Study of Frequency Behavior for a Circular Cylinder behind an Airfoil[END_REF]. The experiments were carried out using the hot-wire anemometry. The findings highlighted the global physical involved phenomena when the aerofoil wake interact with the cylinder body. This mutual interference led to a significant change in Strouhal number values of both obstacles and the vortex patterns depending upon the spacing between them. Also, from an acoustic view, an LES based CFD analysis was performed by [START_REF] Tong | Numerical Analysis of Broadband Noise Reduction with Wavy Leading Edge[END_REF] to investigate the effect of turbulence generated by a rod of 10 mm in diameter place in front of wavy leady edge wing on its broadband noise for Reynolds number of 3.97 × 10 5 . Results showed the fluctuations of lift and drag coefficients are damped out by 65.4% and 71.4% respectively causing a substantial reduction in the mean noise of 9.5 dB. Moreover, this was encouraged recently by [START_REF] Li | Airfoil Noise Reduction using Boundary Layer Control[END_REF] numerical study aiming to supress the far field noise generated by high lift 30P30N aerofoil model. They demonstrated that when the boundary layer con-trol parameters are carefully set, the broadband noise is pronouncedly decreased by 24 dB and the flow field fluctuations are damped. The previously investigated rod scales were mainly above 10% of the aerofoil characteristic length which led to a significant increase in the total drag of these systems. Therefore, a much smaller element scale was selected in case of such slender body's flow control. It was reported from [START_REF] Sefiddashti | Experimental Study of Effects of Circular-Cross-Section Riblets on the Aerodynamic Performance of Risø Airfoil at Transient Flow Regime[END_REF] experimental investigations, that a micro-riblets placed on the suction side of a Riso aerofoil-based wing was used to assess its aerodynamic performances at wide range of angles of attack. They showed a significant drag reduction of 29.7% and 54% for Reynolds numbers of 2.02 × 10 5 and 1.4 × 10 5 respectively, which occurred at 7 • incidence. [START_REF] Wang | Influence of an Off-Surface Small Structure on the Flow Control Effect on Horizontal Axis Wind Turbine at Different Relative Inflow Angles[END_REF] also used the concept of tiny structure to control the separated flow over horizontal axis wind turbine blade. It was found that for the optimum micro-cylinder parameters, the wind turbine power was increased by almost 74% without affecting its stability. Up to this, most of the related works on off-surface vortex generators were focused on either lift enhancement or drag reduction. [START_REF] Michelis | Interaction of an Off-Surface Cylinder with Separated Flow from a Bluff Body Leading Edge[END_REF] have used an off-surface cylindrical rod as a passive device to control the laminar separation bubble on the suction side of an elongated bluff body. The cylinder was placed upstream in the vicinity of the leading edge of the studied body. Their experiments showed that for the optimum control rod configuration, the laminar separation bubble seemed to entirely vanish. With similar experimental approach, [START_REF] Durhasan | Passive Flow Control over an Airfoil by Control Rod at Low Reynolds Number[END_REF] has studied experimentally the laminar separation bubble flow mechanism at pre-stall regime for SD7062 aerofoil controlled by small rod of different diameters at Reynolds number of 3 × 10 4 placed on the suction side at various chordwise locations. The PIV results revealed significant reduction in the size of the laminar separation bubble for proper rod diameter and location, along with decrease in height of the boundary layer by 22% which is coherent with the correlated reduction of 34% in the eddies' turbulent energy in the flow. Considering the mentioned studies, there was an attempt to either suppress the vortex shedding of bluff bodies or improve aerofoils general aerodynamic performance in terms of loading enhancement or noise damping with large scale rod dimensions control element. In this study, a micro cylinder off-surface vortex generator is used to assess the feasibility of such passive technique to control the laminar separation bubble (LSB) over the suction surface of NACA0012 aerofoil at moderate Reynolds number. It aims to establish the optimum size and location of the tiny cylinder to improve its aerodynamic performance. In addition, except for experimental approaches, an issue of using fully turbulent classical RANS models (Spallart-Almaras, k -, k -ω SST) to study flows involving laminar-turbulent separation has emerged as they require LES or experimental backup to identify the separation point. Therefore, detailed parametric investigations are carried out using RANS based simulations with transition sensitive closure model formulation known as Transition SST (γ -Re θ,t ) to extract the optimal configuration related to the rod for laminar separated flow control and emphasize the importance of including transitional effect in the numerical model.

Numerical Method

Case Geometry

In this study, three-dimensional flow over a NACA0012 wing model controlled by a micro-cylindrical rod is considered to assess its effect on the overall aerodynamic performances of this aerofoil. The global dimensions of the case model along with geometrical arrangement adopted for the passive control element are sketched in Fig. 1. A circular cross-section cylinder of various diameters (d/c = 0.67%, 1.33% and 2%) is chosen to act as an off-surface vortex generator and placed at a distance, L upstream from the rod centreline in the perpendicular direction to NACA0012 aerofoil. First, preliminary experiments are performed in this work where the spacing between the two components of the system is scaled according to the corresponding formation length, l f of the vortices behind the tiny cylinder at flow conditions of 43.3 ms -1 freestream velocity and Reynolds number based on the aerofoil chord length, Re c = 4.45 × 10 5 . In accordance with [START_REF] Norberg | LDV-Measurements in the Near Wake of a Circular Cylinder[END_REF] experimental investigation, for a given Reynolds number of Re d = 5931 based on the micro cylinder diameter, the formation length is roughly in the order of 2 × d (Fig. 2). Based on this, the method of fixing rod location and its effect on the wing gain in terms of maximum lift to drag ratio is reported and rod diameter scaling effect on LSB is explored. Four values for the distance, L are explored aiming to determine the optimum value. The spacings are set to be equal to 1 × d, 2 × d, 3 × d and 4 × d which correspond to location states of before, in, after and faraway the formation zone respectively. The sensitivity of varying L on the total aerodynamic loads is explicitly shown in Table 1.

L 0 × d 1 × d 2 × d 3 × d 4 × d c L c D max c L0 c D 0 max
1.00 0.95 1.04 1.12 0.91

For cases where the considered spacing is halved or doubled according to the formation length, a loss of 10% in lift to drag ratio is observed compared to the baseline wing which reflect a severe deterioration of its aerodynamic performances. Thus, these two situations are not preferred as either the interaction is very weak between the rod aerofoil system or higher drag that corresponds to wider system's wake. In one hand, as the gap is close to the formation length, l f , the system recovers its aerodynamic performance with slight increase in the maximum lift to drag ratio. In the other hand, the advantage of setting the micro-cylinder at position of L = 3 × d is demonstrated. As in such circumstances, an important gain in total lift to drag ratio is clearly achieved. Thus, after such comprehensive consideration the spacing adopted for the remaining study is 3 × d as higher performance of the wing can be achieved compared to other sets.

Computational Domain

To examine the sensitivity of the control element diameter scaling on the laminar separated flow and the aerodynamic performance of the cosidered wing model, detailed numerical investigations are presented in this paper. The three dimensions computational domain used to carry out this numerical study is illustrated in Fig. 3. The coordinate system origin is placed at the leading edge of the wing. The overall domain dimensions are 6.3 m long, 3.6 m height and 0.465 m wide. The inlet boundary of the main domain is positioned 12 times the wing's chord length upstream from the centre point of the main aerofoil, and the rear boundary is extended giving a longitudinal length of 30 × c, where c is the wing's chord length. These domain extensions have been chosen in order to maintain a low blockage level and to ensure as little influence of all domain boundaries as possible which is demonstrated in [START_REF] Mitchell | Aerodynamic Characteristics of a Single Airfoil for Vertical Axis Wind Turbine Blades and Performance Prediction of Wind Turbines[END_REF] numerical study. A sensitivity test has been performed to optimize calculations in terms of time, memory and cost. The smallest size domain giving accuracy of less than 1% in terms of aerodynamic forces; with respect to each chosen downstream and lateral domain lengths; has a total height of 24 × c and total length of 42 × c with 2.7% blockage. The spanwise length, L z is carefully chosen as a short spanwise width will fail to deliver results comparable with data from wind tunnels especially at high angles of attack. Thus, domain size dependence is inevitable to capture the 3D flow structures and ensure reliable results. For such analysis, a pair of symmetrical boundary conditions is set for both spanwise sides. Figure 4 shows that for spanwise length less than 2 × c calculations give higher prediction for both coefficients. When further increasing this length, the variation tendency for both lift and drag seem to level off. Hence, a spanwise width of 3 × c is selected for the 3D simulations. In addition, these domain characteristics are in good agreement when compared to literature. Based on the relevant numerical investigations carried out by [START_REF] Lei | Spanwise Length Effects on Modelling of Flow over a Circular Cylinder[END_REF] and [START_REF] Aguedal | 3D Vortex Structure Investigation using Large Eddy Simulation of Flow around a Rotary Oscillating Circular Cylinder[END_REF], in which the effect of spanwise length of a computational domain on the accuracy of predicting the 3D flow structures past a circular cylinder was analysed. They showed that a minimum spanwise length of π times the characteristic length is required to predict a formation length that is in better agreement with experimental value. Thus, as this parameter plays an important role in rod distance parametric investigation, the adopted spanwise length for the current domain lies in the appropriate range.

Mesh Generation

The adopted grid topology for the considered computational domain is generated using ANSYS R -Meshing commercial software as illustrated in Fig. 5. In order to construct a high-quality structured mesh, the domain is divided into several blocks allowing finer mesh in regions of susceptible eddies formation and coarse grid in far field regions. This strategy is performed using interface boundary conditions between different blocks. Furthermore, to reduce the overall skewness and to ensure better cells orthogonality for the near wall grid elements as density is increased for the following mesh dependence analysis, an O-grid type topology is adopted and centred around both the wing and the cylindrical rod. This approach gives a smooth mesh transition in the leeward zone of the rod to the leading edge of the wing and allows a high resolution of the boundary layer while keeping relatively a low cell counts compared to C-grid type topology used in former similar studies which showed a reduction of about 50% in mesh size. The structured quadrilateral elements are chosen for the whole computational domain since they provide a high flexibility to control the cell counts with low memory consumption, accurate flow resolution and good flow alignment with respect to the mesh cells thus the overall rate of solution convergence is increased. Refined quadrilateral meshes are applied to critical flow zones where most important phenomena affecting the flow structure around the wing are recorded. These zones are mainly the leeward side of the rod where the vortex shedding interacts with the leading edge of the wing ensuring thus minimization of the numerical viscosity effect responsible of artificial dissipation of the shedded vortices. In addition, a particular attention is paid to the mesh refinement at the aerofoil suction side to capture accurately flow features (Laminar Separation Bubble) involved in vicinity of the wing aerofoil including leading and trailing edges. In order to assess the significant impact of grid density on accuracy and convergence rate of the numerical solution, three different grid sizes are investigated for the mesh dependency analysis (Table 2). The coarse grid (MESH-1) contains 200 points on the suction side of the aerofoil, the fine mesh (MESH-2) with 450 grid points and the refined grid (MESH-3) has 640 points. It is noticed throughout sensitivity assessment for lift coefficient force that a deviation of about 4.59% is recorded for the coarse (MESH-1) and fine (MESH-2) meshes. While for the refined mesh (MESH-3), this difference in c L drops to reach a value of roughly 0.56%. As a result, MESH-2 is retained and found to be suited to achieve accurately enough results for this numerical investigation. Another important grid independence was found to be in the exposure of turbulence models to y + values. Values of y + less than unity found to be necessary in order to capture transitional features of the flow on the suction side of the aerofoil by accurately modelling the boundary layer. A target value of 0.7 based on the freestream inlet velocity is applied to the entire wing surface and the cylindrical rod during pre-processing simulation case. After that, y+ distribution along the upper and lower surfaces of the aerofoil and around the cylinder circumference is checked in post-processing phase for all simulations and showed to be less than unity. A total of 15 inflation layers were chosen with a stretching ratio of 18% at the wall surface giving a prism layer thickness large enough to fit development of the boundary layer along with the separation bubble that may occur. As a result of this parametric strategy, a reasonable quality and sufficiently fine structured mesh is with solution-mesh independency condition satisfied.

Turbulence Model

The conventional turbulence models available for RANS formulation such as one equation Spalart-Allmaras model and two equations k -and k -ω models, do not include laminar effects especially for low Reynolds number flows, thus considering fully turbulence models in such studies is invalid assumption. Even though, DNS and LES methods solve turbulence with sufficient precision so that laminar and transition effects can be predicted, grid density and time resolution requirements are prohibitive and probably inappropriate to be used for industrial design purposes. This led to an attempt to include laminar transitional effects using damping functions in near wall regions to gain relatively good prediction of transition by diffusing freestream turbulence into the boundary layer. In addition, this strategy offers a good balance between grid density requirement and flow complexity. This model known as Transition Shear Stress Transport (γ -Re θ,t ) introduced by [START_REF] Menter | A Correlation-Based Transition Model using Local Variables-Part I: Model Formulation[END_REF] for low Reynolds number approach introduces a so-called viscosity limiter that reduces prediction of the wall shear stress thus it is more likely to get separation over smooth surface and results will give better agreement with experiments of separated flows. Several RANS turbulence models are used to investigate laminar separated flows at low Reynolds numbers. It is noticed that Transition SST model provided good prediction of pre and post stall behaviour for aerodynamic loads ( [START_REF] Benim | Prediction of Roughness Effects on Wind Turbine Aerodynamics[END_REF]) and is in good agreement with correlating comparative wind tunnel experiments ( [START_REF] Steed | High Lift CFD Simulations with an SST-Based Predictive Laminar to Turbulent Transition Model[END_REF]; [START_REF] Roberts | Boundary-Layer Transition on Wings in Ground Effect[END_REF]). Furthermore, this latter model proved its reliability in capturing flow interaction features between bodies and its high revelation of parietal pressure oscillations. According to what has been discussed so far, it seems that Transition SST turbulence model is capable of capturing accurately laminar separation with transition of the boundary layer on wings operating at low Reynolds flow regime, allowing to solve details inside the boundary layer when blending effectively between the standard k -ω in near wall regions and k -model in the far field. Hence, Transition SST model is adopted for all simulations which eventually increases the CPU time to about 22% whilst transition sensitive SST model formulation is used to model the laminar separation bubble near the leading edge when boundary layer transits form laminar to turbulent state, because this LSB is a dominant phenomenon when characterizing pressure distribution, on-surface flow structure and wall shear stress over the wing and therefore the overall aerodynamic loads of the wing. However, even though SST k -ω is inefficient in correlating a good prediction of flow features at moderate Reynolds, it shows a comparative correlation with wind tunnel experiments at very high angles of attack.

∂ (ρk) ∂t + ∇ • (ρUk) = ∇ • µ + µ t σ k ∇k (1)
+ P k -D k ∂ (ρω) ∂t +∇ • (ρUω) = ∇ • µ + µ t σ k ∇ω + γ ν t P k (2) -βρω 2 + 2 (1 -F 1 )
ρσ ω2 ω ∇k : ∇ω Where,

∇k : ∇ω = ∂k ∂x j ∂ω ∂x j P k → P k γ 0 ≤ γ ≤ 1 (3) D k → D k min [max (γ, 0.1) , 1.0] (4)
F 1 , new = max (F 1 , F 3 ) (5)
Where, F 3 = e -(Rey/120) 3 and Re y = ρy √ k µ

The γ -Re θ,t turbulence model still solves the transport equations for both k and ω quantities as k -ω SST does, which are otherwise identical with an exception for production, P k and dissipation, D k terms in kinetic energy (Eq.( 1)) along with the blending function F 1 in specific dissipation rate (Eq.( 2)) that controls the use of ω near the walls and far away. The production term, P k is correlated with a turbulence intermittency, γ (Eq.( 3)) bounded between values of 0 and 1, indicating physically the state of the flow locally by accounting for the percentage of time corresponding to the presence of turbulent fluctuations in the boundary layer. Where a value of 0 indicates that the flow is locally laminar hence close to the leading edge of the aerofoil and when γ saturates to a value of 1 the flow is fully turbulent. Notice that freestream flow contains turbulent kinetic energy and specific dissipation rate, thus γ should have a value of 1 at inlet boundary condition and freestream so that it will be represented by the fully k -ω SST model. Regarding the dissipation term, Dk it is limited by a min max limiter (Eq.( 4)) calibrated to ensure that its value does not drop below 10% of its fully turbulent value when γ = 0, because even when the boundary layer is laminar with no production of turbulent kinetic energy, if is there any turbulence the wall still dissipate it by damping these turbulent fluctuations of the flow. In other hand, numerical tweak is used to correct the blinding function, F 1 to prevent it to take a value of 0 in the laminar boundary layer near the wall where normally takes a value of 1 allowing the Transition SST model to take its k -ω form of k -ω SST turbulence formulation rather than k -which is ineffective in the near wall regions particularly when the values of y + are very small. Hence, the blending function is limited; to prevent it to switch accidentally to 0 in the laminar region; by a new term defined by an exponential decay of a Reynolds number like term (Eq.( 5)) governed by a square-root of turbulent kinetic energy which is very small in the laminar portion of the boundary layer making this new term to tend to 1.

∂ (ργ) ∂t + ∇ • (ρUγ) = ∇ • µ + µ t σ γ ∇γ (6)
+ P γ -D γ P γ,1 = F length C a1 ρ S (γ F onset ) 0.5 (1 -C e1 γ) (7) ∂ ρ Re θ,t ∂t +∇ • ρU Re θ,t = (8)
∇ • µ + µ t σ θ,t ∇ Re θ,t + P θ,t P θ,t = 0.03 ρ t Re θ,t -Re θ,t (1 -F θ,t ) (9)
In order to compute the values of γ in the flow field an additional standard transport equation for intermittency is solved (Eq.( 6)). The production terms, P γ is used to control the length of the transition region of the boundary layer by producing turbulent fluctuations and pushing γ to saturates at a value of 1 in the full turbulent region. This term takes the form of Eq.( 7) characterized by two main terms F onset and F length , where the former is a switching function for the production of intermittency and the latter controls the rate of this production happens to be when F length is large γ is produced quickly and saturates to 1 rapidly which shortens the length transition region. To assess the F onset , the Transition SST model introduced another transported variable, Re θ,t calculated everywhere in the computational domain (Eq.( 8)). This variable connects the empirical correlations to the onset criterion of intermittency by calculating two local values Re θ,t and Re θ,c where transition occur and fluctuations begin to take place respectively. The source term, P θ,t of Eq. ( 8) is included to force Re θ,t to take its experimental value (proprietary of ANSYS R -Fluent 18.2 CFD code) specified at the inlet except near the wall where it is turned off allowing the freestream value of Re θ,t to convect and diffuse toward the wall by introduction a blending function F θ,t (Eq.( 9)). Once Re θ,t is computed F length and F onset are evaluated using algebraic relationships and when they are beyond certain threshold the source term, P γ is switched on, hence production of intermittency γ will start for each cell ( [START_REF] Langtry | Correlation-Based Transition Modeling for Unstructured Parallelized Computational Fluid Dynamics Codes[END_REF]; [START_REF] Malan | Calibrating the γ -Re θ Transition Model for Commercial CFD[END_REF]; [START_REF] Colonia | Assessment and Calibration of the γ-Equation Transtion Model at Low Mach[END_REF]).

Boundary Conditions

In addition to the domain sensitivity analysis conducted to set right dimensions of upstream, downstream, top and bottom boundaries of the computational domain, an appropriate implementation of boundary conditions needs to be carefully considered. This is achieved by splitting the domain boundaries into three distinct regions, velocity inlet, pressure outlet downstream and symmetry side boundaries. A velocity inlet type was specified as Dirichlet boundary condition at inlet boundaries in such a way that freestream velocity is set to maintain a chord-based Reynolds number at 4.45 × 10 5 . The used method for freestream turbulence is by setting values of turbulent kinetic energy and specific dissipation rate at the inlet boundary condition, unless otherwise stated, as k = 2.53 m 2 s -2 , ω = 151.55 s -1 and an intermittency factor of 0.85 as dependence of flow features on γ values at the inlet was observed. These quantities are calculated based on empirical correlations related to SST k -ω turbulence model. The downstream domain boundary was prescribed as pressure-outlet boundary condition and set to a value of zero-gauge pressure as operating conditions of the working fluid are specified according to standard atmospheric pressure of 101.325 kP a and temperature of 288.15 K. In addition, an estimation of turbulent length scale of 10% the characteristic length is set at the pressureoutlet along with 5% backflow turbulence intensity. A stationary wall condition is adopted for both aerofoil and rod surfaces with no-slip boundary condition, hence, no turbulence production from the wall itself to the freestream flow. The two side boundaries of the computational domain are assigned to symmetry boundary condition which enforces to null the normal component of all flow variables gradient by acting as zero shear slip wall. This condition seems to be satisfactory as the computational domain is large enough to maintain these far boundaries at constant static pressure (Fig. 3).

Solver and Numerical Schemes

[27] have shown that the laminar separation bubble does not involve any vortex shedding via their time-dependent numerical investigation of the laminar separated flow subjected to an adverse pressure gradient. As a motivation, three dimensional, incompressible, steady state, pressure based and segregated double precision numerical solver is adopted to conduct present simulations. All Steady RANS CFD runs were undertaken using the industry standard commercial CFD package known as ANSYS R -Fluent, Release 18.2. The code solves the Reynolds Averaged Navier Stokes equations in a finite volume environment with Message Passing Interface library parallelization. The Transition SST (γ -Re θ,t ) four (4) equations turbulence closure model is adopted to model the Reynolds stresses term. The choice of this model is motivated by its ability to capture up to a good extent the separation bubble that occurs on the suction side of the wing which appeared to be harmful for its aerodynamic performances and also, has proven reliability in flows involving wake interactions. Pressure-Velocity field coupling is ensured using SIMPLE algorithm advantageous for problems requiring boundary layer high resolution. Spatial discretization of pressure term is treated by a standard interpolation scheme. Moreover, momentum diffusive terms, turbulent kinetic energy and specific dissipation rate equations are second-order central difference special discretization scheme and convection terms are treated using second-order upwind discretization scheme in space. Finally, the cell gradients of all variables are computed using weighted least square cell-based scheme. The initial flow field for all simulations is set according to predicted freestream inlet conditions. The different cases are computed for 50000 iterations where scaled residuals of all equations along with lift and drag coefficients of the wing are monitored to ensure adequate convergence. The overall computational time for every simulation to be converged was achieved in about 144 hours (wall clock) that led to reduction in the residuals of O(6) (Fig. 6). Afterwards, to ensure the convergence accuracy, both lift and drag coefficients are monitored to the extent in which they are fully converged as demonstrated in Fig. 7 where solutions for all cases are continued until there is no variation in C L and C D which characterized by the horizontal asymptotic looks like.

Results and Discussions

Simulations of the flow around the wing at various angles of attack ranging between [0 • -24 • ] for a chord-based Reynolds number of 4.45 × 10 5 . are carried out for both baseline and controlled configurations in order to understand how wing surface pressure are altered when a passive control method is introduced to enhance the overall wing aerodynamic performances at pre-stall and stall regimes.

Validation

For reason of results confidence, it is imperative to compare the obtained results from CFD simulations of the baseline wing with Wind Tunnel experiments to validate the numerical model adopted herein. Tests are carried out in an open jet closed circuit wind tunnel having circular test section of 600 mm in diameter. The wing model is NACA0012 aerofoil based with 150 mm chord length and 465 mm span giving a maximum blockage ratio of 6.4%. The forces acting on the wing under the stream are measured by aerodynamic balance composed of three masts D#1, D#2 and D#3. The measurements range of each one is respectively [-100 N -+100 N ] for both D#1 and D#2 and [-50 N -+50 N ] for the last one. A mechanism for angle of attack variation is also incorporated with the

Mast-D#1

Mast-D#2

Mast-D#3

balance allowing α to be set between [-20 • -+40 • ] (Fig. 8).

The experiments for the baseline wing configuration used for CFD validation are obtained for a freestream velocity of 43.32 ms -1 , that corresponds to a Reynolds number based on the chord length of 4.45 × 10 5 . The turbulence intensity for these flow conditions is 1%. Table 3 shows the error in the aerodynamic forces predicted by current numerical investigation when compared to CFD simulations of [START_REF] Huang | Numerical Study of Blowing and Suction Control Mechanism on NACA0012 Airfoil[END_REF] and the wind tunnel results. It is noteworthy that deviations remain less than 8.28% in all instances except at 14 • incidence where the wing stalls and making this region quite difficult to solve with turbulence modelling formulation of Navier-Stokes equations. Although simulations overpredict lift force coefficient the fact that the error is relatively consistent for each case, highlights that the overall trends shown by the numerical results are coherent with those occurring in the wind tunnel. The comparison between CFD, experimental and literature data are shown in Fig. 11. [START_REF] Jacobs | Airfoil Section Characteristics as Affected by Variations of the Reynolds Number[END_REF], [START_REF] Critzos | Aerodynamic Characteristics of a NACA0012 Airfoil Section at Angles of Attack from 0deg to 180deg[END_REF] have presented wind tunnel experiments in which flow characteristics over NACA0012 aerofoil is deeply investigated at various Reynolds numbers and wide range angles of attack. The obtained results from numerical and experimental approach adopted in the current work are validate against these produced experimental data. As it is seen, for angles of attack less than 12 • both CFD and experimental outcomes are in good agreement with the selected experimental data. The highest recorded error is less than 4% for lift coefficient and in the order of 5% in case of drag coefficient. Whereas, after stall angle, the numerical predicted results show a significant deviation from the current experimental data especially for lift coefficient, which can be attributed to the effect of close walls in experiments. On the other hand, lift coefficient from CFD gave satisfactory results as it is in better agreement with [START_REF] Jacobs | Airfoil Section Characteristics as Affected by Variations of the Reynolds Number[END_REF] experimental data.

In addition, statistical analysis between numerical predic- tions and experimental data is performed using Pearson correlation criterion. It consists of assigning a quality level for given data set, by measuring the linear correlation strength between two sets of data. The related coefficient, r is given by Eq.( 10) as follow :

r X,Y = cov(X, Y ) σ X × σ Y (10)
Table 4 summarises the values of the strength for possible linear correlation between the present obtained experimental as well as numerical prediction results and the experimental data taken from literature in terms of lift and drag coefficients. The main interpretation behind, r is that any value ranged between [0.9 -1] indicates the existence of a very strong correlation between the data sets. Thus, overall, CFD results and experimental measurements gave satisfactory prediction for both lift and drag coefficients and are comparative to experiments taken from literature. ]. Thereafter, it increases slightly in the meantime as the angle of attack is further changed possibly because of premature separation of the laminar boundary layer, followed by a significant rise for α = 14 • indicating presence of turbulent stall of the boundary layer. By further increasing the angle of attack, the drag coefficient continues growing till a maximum value is attained at incidence 24 • . The streamline close to the wing surface is presented in Fig. 14 to make a comparison between various angles of attack and the presence of flow separation on the upper surface of the aerofoil. It is believed that the nonlinearity in the lift curve slope is linked to presence of laminar boundary layer separation occurring at this low Reynolds number flow (Re c = 4.45×10 5 ). For an angle of attack between 0 • to 8 • , a separation bubble appears on the aerofoil suction side that travels from the trailing edge to the leading edge as the incidence is increased, leading effectively to a change in location of the transition line from laminar to

Effect of Passive Control on LSB

Three-dimensional development of the LSB on the wing upper surface is inferred from average pressure distribution analysis on suction side presented in Fig. 15. Showing mean pressure rapid increasing along the chordwise direction. It can be seen from colour contrasts of pressure . Comparing the two configurations side to side, baseline case and controlled case, it seems that the rod has significant influence on dimensions of the laminar separation bubble occurring on the suction side as a result of premature laminar boundary layer separation. Both regions, laminar and transitional, are influenced by presence of the rod wake. A significant reduction of approximately 50% in the laminar region is depicted by the three-dimensional average pressure distribution on the wing surface which yielding to a considerable enlargement of the transitional zone by approximately 30%. Thus, this vortex generator placed upstream the wing showed that LSB is extremely sensitive to the interaction between the rod wake and the aerofoil leading to a bubble size reduction that is considered to be harmful for the wing overall aerodynamic performances. Therefore, a promising enhancement in aerodynamic characteristics will certainly contribute to improve efficiency of such wings operating at low Reynolds numbers. 17 demonstrates both friction lines and average pressure distribution on the wing upper surface. When incidence is null (α = 0 • ), friction lines on the wall surface are uniform and flows in the same direction as that of the main incoming flow, starting from the leading edge towards the trailing edge. As the angle of attack is increased (α < 11 • ), friction lines begin to lose their uniformity and smoothness behaviour. An instability is triggered on both ends of the trailing edge, indicating presence of a partial separation of the boundary layer, developing progressively from the trailing edge towards the leading edge. Furthermore, the average pressure gradually increases in the chordwise direction starting from the leading edge. It begins with a zone of depression at the leading-edge level indicated by a uniform and homogeneous coloured zone and keeps growing to a high-pressure region at the trailing edge. For the baseline configuration at α = 11 • , the depression zone is narrow on both ends of the leading edge, because the instability generated between endplates and aerofoil boundary layers mutually interact. This is much clearer through the friction lines distribution showing a pair of counter-rotating vortices swirling on the wing suction side. Consequently, partial turbulent stall of the main wing is triggered close to these regions. In case of controlled wing configuration, it shows, however, important changes noticeable on the mean pressure field as well as on the friction lines (Fig. 17). Mainly, the depression region near the leading edge is getting uniformized in the spanwise direction starting from middle towards endplates. Another important outcome is revealed by new friction line patterns resulting from the control and consisting in the fact that the pair of vortices near endplates shrink in size making the incident flow more dominate on the wing surface. As a result, a reduction of about 56% on the overall drag force is reported for the controlled case. As the angle of attack is increased the instability amplifies more and more to form two large counter-rotating vortices covering a large part of the wing upper surface. A backflow is generated as a consequence of these two vortices. This reveals that energy of the incoming flow is more dominant than that of the backflow. As a result of this dominance the incident flow imposed its presence mainly in the midspan zone, hence, the boundary layer is not completely separated from the wing surface except in the region where the vortices are present. For α = 12 • , direction of friction lines on the midspan zone shows that the backflow becomes dominant and imposes its presence with respect to the incoming flow. The two counter-rotating vortices still exist, but larger in dimensions and intensities when compared to those observed for incidence α = 11 • . Therefore, secondary instabilities appear as a consequence of this interaction between the two flows of opposite sense and the significant energy of this backflow, Fig. 18.In this case, the boundary layer is completely separated from the wall. This irregularity behind appearance of this instability disturbs the uniform distribution of the average pressure in both ends of the wing close to the leading edge. However, when flow control is applied, the average pressure distribution becomes uniform and homogeneous and the disturbed zones near the endplates vanish. In addition, this depression portion grows in size compared to previous angles. Moreover, presence of the rod has a significant impact on the flow developing on the wing suction side. From the Fig. 18 it is easily noticeable that secondary instabilities reported before completely disappear and the incoming flow gained energy from the feeding vortex generator to become dominant over a large portion of the wing because it overpassed the early backflow seen without control. As a consequence, the two counter-rotating vortices died out and straitened in size resulting a reduction in the drag force by 73% including the rod contribution. For in- cidence greater than 12 • , the uniform distribution of the mean pressure over the entire upper surface is completely disturbed for the baseline configuration, Fig. 19. Consequently, the flow within the boundary layer is completely unstable. The intensity of these instabilities increases with the angle of attack causing the so-called turbulent stall to launch. As α increases, these instabilities are intensified and the boundary layer separation is generalized to the entire wing surface. For α = 13 • the flow seems to behave in the same way as for 12 • incidence without control. This is clearly shown by friction lines pattern, where the dominance of the backflow resulting from instabilities generated inside the boundary layer is well represented leading to secondary irregularities as a result of the interaction between the incident and back flows. In this case, the boundary layer is completely detached from the surface and this incident angle corresponds to the turbulent and massive stall angles of the flow. In the other hand, as the control technique is applied, the mean pressure distribution on the suction side of the wing exhibits a different behaviour in comparison to the baseline case. As stated earlier, the mean pressure distribution is random due to partial separation of the boundary layer near endplates caused by the primary and secondary vortices that gain intensity to cover an important portion of the wing surface as illustrated in friction lines patterns showing dominance of the backflow over the incoming one. Thereby, presence of the vortex generator brought the depression zone near the leading edge from its randomness attitude to its uniform and homogenous shape similar to what was observed at lower angles of attack when control is absent. It is worthwhile noting from pressure contours that flow characteristics are delayed by 2 • from the baseline configuration with flow control applied. In other words, the stall is postponed when the wing is operating in the rod wake. In fact, looking at the new pattern of friction lines after low control is set, what have been seen in average pressure distribution is emphasized throughout the assessment of the flow before and after the control. As a matter of fact, the wake produced by the rod evolves into shedded eddies interacting with the boundary layer developing on the wing upper surface in a way that they feed energy to the flow close to the wall from free-stream. Premature separated boundary layer is reinforced to overcome pressure gradients that is subjected to and stick again to the wall surface. Consequently, this phenomenon of energizing the upcoming flow is established by defeating the dominance of the backflow spotted in the baseline case, so that the flow recovered its chordwise direction mostly in the midspan portion of the wing. Although the pair of vortices remain near the endplates a part of consequence of three-dimensionalities attitude of the flow outcomes of interferences between the two boundary layers developing on either endplates and wing body, the secondary instabilities are completely vanished for the controlled case. To sum up this qualitative analysis that provides main differences on the flow structures before and after control, it is important to link these essential findings to their impact on aerodynamic loads specifically lift and drag forces. It is found that actual aerodynamic coefficients are highly sensitive to flow control in an extent of a lift enhancement of 17% to gather with a maximum drag reduction of 49% which is promising in a perspective of fineness improvement. depicts mean surface pressure distribution on the suction side generated by the wing operating in the low Reynolds number condition and at 15 • incidence angle. For the noncontrolled case the surface pressures appear very instable, as expected, whilst at such high angle of attack a presence of an intense turbulent agitation of the flow within the boundary layer is naturally understandable as turbulent stall of the boundary layer has taken place. This is shown on the friction lines illustration (bottom right-hand column), which reveal in the low-pressure regions near endplates formation of a pair of counter rotating vortices due to circulation of flow from higher to lower pressure regions. Also, when looking at the evolution way of these lines going from the trailing edge towards the leading edge over the entire wing surface is an evidence on the backflow dominance over the incoming stream. Similarly, to what have been described in the previous angle of attack, the boundary layer is fully separated from the wing upper surface giving rise to a wing stalled state. In other hand, a micro rod is used for prospective wing performances enhancement. Some relevant findings are presented in the left-hand column of the Fig. 20. The instabilities in the average pressure distribution seem to be damped down to an extent where uniformity and homogenous appearance is recovered close to the leading edge. An analogy of this pressure contour with the baseline case, turns up this flow pattern resulting from the control device is already seen at earlier angle of attack, 11 • . This inferred that this passive control delays heavy stall by few degrees. There appear to be larger area of low pressure in the mid-span region of the wing as the main flow gain control over the backflow and sticks to the wing again through involvement of rod eddies that convey energy from free-stream to boundary layer. Whereas, close to the endplates, less suction is generated as the three-dimensional interfaces in the boundary layers of these two regions arise tiny secondary instabilities on both sides that block the flow from it. To complete what is said so far, friction lines on the upper surface of the wing at 15 • incidence with flow control are also considered. It is found that the flow got back its streamwise direction in comparison to the baseline case. This is because the boundary layer is completely separated from the wall and reattach to the wing particularly in the mid-span region constraining the stream to flow downward to the trailing edge. Two thirds of the wing suction side are occupied by a pair of counter rotating vortices resulting from recirculation of the flow from high to low pressure. In addition of secondary instabilities close to the leading edge on both sides are generated the flow from the endplates interferes with the wing boundary layer, which yield it to partially stalled state. Finally, all these results conclusively show that passive control as uses here showed capability to delay full stall by 3 • . Furthermore, an enhancement in aerodynamic loads is also ensured with control. 21% of drag reduction and 23% in lift growth are effectively obtained.

Conclusion

In this paper a passive flow control approach for a pre-stalled and post-stalled aerofoil is numerically investigated. Flow response to placing a microcylindrical shaped rod near NACA0012 wing-based profile at low chord-based Reynolds number of Re c = 4.45 × 10 5 is carefully examined. Particular focus is put on parietal pressure distribution on suction side surface along with flow field patterns with corresponding impart on the overall aerodynamic loads particularly high angles of attack ranging from 0 • to 24 • . The longitudinal gap between the wing and rod together with rod's diameter are key parameters carefully considered in this study. Detailed RANS CFD calculations are carried out for both controlled and baseline cases to show the influence of the joint device on the flow control effectiveness. Friction lines are further analysed to give information on the flow field behaviour on wall surface and correlative pressure distribution. The diameter of the microcylindrical rod considered for this current study is d/c = 0.013, placed away from the wing leading edge at 3 times the diameter in the chordwise direction. For the baseline case, flow fields and average pressure coefficient distribution analysis made it possible to highlight the underlying physical phenomenology involved in laminar and turbulent boundary layer separation. It is shown that laminar boundary layer is particularly sensitive to pressure gradients induced by the wing at incidence. This caused premature separation of the flow followed by reattachment after transition to turbulent state accompanied by formation of the so-called laminar separation bubble (LSB). At the pre-stall flow condition (α < 13 • ), the bubble forms near the trailing edge of the wing. As the incidence increases the LSB significantly decreases in length while migrating upstream to stabilize in the vicinity of the leading edge. This feature is encountered between α = 2 • and α = 8 • . When the angle of attack is greater than 8 • , the aerofoil exhibits presence of a second conventional bubble near the trailing edge as well as the LSB near the leading edge. At α = 12 • the laminar flow completely separates from the wing upper surface and forming a large recirculation zone; hence the wing is stalled. Any further increase in incidence results in the laminar separation bubble sudden burst leading to a heavy turbulent separation of the boundary layer subsequently to the important pressure gradient. Analysis of lift and drag forces with the LBS formed showed significant degradation in the overall aerodynamic loads pushing the wing to an abrupt stall. When the microcylindrical control device is introduced into the flow, a considerable influence on the behaviour of the laminar separation bubble is induced. Presence of the rod revealed ability to partially eliminate the bubble by causing an early reattachment of the premature laminar separated boundary layer over a significant extent. The total length of the LSB is found to reduce by more than 75% and almost vanish for an incidence of 10 • (Table 5). Furthermore, the vortex like generator deeply impacts the flow structure especially at stall and post-stall conditions. At AOAs ranged between 10 • to 15 • , the control device acts similarly to a pumping system energizing the boundary layer on the aerofoil counter-balancing the adverse pressure gradients. In fact, eddies generated in the cylinder wake interact with the boundary layer of the wing, ensuring convection of energy from freestream to near wall flow. The size of the downstream recirculation zone is significantly reduced for post-stall conditions. In other words, mean pressure contours and corresponding friction lines on the suction side are corroborating results for control effectiveness of the rod beyond stall conditions of the baseline wing. It was capable of maintaining attached flow for incidence of 15 • till 33% of the whole wing upper surface close to the midspan region, while the backflow took ad-vantage in non-controlled configuration, hence a complete separation from the leading edge was noticed. Finally, for pre-stall regime, the controlled flow doesn't provide significant improvement in either lift and drag, but does show effectiveness in reducing the length of the laminar separation bubble leading to further understanding of the intricate involved dynamics influencing the wing aerodynamic performances. Meanwhile, at post-stall regime presence of the rod close to the leading edge seems to effectively suppress or diminish and in most favourable case delay flow separation to higher angles of attack. This is advantageous to either lift enhancement or drag reduction, resulting respectively in a gain of 23% and a reduction of 73% (Fig. 24 and Fig. 25)as evaluated in this study (Fig. 22 and Fig. 23).

Nomenclature

Latin

Figure 1 :

 1 Figure 1: (a) 3D representation of the control set-up; (b) Schematic of the wing controlled by micro-cylinder placed at distance, L from the leading edge (flow moving left to right).

Figure 2 :

 2 Figure 2: Vortex formation length, l f as function of Reynolds number for a cylinder [26].

Figure 3 :

 3 Figure 3: Computational domain topology and boundary conditions.

Figure 4 :

 4 Figure 4: Variation of Lift and Drag coefficients for different spanwise lengths.

Figure 5 :

 5 Figure 5: (a) Schematic of the adopted mesh topology of the CFD domain (flow moving left to right); (b) Structured mesh around the rod close to the aerofoil; (c) Structured mesh close to the trailing edge of the aerofoil.Y

Figure 6 :

 6 Figure 6: Convergence history for scaled residuals.

Figure 7 :Figure 8 :

 78 Figure 7: Lift and Drag coefficients convergence for angle of attack 14 • .

Figures 9

 9 Figures 9 and 10 depict distribution of the parietal dimensionless distance, y + , on both wing surfaces and the rod circumference. The average value seems to be largely less than unity in order of 10 -1 . Thus, resolution of the boundary layer of the flow near both bodies is performed accurately enough throughout the grid used to conduct the different CFD simulations. Hence, spacing of the first cell from the wall is validated for this numerical study.

Figure 9 :

 9 Figure 9: Wall y + distribution on the wing at 10 • and 18 • of incidence.

Figure 10 :

 10 Figure 10: Wall y + distribution over the micro-cylindrical rod at angle of attack of 16 • .

  show variation of the aerodynamic coefficients, c L and c D , at low Reynolds number and various angles of attack for the baseline NACA0012 wing simulations along with a schematic of the flow patterns observed

Figure 11 : 6 D

 116 Figure 11: Comparison between the current CFD and experimental results with previous experimental data.CL_CFD_γ-Reθ,t, Present CL_Experiment, Present CL_Experiment, Jacobs et al.(1937) CD_CFD_γ-Reθ,t, Present CD_Experiment, Present CD_Experiment, Critzos et al.(1955)

Figure 12 :Figure 13 :

 1213 Figure 12: CFD lift coefficient for baseline case at varying AOAs.

Figure 14 :

 14 Figure 14: Process of LSB migration at low Reynolds number of 4.45 × 10 5 for NACA0012 wing; arrows highlights the LSB location (flow moving left to right).= 2 = 2°=

Figure 15 :

 15 Figure 15: LSB characterisation at flow conditions of Rec = 4.45 × 10 5 and α = 8 • (flow moving top to bottom).

Figure 16 : 4 .

 164 Figure 16: Effect of micro-cylinder rod on LSB dimensions at α = 10 • ; (a) baseline case; (b) controlled case (flow moving left to right).

Figure 17 :

 17 Figure 17: Contours of pressure coefficient (upper column) and friction lines distribution (lower column) on wing suction side surface at flow conditions of α = 11 • and Rec = 4.45 × 10 5 ; (a) baseline case; (b) controlled case (flow moving top to bottom).

Figure 18 :

 18 Figure 18: Contours of pressure coefficient (upper column) and friction lines distribution (lower column) on wing suction side surface at flow conditions of α = 12 • and Rec = 4.45 × 10 5 ; (a) baseline case; (b) controlled case (flow moving top to bottom).

Figure 19 :

 19 Figure 19: Contours of pressure coefficient (upper column) and friction lines distribution (lower column) on wing suction side surface at flow conditions of α = 13 • and Rec = 4.45 × 10 5 ; (a) baseline case; (b) controlled case (flow moving top to bottom).

Figure 20 :

 20 Figure 20: Contours of pressure coefficient (upper column) and friction lines distribution (lower column) on wing suction side surface at flow conditions of α = 15 • and Rec = 4.45 × 10 5 ; (a) baseline case; (b) controlled case (flow moving top to bottom).

Figure 21 :

 21 Figure 21: Contours of pressure coefficient (upper column) and friction lines distribution (lower column) on wing suction side surface at flow conditions of α = 10 • ; Rec = 4.45 × 10 5 and diffrent rods' size (flow moving top to bottom).

Figure 22 :

 22 Figure 22: Contours of pressure coefficient (upper column) and friction lines distribution (lower column) on wing suction side surface at flow conditions of α = 11 • ; Rec = 4.45 × 10 5 and diffrent rods' size (flow moving top to bottom).

Figure 23 :

 23 Figure 23: Contours of pressure coefficient (upper column) and friction lines distribution (lower column) on wing suction side surface at flow conditions of α = 12 • ; Rec = 4.45 × 10 5 and diffrent rods' size (flow moving top to bottom).

Figure 24 :

 24 Figure 24: Contours of pressure coefficient (upper column) and friction lines distribution (lower column) on wing suction side surface at flow conditions of α = 13 • ; Rec = 4.45 × 10 5 and diffrent rods' size (flow moving top to bottom).

Figure 25 :

 25 Figure 25: Contours of pressure coefficient (upper column) and friction lines distribution (lower column) on wing suction side surface at flow conditions of α = 15 • ; Rec = 4.45 × 10 5 and diffrent rods' size (flow moving top to bottom).

Table 1

 1 Rod location analysis

		Baseline Before	In for-	After	Faraway
	Location	forma-	mation	forma-	forma-
		tion	zone	tion	tion
		zone		zone	zone

Table 2

 2

		Grid independence analysis
	Mesh	Nodes	c L at	Error (%) with
		on	10 •	previous mesh
		Aerofoil
	MESH-1	200	0.8636	-
	MESH-2	450	0.9035	4.59
	MESH-3	640	0.9086	0.56

Table 3

 3

	Lift coefficient from CFD simulations and
	wind tunnel

Table 4

 4 Pearson's coefficient for different data sets

		Present		Present		Present
		CFD	vs.	CFD	vs.	Experi-
		Present		Literature	ment	vs.
		Experiment	Experiment	Literature
						Experiment
		C L	C D	C L	C D	C L	C D
	r	0.9743 0.9581 0.9935 0.9463 0.9797 0.9931
	3.2. Baseline Flow Analysis
	Figures 12 and 13

Table 5

 5 Rod effects on LSB size

	Angle of attack	5 •	8 •	10 •
	Size Reduction [%]	31.29	15.57	75


F Fabre

A Finco

A Purbawati

A Hadj-Azzem

N Rougemaille

Johann Coraux

Isabelle Philip

Vincent Jacques

Characterization of room-temperature in-plane magnetization in thin

Characterization of room-temperature in-plane magnetization in thin flakes of CrTe 2 with a single spin magnetometer F. Fabre, 1 A. Finco, 1 A. Purbawati, 2 A. Hadj-Azzem, 2 N. Rougemaille, 2 J. Coraux, 2 I. Philip, 1 and V. Jacques We demonstrate room-temperature ferromagnetism with in-plane magnetic anisotropy in thin flakes of the CrTe2 van der Waals ferromagnet. Using quantitative magnetic imaging with a single spin magnetometer based on a nitrogen-vacancy defect in diamond, we infer a room-temperature in-plane magnetization in the range of M ∼ 27 kA/m for flakes with thicknesses down to 20 nm. In addition, our measurements indicate that the orientation of the magnetization is not determined solely by shape anisotropy in micron-sized CrTe2 flakes, which suggest the existence of a nonnegligible magnetocrystalline anisotropy. These results make CrTe2 a unique system in the growing family of van der Waals ferromagnets, as it is the only material platform known to date which offers an intrinsic in-plane magnetization and a Curie temperature above 300 K in thin flakes.

I. INTRODUCTION

Ferromagnetic van der Waals (vdW) crystals offer numerous opportunities both for the study of exotic magnetic phase transitions in low-dimensional systems [1] and for the design of innovative, atomically-thin spintronic devices [2,3]. Since the discovery of a two-dimensional (2D) magnetic order in monolayers of CrI 3 [4] and Cr 2 Ge 2 Te 6 [5] crystals, the family of vdW ferromagnets has expanded very rapidly [6][7][8]. However, most of these compounds have a Curie temperature (T c ) well below 300 K, which appears as an important drawback for future technological applications. An intense research effort is therefore currently devoted to the identification of high-T c 2D magnets [3].

In this context, the vdW crystal Fe 3 GeTe 2 appears as a serious candidate because it can be grown in waferscale through molecular beam epitaxy and it exhibits a strong perpendicular magnetic anisotropy [9,10]. Although its intrinsic T c drops to 130 K in the monolayer limit [11], it might be raised above room-temperature either by ionic gating [12,13], interfacial engineering [14], or by micro-patterning, as demonstrated so far for rather thick films [15]. In addition, other FeGeTe alloys, such as Fe 4 GeTe 2 and Fe 5 GeTe 2 , exhibit high T c , still lower than room temperature but close to it [16,17]. Another promising strategy consists in incorporating magnetic dopants into 2D materials to form dilute magnetic semiconductors [18]. This approach was recently employed to induce room-temperature ferromagnetism in WSe 2 monolayers doped with vanadium [19,20]. Finally, an intrinsic ferromagnetic order was reported in epitaxial layers of VSe 2 [21], MnSe 2 [22] and VTe 2 [23] under ambient conditions, although the interpretation of these experiments still remains debated [24,25].

In this work, we follow an alternative research direction by studying the room-temperature magnetic properties of micron-sized flakes exfoliated from a CrTe 2 crystal with 1T structure. In its bulk form, this layered transition metal dichalcogenide is a ferromagnet with in-plane magnetization, i.e. pointing perpendicular to the c axis, and a T c around 320 K [26]. This combination of properties is unique in the growing family of vdW ferromagnets. Recent studies have reported that the magnetic order is preserved at room temperature in exfoliated CrTe 2 flakes with thicknesses in the range of a few tens of nanometers [27,28]. However, obtaining quantitative estimates of the magnetization in such micron-sized flakes remains a difficult task, which requires the use of non-invasive magnetic microscopy techniques combining high sensitivity with high spatial resolution. These performances are offered by magnetometers employing a single nitrogenvacancy (NV) defect in diamond as an atomic-size quantum sensor [29][30][31]. In recent years, this microscopy technique has found many applications in condensed matter physics [32], including the study of chiral spin textures in ultrathin magnetic materials [33][34][35], current flow imaging in graphene [36] and the analysis of the magnetic order in vdW magnets down to the monolayer limit [37][38][39].

Here we use scanning-NV magnetometry to infer quantitatively the in-plane magnetization in exfoliated CrTe 2 flakes under ambient conditions. Our measurements confirm that the ferromagnetic order is preserved in few tens of nanometers thick flakes, although with a low roomtemperature magnetization M ∼ 27 kA/m. This value is five times smaller that the one measured in a bulk CrTe 2 crystal. Such a reduction of the magnetization is attributed to a decreased Curie temperature in exfoliated flakes. Moreover, our results show that shape anisotropy alone does not fix the in-plane orientation of the magnetization in micron-sized CrTe 2 flakes, pointing out the existence of a substantial magnetocrystalline anisotropy.

II. MATERIALS AND METHODS

A bulk 1T -CrTe 2 crystal was synthesized following the procedure described in Ref. [26]. The in-plane magnetization of this layered ferromagnet was first characterized as a function of temperature through vibrating sample mag-netometry under a magnetic field of 500 mT. The results shown in Fig. 1(a) indicate a Curie temperature around 320 K and a magnetization reaching M ∼ 120 kA/m under ambient conditions. CrTe 2 flakes with thicknesses ranging from a few tens to a hundred of nanometers were then obtained by mechanical exfoliation and transferred on a SiO 2 /Si substrate. We note that the probability to obtain thin CrTe 2 flakes through mechanical exfoliation is still very low compared to other layered transition metal dichalcogenides, such as MoS 2 or WSe 2 [28]. The thinnest flake studied in this work has a thickness of ∼ 20 nm. Like all van der Waals ferromagnets known to date, CrTe 2 flakes are unstable under oxygen atmosphere. However, a recent study combining X-ray and Raman spectroscopy has shown that oxidation of CrTe 2 flakes occurs typically within a day scale under ambient conditions and is limited to the very first outer layers [28]. In this work, CrTe 2 flakes were not encapsulated and all the measurements were done within a day after exfoliation to mitigate oxidation.

Magnetic imaging was performed with a scanning-NV magnetometer operating under ambient conditions [31]. As sketched in Fig. 1(b), a single NV defect integrated into the tip of an atomic force microscope (AFM) was scanned above CrTe 2 flakes to probe their stray magnetic fields. At each point of the scan, a confocal optical microscope placed above the tip was used to monitor the magnetic-field-dependent photoluminescence (PL) properties of the NV defect under green laser illumination. In this work, we employed a commercial diamond tip (Qnami, Quantilever MX) with a characteristic NV-tosample distance d NV = 80 ± 10 nm, as measured through an independent calibration procedure [40]. Two different magnetic imaging modes were used. In the limit of weak stray fields (< 5 mT), quantitative magnetic field mapping was obtained by recording the Zeeman shift of the NV defect electron spin sublevels through optical detection of the electron spin resonance (ESR). This method relies on microwave driving of the NV spin transition combined with the detection of the spin-dependent PL intensity of the NV defect [41]. For stronger magnetic fields (> 5 mT), the scanning-NV magnetometer was rather used in all-optical, PL quenching mode [42,43]. In this case, localized regions of the sample producing large stray fields are simply revealed by an overall reduction of the PL signal induced by a mixing of the NV defect spin sublevels [44]. We note that the diameter of the scanning diamond tip is around 200 nm in order to act as an efficient waveguide for the PL emission of the NV defect [45,46]. As a result, such tips cannot provide precise topographic information of the sample. For the thinnest flakes, the topography was thus imaged by using conventional, sharp AFM tips.

For a ferromagnetic material, stray magnetic fields can be produced by magnetization patterns presenting a nonzero divergence ∇ • M = 0. Considering a CrTe 2 flake, this can occur (i) at the edges of the flake, (ii) at the location of non-collinear spin textures such as domain walls, and (iii) at the position of thickness steps on the flake. In this work, the magnetization is inferred from the measurement of the stray field produced at the edges of the flake [Fig. 1(b)]. Therefore, the analysis is made simpler for uniformly magnetized flakes with a homogeneous thickness. While uniform spin textures can be readily obtained by applying a bias magnetic field, obtaining thin CrTe 2 flakes with a uniform thickness through mechanical exfoliation is currently achieved with a low probability. Consequently, variations in thickness must be carefully taken into account when analyzing the experimental results.

III. RESULTS AND DISCUSSION

In a first experiment, we imaged a CrTe 2 flake with a large thickness t ∼ 150 nm [Fig. 2(a)]. Considering the saturation magnetization of the bulk CrTe 2 crystal, magnetic simulations predict stray field amplitudes larger than 10 mT at a distance d NV ∼ 80 nm above the edges of the 150-nm-thick flake. The scanning-NV magnetometer was thus operated in the PL quenching mode for such a thick flake. Fig. 2(b) displays the resulting PL map. Several observations can be made. First, a quenching of the PL signal is observed when the NV defect is placed above the flake. This quenching effect, which is constant all over the flake, does not have a magnetic origin. It is linked to the metallic character of CrTe 2 [47,48]. Second, a stronger PL quenching is obtained at the flake edges, as expected for a single ferromagnetic domain. We note that the additional quenching spot observed at the topright edge of the flake results from the stray field produced by a large variation of the thickness [Fig. 2(a)]. Given the different sources of PL quenching involved in this experiment, quantitative estimates of the magnetization cannot be obtained. However, the analysis of the PL quenching distribution can be used to extract qualitative information about the orientation of the magnetization. To this end, simulations of the PL map were carried out considering only the effect of magnetic fields. The geometry of the flake used for the simulation was inferred from the AFM image simultaneously recorded with the NV microscope [Fig. 2(a)]. Assuming a uniform in-plane magnetization M with an azimuthal angle φ M in the (x, y) plane, the stray field was first calculated at a distance d NV above the flake. The resulting PL quenching map was then simulated by using the model of magnetic-field-dependent photodynamics of the NV defect described in Ref. [44]. Simulated PL maps obtained for two different orientations of the magnetization are shown in Figs. 2(c,d). In micron-sized flakes, shape anisotropy should favor a magnetization pointing along the long axis of the flake. Interestingly, the PL map simulated for such a magnetization direction disagrees with the experimental data [Fig. 2(c)], which instead suggest a magnetization pointing perpendicular to the long axis of the flake [Fig. 2(d)]. This result is a first indication of a non-negligible, in-plane magnetocrystalline anisotropy. Thinner CrTe 2 flakes, i.e. producing less stray field, were then studied through quantitative magnetic field imaging. To this end, a microwave excitation was applied through a copper microwire directly spanned on the sample surface and the Zeeman shift of the NV defect's ESR frequency was recorded at each point of the scan. In the weak field regime, the ESR frequency is shifted linearly with the projection B NV of the stray magnetic field along the NV defect quantization axis [31]. This axis was precisely measured by applying a calibrated magnetic field [49], leading to the spherical angles θ NV = 58 • and φ NV = 103 • in the laboratory frame (x, y, z), as illustrated in Fig. 1(b). For these measurements, a bias field of 3 mT was applied along the NV defect axis in order to infer the sign of the stray magnetic field [31]. Furthermore, this bias field is strong enough to erase spin textures such as domain walls, given the weak coercitive field of CrTe 2 flakes [28].

A map of B NV recorded above a 35-nm thick CrTe 2 flake is shown in Fig. 3(b). A stray magnetic field around ±1.5 mT is detected at two opposite edges of the flake. In principle, the underlying sample magnetization can be retrieved from such a quantitative magnetic field map by using reverse propagation methods with well-adjusted filters in Fourier-space [32]. Under several assumptions, this method can be quite robust for magnetic materials with out-of-plane magnetization [37][38][39]. For in plane magnets, however, the reconstruction procedure amplifies noise and is thus much less efficient [50]. As a result, the recorded magnetic field distribution was rather directly compared to magnetic calculations in order to extract quantitative information on the sample magnetization.

To obtain precise information on the geometry of the flake, the topography image was here measured with a sharp AFM tip [Fig. 3(a)]. The observed variations in thickness were included in the magnetic calculation (see Appendix A). Considering a uniformly magnetized flake with an azimuthal angle φ M , the stray field was calculated at a distance d NV above the flake and then projected along the NV quantization axis in order to simulate a map of B NV . By comparing experimental data with magnetic maps simulated for different values of the angle φ M , the orientation of the in-plane magnetization was first identified, leading to φ M = -145 ± 5 • , i.e. pointing along the short axis of the flake [Fig. 3(c)]. Once again, this result suggests the presence of a nonnegligible magnetocrystalline anisotropy, in agreement with recent works [28]. The norm M of the magnetization vector was then estimated by fitting experimental line profiles with the result of the calculation, leading to M = 27±4 kA/m [Fig. 3(d)]. We note that the stray field amplitude depends on several parameters including d NV , φ M , the flake thickness t and the NV defect quantization axis (θ NV , φ NV ). Any imprecisions on these parameters directly translate into uncertainties on the evaluation of the magnetization M . A detailed analysis of uncertainties is provided in Appendix B.

The room-temperature magnetization measured in the exfoliated CrTe 2 flake is almost five times smaller than the one obtained in the bulk crystal. This observation could be explained by a degradation of the sample surface through oxidation processes, leading to a thinner effective magnetic thickness. However, recent experiments relying on X-ray magnetic circular dichroism coupled to photoemission electron microscopy (XMCD-PEEM) have shown that oxidation is limited to the first outer layers of the CrTe 2 flake [28]. Considering a 35-nm thick CrTe 2 flake, surface oxidation can thus be safely neglected, and cannot explain the observed reduction of the magnetization. This effect is rather attributed to a decrease of the Curie temperature in exfoliated flakes, a phenomenon that has been observed in other vdW magnets such as Fe 3 GeTe 2 below 5-10 nm thicknesses [11,51], and in more traditional ferromagnetic thin films below few nanometers thicknesses [52,53]. The thickness marking the crossover from a thin film (2D-like) to a bulk magnetism (3D) is expected to be strongly material-dependent and cannot be predicted a priori in the case of CrTe 2 [8]. However, for a 35-nm thick flake, bulk-like magnetism is a reasonable assumption and the magnetization's amplitude should not be altered by the film thickness. In this thickness regime, the main parameter altering the magnetization's amplitude is temperature, and how close or far it is from T c . The estimation of the reduction in Curie temperature was tentatively inferred by translating the Curie law measured for a bulk CrTe 2 crystal [Fig. 1(a)] in order to reach M ∼ 27 kA/m at room temperature. This is achieved for a reduction of T c by ∼ 20 K. This value is in line with a recent study, which estimates a Curie temperature around 300 K for CrTe 2 flakes with a few tens of nanometers thickness [27].

To support this finding, a similar analysis was performed for a 20-nm thick flake [Fig. 4(a)]. Here, a stray magnetic field is mainly detected along the bottom edge of the flake. On the opposite edge, the stray field distribution is quite inhomogeneous [Fig. 4(b)]. This observation is attributed to damages of the flake, which can be observed in the AFM image. It is difficult to describe the corresponding complex thickness variations, and hence to take them into account in the simulations. These height variations seem to correspond to several CrTe 2 islands, whose very small sizes could make them more prone to oxidation than larger flakes, and which may exhibit random magnetization orientations. We hence disregarded the thickness variations in our structural model and the magnetic calculation was performed for a flake with uniform thickness, which is a good approximation for the bottom and left edges of the flake. A simulation of the stray field distribution for a magnetization with an azimuthal angle φ M = -100±10 • reproduces fairly well the experimental results [Fig. 4(c)] and a quantitative analysis of line profiles across the bottom edge of the flake leads to M = 26 ± 4.0 kA/m [Fig. 4(d)], a similar value to that obtained for the 35-nm thick flake. This observa-tion indicates that the magnetization is not significantly modified for thicknesses lying in the few tens of nanometers range. A more in-depth analysis of the dependence of magnetization on thickness could not be carried out at this stage, given the very low probability to obtain thin CrTe 2 flakes of homogeneous thickness through mechanical exfoliation.

IV. CONCLUSION

To conclude, we have used quantitative magnetic imaging with a scanning-NV magnetometer to demonstrate that exfoliated CrTe 2 flakes with thicknesses down to 20 nm exhibit an in-plane ferromagnetic order at room temperature with a typical magnetization in the range of M ∼ 27 kA/m. These results make CrTe 2 a unique system in the growing family of vdW ferromagnets, because it is the only material platform known to date which offers an intrinsic in-plane magnetization and a T c above room temperature in thin flakes. These properties might offer several opportunities for studying magnetic phase transition in 2D-XY systems [54] and to design spintronic devices based on vdW magnets. The next challenge will be to assess if the ferromagnetic order is preserved at room temperature in the few layers limit.

APPENDIX A: MAGNETIC FIELD SIMULATION

As indicated in the main text, thickness variations within the flake can result in a magnetization pattern with a non-zero divergence which produces a stray magnetic field. When possible, these variations were taken into account in the magnetic calculation. This is illustrated in Fig. 5(a,b), where the geometry of the flake used for the calculation includes a thickness step observed in the topography image. In Fig. 5(c-e), we show the magnetic field distributions produced at a distance d NV = 80 nm for three different magnetization orientations of the flake. First considering an out-of-plane magnetization [Fig. 5(c)], the simulated magnetic field distribution does not agree with the experimental data shown in Fig. 3. This confirms that the magnetization in lying in the plane of the CrTe 2 flake, i.e. perpendicular

This analysis was performed for each parameter p i . The cumulative uncertainty is finally given by =

2 fit + i 2 pi , (2)
where all errors are assumed to be independent.

A summary of the uncertainties is given in Table I. We obtain M = 27 ± 4 kA/m for the flake shown in Fig. 3 and M = 26 ± 4 kA/m for the flake shown in Fig. 4.

2 FIG. 1 .

 21 FIG. 1. (a)Temperature dependence of the magnetization M in a bulk crystal of CrTe2 with 1T polytype. The measurement is performed by vibrating sample magnetometry under a magnetic field of 500 mT. At room temperature, the magnetization is around 120 kA/m (black dashed lines). The inset shows the layered crystal structure of 1T -CrTe2 . (b) A single NV defect (red arrow) localized at the apex of a diamond tip is scanned above exfoliated CrTe2 flakes. A microscope objective placed above the tip is used to collect the photoluminescence (PL) of the NV defect under green laser excitation. In this work, the NV-to-sample distance is dNV ∼ 80 nm and the quantization axis of the NV defect is characterized by the spherical angles (θNV = 58 • , φNV = 103 • ) in the laboratory frame of reference (x, y, z). The black arrows indicate the magnetic field lines generated at the edges of a CrTe2 flake with uniform in-plane magnetization M.

FIG. 2 .

 2 FIG. 2. (a) AFM image of a 150 nm-thick CrTe2 flake recorded with the scanning-NV magnetometer. (b) Corresponding PL map showing a quenching contour at the edges of the flake. The shape of the flake is indicated by the black dashed lines. The recorded PL signal is normalized with the one measured far from the flake. The experiment is performed at zero external magnetic field. (c,d) Simulated PL maps for two different orientations φM of the in-plane magnetization: (c) φM = 20 • and (d) φM = -80 • . Here, the norm of the magnetization is fixed arbitrarily to its bulk value M = 120 kA/m. The PL quenching induced by energy transfer between the NV sensor and the metallic sample is not included in the simulations.

FIG. 3 .

 3 FIG. 3. (a) AFM image of a 35 nm-thick CrTe2 flake. (b) Corresponding distribution of the magnetic field component BNV. (c) Simulated magnetic field distribution BNV for a uniform in-plane magnetization M (black arrow) with an azimuthal angle φM = -145 • in the (x, y) plane and a norm M = 27 kA/m. The calculation is done at a distance dNV = 80 nm above the flake, whose shape is extracted from the topography image [black dashed lined in (a)]. (d) Fitting of an experimental line profile (black markers) with the magnetic calculation (blue solid line). The position of the linecuts are indicated by the white dashed lines in (b) and (c). A magnetization M = 27 ± 4 kA/m is obtained, the uncertainty being illustrated by the blue shaded area. The red and green solid lines indicate the result of the calculation for M = 40 kA/m and M = 15 kA/m, respectively.

FIG. 4 .

 4 FIG. 4. (a) AFM image of a 20-nm thick CrTe2 flake. (b) Corresponding map of the magnetic field component BNV. (c) Simulated map of BNV for a uniform in-plane magnetization M (black arrow) with an azimuthal angle φM = -100 • in the (x, y) plane and a norm M = 26 kA/m. The black dashed line indicate the shape of the flake used for the simulation. (d) Fitting of an experimental line profile (black markers) with the magnetic simulation (blue solid line). The linecuts are indicated by the white dashed lines in (b) and (c). A magnetization M = 26 ± 4.0 kA/m is obtained, the uncertainty being illustrated by the blue shaded area. The inset shows a linecut across the flake (white dashed lines in (a)).

FIG. 5 .

 5 FIG. 5. (a) AFM image of CrTe2 flake shown in Fig. 3. (b) Geometry of the flake used for the magnetic calculation, which includes the thickness step observed in the AFM imafge. (c-e) Calculated maps of BNV for a uniform magnetization M pointing out-of plane (c), and in-plane with an azimuthal angle φM = -30 • (d) or φM = -145 • (e). The norm of the magnetization is fixed to M = 27 kA/m.

1 1

 1 Laboratoire Charles Coulomb, Université de Montpellier and CNRS, 34095 Montpellier, France 2 Univ. Grenoble Alpes, CNRS, Grenoble INP, Institut NEEL, 38000 Grenoble, France (Dated: February 12, 2021)

TABLE I .

 I Analysis of uncertainties in the measurement of the magnetization (a) CrTe2 flake shown in Figure3parameter pi nominal value pi uncertainty σp i p i (%)

	dNV	80 nm	10 nm	8
	θNV	58 •	2 •	7
	φNV	103 •	2 •	2
	t	35 nm	2 nm	6
	φM	-145 •	5 •	5
	=	2 fit + i	2 p i	14
	(b) CrTe2 flake shown in Figure 4
	parameter pi nominal value pi uncertainty σp i	p i (%)
	dNV	80 nm	10 nm	8
	θNV	58 •	2 •	7
	φNV	103 •	2 •	2
	t	24 nm	2 nm	8
	φM	-100 •	10 •	12
	=	2 fit + i	2 p i

ACKNOWLEDGEMENTS

The authors warmly thank J. Vogel and M. Núñez-Regueiro for fruitful discussions. This research has received funding from the European Union H2020 Program under Grant Agreement No. 820394 (ASTERIQs), the DARPA TEE program, and the Flag-ERA JTC 2017 project MORE-MXenes. A.F. acknowledges financial support from the EU Horizon 2020 Research and Innovation program under the Marie Sklodowska-Curie Grant Agreement No. 846597 (DIMAF).

APPENDIX B: ANALYSIS OF UNCERTAINTIES

The norm of the magnetization M is obtained by fitting line profiles of the measured stray field distribution with the result of the magnetic calculation [see Fig. 3(d)]. In this section, we analyze the uncertainty of this measurement by using the methodology described in Ref. [55]. The uncertainties result (i) from the fitting procedure itself and (ii) from those on the parameters p i = {d NV , θ NV , φ NV , t, φ M } which are all involved in the magnetic calculation. In the following, the parameters p i are expressed as p i = pi + σ pi , where pi denotes the nominal value of parameter p i and σ pi its standard error. These parameters are independently evaluated as follows:

• The probe-to-sample distance d NV is inferred through a calibration measurement, following the procedure described in Ref. [40], leading to d NV = 80 ± 10 nm.

• The NV defect quantization axis is measured by recording the ESR frequency as a function of the • The thickness t of the CrTe 2 flake is extracted from line profiles of the AFM image with an uncertainty of ±2 nm.

• The azimuthal angle of the in-plane magnetization φ M is obtained through the comparison between experimental data and simulated magnetic field maps.

We first evaluate the uncertainty of the fitting procedure. To this end, the line profile is fitted with the result of the magnetic calculation while fixing all the parameters p i to their nominal values pi , leading to M = 26.9 ± 0.7 kA/m. The relative uncertainty linked to the fitting procedure is therefore given by fit = 3%. We note that the intrinsic accuracy of the magnetic field measurement is in the range of δB NV ∼ 5 µT. The resulting uncertainty can be safely neglected.

In order to estimate the relative uncertainty pi introduced by each parameter p i , the fit was performed with one parameter p i fixed at p i = pi ± σ pi , all the other parameters remaining fixed at their nominal values. The corresponding fit outcomes are denoted M ( pi + σ pi ) and M ( pi -σ pi ). The relative uncertainty pi introduced by