Shear Layer Research Articles

Effect of Cavity Geometry and Location on Base Pressure in a Suddenly Expanded Flow at Mach 2.0 for Area Ratio 3.24

The primary goal of this study is to use CFD analysis to investigate the impact that a cavity has on the pressure at the base of a structure. In this analysis, we took into account the NPR, the cavity aspect ratio, and the cavity position. In this case, the area ratio is 3.24 and the Mach number is 2.0. Simulations were run with L/D ratios between 1 and 6, and NPRs of 3, 5, 7.8, 9, and 12. The 2-dimensional model was developed using ANSYS Fluent's Design Modeler. The nozzle is operating at Mach 2.0. Base pressure and wall pressure in the duct were the primary research foci. The C-D nozzle was created for this research. ANSYS Fluent was used to verify the CFD findings. When the nozzles are under-expanded and the cavity is at 0.5D, passive control as a cavity is shown to be effective. It appears that 1D is the bare minimum for duct length. Because the shear layer gets reattached to the duct wall at 1D and the boundary layer grows after reattachment, passive control is not observed in the flow process regardless of whether the cavity is located at 1D, 1.5D, 2D, or 3D. An oscillating base pressure is seen at shorter duct lengths. This phenomenon does not occur at longer duct lengths. Whether or not there are cavities in the duct, the flow field is the same.

A Study of Cavitation Erosion in Artificial Submerged Water Jets

The artificially submerged cavitation water jet is effectively utilized by ejecting a high-pressure water stream into a low-pressure water stream through concentric nozzles and utilizing the cavitation phenomenon generated by the shear layer formed between the two streams. In this study, we investigated the cavitation characteristics of artificially submerged cavitation water jets by combining numerical simulations and erosion experiments. The results indicate that an appropriate standoff distance can generate more cavitation clouds on the workpiece surface, and the erosion characteristics of the artificially submerged cavitation water jet are most pronounced at a dimensionless standoff distance of SD = 30. The shear effect formed between the two jets plays a crucial role in generating initial cavitation bubbles within the flow field of the artificially submerged cavitation water jet. Moreover, increasing the convergent angle between the two jets can significantly enhance the cavitation effect between them and lead to a more substantial cavitation effect. Simultaneously, increasing the pressure of the high-pressure inner nozzle also contributes to enhancing the cavitation effect of the artificially submerged cavitation water jet.

Turbulence generation in transonic shock diffraction

In this Letter, we have quantified the zones of turbulence in transonic shock diffraction across a 90° step corner using Lamb vector analysis. We have analyzed the nature of the three-dimensional turbulence structures using invariants of the velocity gradient tensor. From our results, we conclude that the vortex-induced shock, the lambda shock on the shear layer, and the Kelvin–Helmholtz instability forming on the shear layer are the main contributors in turbulence generation. In our study, we have shown that two-dimensional simulations can resolve the gasdynamic behavior accurately, though three-dimensional simulations are required to understand turbulent structures.

Development and Verification of an Improved Wake-Added Turbulence Model in FAST.Farm

We introduce a generalized wake-added turbulence (WAT) model in the multiphysics, multiturbine simulation tool FAST.Farm. The WAT model introduces additional small-scale turbulence that represents the breakdown of vortical structures and shear layers in the wake. The article describes the development, implementation, calibration, and verification of the model. The novelties of the model include support for wake asymmetry, buildup of WAT across the wind farm, and secondary effects of wake-induced turbulence (e.g., wake meandering) driven by smaller-scale turbulence structures that arise from wake breakdown. Large-eddy simulations were run to support the calibration of the WAT parameters and verification of the model. Previous studies hypothesized that the lack of WAT modeling was the source of underprediction of fatigue loads, in particular for cases at low turbulence intensities and/or stable atmospheric boundary layers. This study confirms that the newly implemented WAT model enhances the loads predictions in these cases.

Two-dimensional direct numerical simulation study of multicomponent mixing with phase transition in a transcritical shear layer

In this paper, we investigate two-dimensional direct numerical simulations of a transcritical shear layer. Three configurations are chosen, which are distinguished by the level of presence of two-phase phenomena. The thermodynamic model is based on a cubic equation of state. It was extended for multicomponent mixtures, and it is able to account for vapor–liquid equilibrium. The thermodynamic modeling with phase-transition is validated using experimental data from the literature. Special focus is put on the effect of the density gradient and the density changes caused by phase-transition on the development of the turbulent shear layer and the associated mixing. In addition to this, the vorticity distribution and the components of its transport equation are analyzed and compared for the different configurations.

Vortex-induced noise suppression of a cylinder with blowing through porous media

To mitigate vortex shedding for flow and noise control of a circular cylinder, an experimental approach combining air blowing and porous coating was implemented simultaneously as a hybrid method. Localized air blowing was symmetrically applied through structured porous media at four angles, corresponding to different regions of the flow field: boundary layers, shear layers on the cylinder, separated shear layers, and the cylinder's base. The study involved synchronizing near-field pressure fluctuation and far-field noise measurements with flow field measurements obtained via particle image velocimetry. Near-field pressure measurements were taken around the cylinder's circumference using a remote sensing method. This comprehensive investigation revealed that vortex shedding primarily induced pressure fluctuations at the cylinder's shoulders, resulting in the propagation of acoustic waves to the far field. The hybrid method, alongside the separate application of porous coating and local blowing, showcased substantial efficacy in mitigating near-field pressure, consequently leading to a reduction in far-field noise. These techniques achieved this by strategically shifting the vortex formation region further downstream and expanding the wake region compared to the baseline. Notably, the hybrid method, particularly when local blowing was applied at the base of the porous coated cylinder, exhibited a significantly enhanced impact in this regard, resembling the behavior observed with the individual application of porous coating.

Influence of freestream turbulence on boundary layer transition over a controlled-diffusion compressor blade

The influence of inlet freestream turbulence (FST) on the boundary layer transition over the suction surface of a controlled-diffusion compressor blade is demonstrated here by employing a well-resolved large-eddy simulation. Inherent to low Reynolds number conditions, a laminar separation bubble (LSB) forms on the suction surface, attributing to substantial flow diffusion. Inlet FST levels ranging from 1.5% to 7.6% are systematically varied, while maintaining a constant Reynolds number based on axial chord and inlet velocity at 2.1 × 105. Transition of the shear layer is initiated via Kelvin–Helmholtz instability with the amplification of selective frequencies until an inlet FST of 2.3%. Secondary instability emerges in the second half of the LSB, attributed to the amplification of perturbations in the braid region, ultimately leading to breakdown near the reattachment. At a moderate FST level of 4.2%, longitudinal streaks in the first half of the blade elongate downstream, causing the LSB to disappear, while the flow becomes inflectional at the mid-chord. Thus, the boundary layer transition in the second half of the blade is attributed to the high receptivity of the inflectional layer and breakdown of streaks, leading to an exponential growth of disturbances. Finally, at an inlet FST of 7.6%, the boundary layer appears pre-transitional in the first half of the blade, exhibiting significant turbulence levels. In the latter half, excitation occurs primarily through the breakdown of streaks, reflecting an algebraic growth of disturbances. Flow features and oscillations in the Nusselt number in this case suggest the outer mode of streak-induced instability.

A high-order upwind compact difference scheme for solving the streamfunction-velocity formulation of the unsteady incompressible Navier–Stokes equations

In this paper, we propose an upwind compact difference method with fourth-order spatial accuracy and second-order temporal accuracy for solving the streamfunction-velocity formulation of the two-dimensional unsteady incompressible Navier–Stokes equations. The streamfunction and its first-order partial derivatives (velocities) are treated as unknown variables. Three types of compact difference schemes are employed to discretize the first-order partial derivatives of the streamfunction. Specifically, these schemes include the fourth-order symmetric scheme, the fifth-order upwind scheme, and the sixth-order symmetric scheme derived by combining the two parts of the fifth-order upwind scheme. As a result, the fourth-order spatial discretization schemes are established for the Laplacian term, the biharmonic term, and the nonlinear convective term, along with the Crank–Nicolson scheme for the temporal discretization. The unconditional stability characteristic of the scheme for the linear model is proved by discrete von Neumann analysis. Moreover, six numerical experiments involving three test problems with the analytic solutions, and three flow problems including doubly periodic double shear layer, lid-driven cavity flow, and dipole-wall interaction are carried out to demonstrate the accuracy, robustness, and efficiency of the present method. The results indicate that the present method not only has good numerical performance but also exhibits quite efficiency.

Flow-induced vibrations of a twistable two-tandem cylinder system across varied gap ratios

This paper presents a simulation study on the flow-induced vibrations (FIVs) of a two-tandem cylinder system with three degrees of freedom, encompassing in-flow, cross-flow, and torsional motions at a Reynolds number of 150. We specifically focus on investigating the impact of gap ratios (G*) between the cylinders, considering values of 0.5, 2, and 3. The ratio of torsional natural frequency to translational natural frequency was maintained at 12.5, while the mass ratio between the cylinder system and the fluid was set to 2. Our findings reveal a notable sensitivity of the FIV responses to the gap ratio between the cylinders. It was observed that the cylinder system's vibrations can undergo locking phenomena at both the translational and torsional natural frequencies, with maximal cross-flow amplitudes occurring within the torsional locked region. Particularly, at G* = 0.5, the maximal cross-flow amplitude reached approximately 21.4 times that of the cylinder system without torsional vibration, posing potential risks of fatigue failure in offshore equipment. Conversely, at G* = 2, an unusual ultra-low frequency was detected in the in-flow vibration. This phenomenon significantly magnified the FIV of the cylinder system, despite the absence of frequency locking. Furthermore, at G* = 3, we observed an expanded range of vibrations with ultra-low frequency compared to G* = 2. Our analysis identified the differential development speed of the shear layers around the upstream and downstream cylinders as a key factor contributing to varying vortex shedding frequencies. This discrepancy led to periodic changes in the vortex shedding mode of the cylinder system over time, thereby introducing low-frequency components into the flow-induced vibration. These insights deepen our understanding of the complex dynamics governing the FIV of tandem cylinder systems, which holds implications for the design and maintenance of offshore structures.

Entropy damping and Bulk Viscosity based artificial compressibility methods on dynamically distorting grids

Artificial Compressibility Methods (ACM) rely on an artificial equation that links the pressure and velocity fields to model incompressible flows. These hyperbolic/parabolic equations can rapidly converge to a ‘nearly’ divergence-free flow field in contrast to the methods based on the elliptic pressure Poisson equation. We compare the computational efficacy of two ACMs, namely, the Bulk Viscosity ACM (BVACM) and Entropically Damped Artificial Compressibility (EDAC) recently proposed in the literature. The methods implemented in the in-house high-order finite difference solver, COMPSQUARE, are validated for the test cases of a 2D doubly periodic shear layer (DPSL), a 3D Taylor Green Vortex (TGV), and 2D/3D NACA0012 airfoil pitching about the quarter chord. The efficacy of these methods was also tested on static and dynamic grids using conservative metrics. Although both ACMs yield competitive results, the divergence of the velocity field is found to be more prominent in the highly unsteady regions. BVACM resulted in (a) a superior divergence-free velocity field and (b) 20−38% higher maximum stable time than the EDAC, thereby increasing the computational speed. A higher value of the bulk viscosity coefficient, A, although ensures a stringent divergence-free velocity field, is shown to have minimal effect on the flow statistics and reduce the maximum stable time step. The parabolic–hyperbolic nature of the governing equations and the lack of dual time-stepping in BVACM and EDAC ensures that both these methods are highly scalable on massively parallel architectures. Since the energy equation is no longer required to compute the velocity field, both EDAC and BVACM approaches are found to be 8−10% faster than the weakly compressible Navier–Stokes simulations under the low-Mach number limit.

Exploring the Substructure of the Near-surface Shear Layer of the Sun

The gradient of rotation in the near-surface shear layer (NSSL) of the Sun provides valuable insights into the dynamics associated with the solar activity cycle and the dynamo. Results obtained with global oscillation mode splittings lack resolution near the surface, prompting the use of the local helioseismic ring-diagram method. While the Helioseismic and Magnetic Imager ring-analysis pipeline has been used previously for analyzing this layer, default pipeline parameters limit the accuracy of the near-surface gradients. To address these challenges, we fitted the flow parameters to power spectra averaged over one-year periods at each location, followed by additional averaging over 12 yr. We find that the NSSL can be divided into three fairly distinct regions: a deeper, larger region with a small shear, steepening toward the surface; a narrow middle layer with a strong shear, with a gradient approximately 3 times larger; and a layer very close to the surface, where the logarithmic gradient is close to zero but becomes steeper again toward the surface. The middle layer appears to be centered at 3 Mm, but the poor resolution in these layers implies that it is potentially located closer to the surface, around 1.5 Mm deep. While our analysis primarily focused on regions along the central meridian, we also investigated systematic errors at longitudes off the center. The east–west antisymmetric component of the gradient reveals a layer of substantial differences between the east and west longitude of around 1.7 Mm, and the amplitude of the differences increases with the longitude.

Exact three-dimensional elasticity analysis for buckling of composite laminated plates resting on viscoelastic foundation

This paper aims to present novel exact solutions for the buckling of a laminated plate resting on the viscoelastic foundation with both normal and shear viscoelastic layers. The governing equations of plate buckling are derived using three-dimensional elasticity theory and state-space formulation. The normal and shear layers of the viscoelastic foundations are modeled using the generalized Maxwell model to represent both the elastic and viscose properties of the foundation. To couple the viscoelastic foundation equation with the buckling equation, Boltzmann’s superposition principle along with the Laplace transform is utilized. Then, the effects of geometry, relaxation modulus of normal and shear layers, viscosity, and time are investigated on the buckling load. The results reveal that the higher viscosity coefficient leads to a slower rate of change in the buckling loads. In addition, the viscoelastic properties have a significant impact on the buckling behavior of the plate. In this regard, the results show that instead of the expected second mode at a constant aspect ratio, the plate experiences the first mode as time passes. The computed results also show that there is a critical threshold. When the foundation stiffness exceeds this threshold, the conventional method of reducing the aspect ratio to prevent buckling not only proves ineffective in reducing the probability of buckling but also, in fact, leads to an increase in buckling occurrences. In addition to the analytical investigation, a finite element (FE) analysis is carried out to study the buckling response of the composite plate. The finite element results also show a reasonably good agreement with those of the analytical method.

Experimental Parametric Study on Flow Separation Control Mechanisms around NACA0015 Airfoil Using a Plasma Actuator with Burst Actuation over Reynolds Numbers of 105–106

Dielectric barrier discharge plasma actuators (DBD-PAs) have the potential to improve the performance of fluid machineries such as aircrafts and wind turbines by preventing flow separation. In this study, to identify the multiple flow control mechanisms in high Reynolds number flow, parametric experiments for an actuation parameter F+ with a wide range of Re values (105–106) for NACA0015 airfoil was conducted. We conducted wind tunnel tests by applying a DBD-PA to the flow field around a wing model at the leading edge. Lift characteristics, turbulent kinetic energy in the flow field, shear layer height, and the separation point of the boundary layer were evaluated based on pressure distributions on the wing surface and velocity of the flow field, with the effect of DBD-PA on the post-stall flow around the wing and the mechanism behind the increase in the lift coefficient CL analyzed based on these evaluation results. The following phenomena contributed to the increase in CL: (1) increase in turbulent kinetic energy; (2) increase in circulation; and (3) acceleration of the flow near the leading edge. Thus, this study effectively investigated the dependence of the increase in lift on F+ and the lift-increasing mechanism for a wide range of Re values.

Scalar transport after a high-resolution solitary fractal tree based on large-eddy simulation: Implication to urban green infrastructure

Scalar transport is one of the most concerned issues in the urban environment (e.g., pollutants, pollens, and aeolian transport). Vegetations are commonly regarded as scalar-moderation structures in urban areas by influencing airflows, avoiding aeolian erosion and pollution. The morphological multi-scale branch structures of the single tree were usually neglected in previous studies, which, however, can lead to highly inhomogeneous wake. The turbulence generation through the multi-scale branches is not yet clear. The corresponding spatial details of scalar distribution and its mechanism need to be further examined. In this study, the wake after a single fractal tree is investigated using the large-eddy simulation (LES). The spatial characteristics of the wakes are illustrated, indicating susceptible zones in the canopy, first-generation levels of branches, and recirculation regions. Higher-order moments signify the rare, extreme events with strong upward acceleration, removing air masses likely to occur at the canopy shear layer within 3h downstream where h is the characteristic tree height. The quadrant analysis elucidates the scalar transport potential. The scalar distribution and its transport mechanism are visualized, showing dilution and accumulation by the streamwise and vertical velocities. Involving the multi-scale geometries of the tree can provide guidance for urban green infrastructure (UGI) planning, advancing environmental and sustainability development.

High-order discretization–based self-adaptive turbulence eddy simulation for supersonic base flow with PHengLEI software

PurposeThe purpose of the present study is to develop a new numerical framework that can predict the supersonic base flow more accurately, including the development of axisymmetrically separated shear layer and recompression shock. To this end, two aspects are improved and combined, i.e. a newly self-adaptive turbulence eddy simulation (SATES) turbulence modeling method and a high-order discretization numerical scheme. Furthermore, the performance of the new numerical framework within a general-purpose PHengLEI software is assessed in detail.Design/methodology/approachSatisfactory prediction of the supersonic separated shear layer with unsteady wake flow is quite challenging. By using a unified turbulence model called SATES combining high-order accurate discretization numerical schemes, the present study first assesses the performance of newly developed SATES for supersonic axisymmetric separation flows. A high-order finite differencing-based compressible computational fluid dynamics (CFD) code called PHengLEI is developed and several different numerical schemes are used to investigate the effects on shock-turbulence interactions, which include the monotonic upstream-centered scheme for conservation laws (MUSCL), weighted compact nonlinear scheme (WCNS) and hybrid cell-edge and cell-node dissipative compact scheme (HDCS).FindingsCompared with the available experimental data and the numerical predictions, the results of SATES by using high-order accurate WCNS or HDCS schemes agree better with the experiments than the results by using the MUSCL scheme. The WCNS and HDCS can also significantly improve the prediction of flow physics in terms of the instability of the annular shear layer and the evolution of the turbulent wake.Research limitations/implicationsThe small deviations in the recirculation region can be found between the present numerical results and experimental data, which could be caused by the inaccurate incoming boundary layer condition and compressible effects. Therefore, a proper incoming boundary layer condition with turbulent fluctuations and compressibility effects need to be considered to further improve the accuracy of simulations.Practical implicationsThe present study evaluates a high-order discretization-based SATES turbulence model for supersonic separation flows, which is quite valuable for improving the calculation accuracy of aeronautics applications, especially in supersonic conditions.Originality/valueFor the first time, the newly developed SATES turbulence modeling method combining the high-order accurate WCNS or HDCS numerical schemes is implemented on the PHengLEI software and successfully applied for the simulations of supersonic separation flows, and satisfactory results are obtained. The unsteady evolutions of the supersonic annular shear layer are analyzed, and the hairpin vortex structures are found in the simulation.

Spray characteristics of different regions downstream of a swirl cup

The spray characteristics of different regions downstream of swirl cups play a critical role in cold start and re-ignition of gas turbines. The spray measurements were performed at the fuel pressures of 0.5, 0.8, 1.0, 1.5, and 2.0 MPa and the fuel temperatures of −23, −13, −3, 7, 17 and 27 ℃, respectively. The droplet size, droplet velocity, droplet number, and instantaneous spatial spray image of sprays from an aviation kerosene Jet-A were measured using a two-component phase Doppler particle analyzer and a digital off-axis holography system. As the fuel pressure and temperature increase, the Sauter Mean Diameter (SMD) and spray non-uniformity of the Spray Shear Layer (SSL) gradually decrease. As the fuel pressure increases, the SMD and spray non-uniformity of the Central Toroidal Recirculation Zone (CTRZ) gradually decrease, and the slopes of these curves both decrease. As the fuel pressure increases, the SMD and spray non-uniformity of the CTRZ rapidly decrease at the fuel temperature of −23 ℃, while slightly decrease at the fuel temperature of 27 ℃. The droplets in the CTRZ come from 3 different sources: simplex nozzle, venturi, and outside the CTRZ. As the fuel pressure increases, the proportion of droplets recirculated from outside the CTRZ decreases. This study proposed the concept of the “pressure critical point” for the swirl cups. As the fuel temperature decreases, the proportion of droplets recirculated from outside the CTRZ increases below the critical pressure, while decreases above the critical pressure. In addition, through the models of liquid film formation and breakup on the curved cylindrical wall, a semi-theoretical model was established to predict the SMD of SSL for swirl cups. The prediction uncertainty of this model is less than 6% for all 14 conditions in this paper.

Numerical investigation on the transverse jet into a supersonic crossflow with different pressure ratios

The jet can be applied to the yaw control of a projectile, however, the complex interaction of the jet with the supersonic mainstream makes the flow field complex and yaw force unpredictable. To reveal the evolution of flow structures under different pressure ratios (PRs), or momentum flux ratios, a transverse sonic jet injected into a supersonic laminar crossflow has been studied numerically. Large-eddy simulations are employed to simulate the flow fields and evolution tendency of flow structures under different PRs of 10, 50, 100, 300, and 500. Our results show clearly the shock and flow structures of the jet interaction with crossflow under different PRs. Moreover, we find that, with the increase of PR, a larger upstream recirculation zone (RZ) and jet shock core appear, which accelerates the transformation of the bow shock (BoS) and the instability of the jet shear layer due to its stronger interaction with the crossflow. In addition, a high PR also accelerates RZ instability and produces a strong compressing effect on the major counter-rotating vortex pair in the jet flow, which makes the streamwise vortex tube stronger in the wake. These findings provide important information for applications of jet control of projectiles.

Aerodynamic flow field analysis of vortex interactions on multi-swept wing configurations

The underlying flow physics of vortex breakdown and vortex-vortex interactions at moderate to high angles of attack is not well understood. This lack of understanding is a challenge to the development of optimized wing designs for fighter aircraft. The goal of this study is to gain insight into the breakdown of vortices on multi-swept wing configurations. Two generic models were considered, one with a gradual vortex breakdown and one with a rapid (abrupt) breakdown. Aerodynamic forces were calculated from the US Air Force Academy (USAFA) subsonic wind tunnel to validate the experimental setup using computational fluid dynamics (CFD) and previous work. Stereoscopic particle imaging velocimetry (SPIV) was used to visualize vortex-vortex interactions, measure vortex strengths, and identify vortex breakdown locations. The measured lift and drag forces showed flow physics that agreed with previous experiments and CFD, supporting the validity of the new experimental setup. Analysis of SPIV and CFD data improved our understanding of the different types of vortex interactions observed for both models. The only difference between the models was the location of the secondary sweep angle on the chord. The C4L80 model, with a secondary sweep angle at 40% of the root chord, had a large shear layer along the leading edge, which led to vortex merging and gradual vortex breakdown at lower angles of attack than its counterpart. The C6L80 model, with a secondary sweep angle at 60% of the root chord, had no such shear layer and instead developed strong, isolated vortices that were characterized by vortex braiding and abrupt breakdown at high angles of attack. These vortices were entrained into the main vortex at the leading edge of the strake (inboard vortex), creating a balanced vortex system that enhanced lift and delayed the onset of gradual vortex breakdown seen in the C4L80 model.

Exploring the influence of vegetated mid-channel bar on flow and turbulence in bifurcated channels: An experimental approach

The evolution of fluvial systems is greatly impacted by mid-channel bars, a typical morphodynamic process in natural rivers. Sometimes, the growth of vegetation over these bars complicates the morphological behaviour by interacting with the flow. It is therefore necessary to have a fundamental interpretation of the flow-turbulence structure around the mid-bar in presence of vegetation cover in order to understand braiding dynamics, still studies in this area are scarce. The present study investigates the process-form-vegetation-interaction through experimental investigation at a flume scale mid-channel bar model with different natural vegetation cover arrangements (paddy, leafy, and rigid stem). The flow-turbulence behaviour has been observed through the bifurcated channel using the three-dimensional Acoustic Doppler Velocimeter (ADV). Results showed that the longitudinal velocity component varies with the different vegetation cover, and it was highest with leafy vegetation (about 32%). Similarly, the Reynolds Stress and Turbulence Intensity were also observed to be higher in case of leafy vegetation. A unique pattern of flow-turbulence parameters was observed near the bar level, the lower canopy level, and the upper canopy level. Moreover, it was found that vegetation structure and its flexible nature influence both longitudinal velocity reduction and momentum transfer at and over the canopy, as well as the thickness of the shear layer region.

DeHNSSo: The Delft Harmonic Navier-Stokes Solver for nonlinear stability problems with complex geometric features

A nonlinear Harmonic Navier-Stokes (HNS) framework is introduced for simulating instabilities in laminar spanwise-invariant shear layers, featuring sharp and smooth wall surface protuberances. While such cases play a critical role in the process of laminar-to-turbulent transition, classical stability theory analyses such as parabolized or local stability methods fail to provide (accurate) results, due to their underlying assumptions. The generalized incompressible Navier-Stokes (NS) equations are expanded in perturbed form, using a spanwise and temporal Fourier ansatz for flow perturbations. The resulting equations are discretized using spectral collocation in the wall-normal direction and finite-difference methods in the streamwise direction. The equations are then solved using a direct sparse-matrix solver. The nonlinear mode interaction terms are converged iteratively. The solution implementation makes use of a generalized domain transformation to account for geometrical smooth surface features, such as humps. No-slip conditions can be embedded in the interior domain to account for the presence of sharp surface features such as forward- or backward-facing steps. Common difficulties with Navier-Stokes solvers, such as the treatment of the outflow boundary and convergence of nonlinear terms, are considered in detail. The performance of the developed solver is evaluated against several cases of representative boundary layer instability growth, including linear and nonlinear growth of Tollmien-Schlichting waves in a Blasius boundary layer and stationary crossflow instabilities in a swept flat-plate boundary layer. The latter problem is also treated in the presence of a geometrical smooth hump and a sharp forward-facing step at the wall. HNS simulation results, such as perturbation amplitudes, growth rates, and shape functions, are compared to benchmark flow stability analysis methods such as Parabolized Stability Equations (PSE), Adaptive Harmonic Linearized Navier-Stokes (AHLNS), or Direct Numerical Simulations (DNS). Good agreement is observed in all cases. The HNS solver is subjected to a grid convergence study and a simple performance benchmark, namely memory usage and computational cost. The computational cost is found to be considerably lower than high-fidelity DNS at comparable grid resolutions. Program summaryProgram Title: DeHNSSoCPC Library link to program files:https://doi.org/10.17632/9bnms99kk2.1Developer's repository link:https://github.com/SvenWesterbeek/DeHNSSoLicensing provisions: GPLv3Programming language: MatlabSupplementary material: The supplementary material contains the code as well as a user manual.Nature of problem: Fluid flows are subject to laminar-to-turbulent transition following the growth of instabilities. To avoid computationally demanding Direct Numerical Simulations (DNS), perturbation theory is often applied to their analysis. However, classical stability methods based on the Orr-Sommerfeld equation or the Parabolized Stability Equations neglect the influence of streamwise gradients in varying degrees. The validity of these assumptions is difficult to estimate a priori.Solution method: The Delft Harmonic Navier-Stokes Solver (DeHNSSo) solves the harmonic Navier-Stokes equations nonlinearly on domains featuring sharp and smooth spanwise invariant surface features using a generalized grid approach in combination with an embedded boundary method. This allows the user to include the effects of streamwise gradients on flow at a fraction of the cost of DNS.Additional comments including restrictions and unusual features: In DeHNSSo, the equations are solved using direct matrix solvers. As such, memory is treated as a scarce resource. The problem is formulated in a mode-independent manner such that left-hand side matrices need only be computed and stored once despite incorporating nonlinear terms. Additionally, DeHNSSo offers the user the possibility to prescribe inhomogeneous boundary conditions to introduce and solve the receptivity problem or manipulate instabilities. Due to the double Fourier expansion, the solver is restricted to spanwise and temporally periodic problems.