Turbulence Closure Model Research Articles

Numerical simulation on the influence of artificial island on reef hydrodynamics

The construction of artificial island can greatly change the reef hydrodynamics, leading to increased water level and wave height as well as changes in flow field distribution. These alterations can affect sediment transport on the reef, and increase the risk of overtopping and structure instability. A numerical model based on the Reynolds Averaged Navier-Stokes (RANS) equations and k-ε turbulence closure model was developed to investigate the influence of artificial island on reef hydrodynamics. The numerical model was validated against the experimental results of wave height, mean water level, and wave breaking morphology. Detailed flow field, wave height, and wave set-up in front of artificial island were further analyzed based on the validated model. After building the reef-top structure, the wave breaking and offshore currents at reef edge were amplified. The flow stratification and increase of the wave set-up were also found on the reef flat. Furthermore, we found the relationship between maximum flow velocities on the reef flat and incoming wave conditions could be characterized by two non-dimensional parameters: |u¯|max/g(η¯+hr), (η¯+hr)/Hi.

Spatially Local Surrogate Modeling of Subgrid-Scale Effects in Idealized Atmospheric Flows: A Deep Learned Approach Using High-Resolution Simulation Data

Abstract We introduce a machine learned surrogate model from high-resolution simulation data to capture the subgrid-scale effects in dry, stratified atmospheric flows. We use deep neural networks (NNs) to model the spatially local state differences between a coarse-resolution simulation and a high-resolution simulation. The setup enables the capture of both dissipative and antidissipative effects in the state differences. The NN model is able to accurately capture the state differences in offline tests outside the training regime. In online tests intended for production use, the NN-coupled coarse simulation has higher accuracy over a significant period of time compared to the coarse-resolution simulation without any correction. We provide evidence of the capability of the NN model to accurately capture high-gradient regions in the flow field. With the accumulation of the errors, the NN-coupled simulation becomes computationally unstable after approximately 90 coarse simulation time steps. Insights gained from these surrogate models further pave the way for formulating stable, complex, physics-based spatially local NN models which are driven by traditional subgrid-scale turbulence closure models. Significance Statement Flows in the atmosphere are highly chaotic and turbulent, comprising flow structures of broad scales. For effective computational modeling of atmospheric flows, the effects of the small- and large-scale structures need to be captured by the simulations. Capturing the small-scale structures requires fine-resolution simulations. Even with the current state-of-the-art supercomputers, it can be prohibitively expensive to simulate these flows when computed for the entire earth over climate time scales. Thus, it is necessary to focus on the larger-scale structures using a coarse-resolution simulation while capturing the effects of the smaller-scale structures using some parameterization (approximation) scheme and incorporating it into the coarse-resolution simulation. We use machine learning to model the effects of the small-scale structures (subgrid-scale effects) in atmospheric flows. Data from a fine-resolution simulation is used to compute the missing subgrid-scale effects in coarse-resolution simulations. We then use machine learning models to approximate these differences between the coarse- and fine-resolution simulations. We see improved accuracy for the coarse-resolution simulations when corrected using these machine learned models.

Generalized field inversion strategies for data-driven turbulence closure modeling

Most data-driven turbulence closures are based on the general structure of nonlinear eddy viscosity models. Although this structure can be embedded into the machine learning algorithm and the Reynolds stress tensor itself can be fit as a function of scalar- and tensor-valued inputs, there exists an alternative two-step approach. First, the spatial distributions of the optimal closure coefficients are computed by solving an inverse problem. Subsequently, these are expressed as functions of solely scalar-valued invariants of the flow field by virtue of an arbitrary regression algorithm. In this paper, we present two general inversion strategies that overcome the limitation of being applicable only when all closure tensors are linearly independent. We propose to either cast the inversion into a constrained and regularized optimization problem or project the anisotropy tensor onto a set of previously orthogonalized closure tensors. Using the two-step approach together with either of these strategies then enables us to quantify the model-form error associated with the closure structure independent of a particular regression algorithm. Eventually, this allows for the selection of the a priori optimal set of closure tensors for a given, arbitrary complex test case.

Hydrodynamic force coefficients for spherical triangle shell fragments: Dependence on the aspect ratio and flatness

Euler–Lagrange simulations of particle-laden flow require hydrodynamic models of drag and lift forces for individual particles. Our goal is to develop models that can prescribe these forces for arbitrarily orientated shell objects. Here, we use computational fluid dynamics simulations of steady bottom-boundary layer flow over a series of spherical triangle shell fragments to calculate the hydrodynamic forces. The simulations explicitly resolve the wall boundary layers using grid resolution on the order of y+=1 at the shell fragment surface and use the SST k-omega turbulence closure model. These fragments cover a range of aspect ratio and flatness characteristics. The shell fragments are generated as triangular selections of a spherical shell with azimuthal and longitudinal angles proscribed based on elongation and flatness parameters (varying between 1 to 5, and 0.02 to 0.2 respectively), while characteristic length of the fragment is held constant to define the overall fragment size. Fragment orientations are considered with independently varying pitch, roll, and yaw each ranging from 0 to 180 degrees. The numerical estimates for the forces from all simulations were used to develop robust parameterizations of the drag and lift as a function of aspect ratio and flatness characteristics, as well as orientation of the shell fragments.

Response of Upper Ocean to Parameterized Schemes of Wave Breaking under Typhoon Condition

The study of upper ocean mixing processes, including their dynamics and thermodynamics, has been a primary focus for oceanographers and meteorologists. Wave breaking in deep water is believed to play a significant role in these processes, affecting air–sea interactions and contributing to the energy dissipation of surface waves. This, in turn, enhances the transfer of gas, heat, and mass at the ocean surface. In this paper, we use the FVCOM-SWAVE coupled wave and current model, which is based on the MY-2.5 turbulent closure model, to examine the response of upper ocean turbulent kinetic energy (TKE) and temperature to various wave breaking parametric schemes. We propose a new parametric scheme for wave breaking energy at the sea surface, which is based on the correlation between breaking wave parameter RB and whitecap coverage. The impact of this new wave breaking parametric scheme on the upper ocean under typhoon conditions is analyzed by comparing it with the original parametric scheme that is primarily influenced by wave age. The wave field simulated by SWAVE was verified using Jason-3 satellite altimeter data, confirming the effectiveness of the simulation. The simulation results for upper ocean temperature were also validated using OISST data and Argo float observational data. Our findings indicate that, under the influence of Typhoon Nanmadol, both parametric schemes can transfer the energy of sea surface wave breaking into the seawater. The new wave breaking parameter RB scheme effectively enhances turbulent mixing at the ocean surface, leading to a decrease in sea surface temperature (SST) and an increase in mixed layer depth (MLD). This further improves upon the issue of uneven mixing of seawater at the air–sea interface in the MY-2.5 turbulent closure model. However, it is important to note that wave breaking under typhoon conditions is only one aspect of wave impact on ocean disturbances. Therefore, further research is needed to fully understand the impact of waves on upper ocean mixing, including the consideration of other wave mechanisms.

Numerical investigation of premixed syngas/air flame evolution in a closed duct with tulip flame formation

In this study, the dynamic behavior of the premixed syngas/air flames in enclosed ducts with obstacles at the blockage ratio of 0.5, as well as without obstacles, was studied via large eddy simulation and the turbulent flame speed closure model. The hydrogen volume fraction (hydrogen content in fuel) was set to 50 %. Results indicate that the simulation reproduced the flame shapes, flame speed and explosion overpressure observed in the experiments. In cases without obstacles, a “tulip” flame with a tulip cusp can be seen, and secondary cusps emerge on tulip lips due to the high-pressure zone ahead of “tulip” flame lips, leading to in a distorted “tulip” flame. In cases with obstacles, the flame develops hemispherical and finger shapes upstream of the obstacles and evolves into a “tulip” flame downstream of the obstacles due to vortices and adverse pressure gradients. The alteration of the flow direction within the unburned and burned zones results in intricate flame shapes. Flame speed and overpressure show a close relationship and exhibit oscillations following the generation of the “tulip” flame. Typically, the velocity, pressure and flow fields adjacent to the flame front can affect dynamic behaviors of premixed syngas/air flames.

Grey‐zone simulations of shallow‐to‐deep convection transition using dynamic subgrid‐scale turbulence models

AbstractWe examine the ability of two dynamic turbulence closure models to simulate the diurnal development of convection and the transition from dry to shallow cumuli and then to deep convection. The dynamic models are compared with the conventional Smagorinsky scheme at a range of cloud‐resolving and grey‐zone resolutions. The dynamic schemes include the Lagrangian‐averaged, scale‐dependent dynamic Smagorinsky model and a Lagrangian‐averaged, dynamic mixed model. The conventional Smagorinsky model fails to reproduce the shallow convection stage beyond the large‐eddy simulation regime, continuously building up the convective available potential energy that eventually leads to an unrealistic deep convection phase. The dynamic Smagorinsky model significantly improves the representation of shallow and deep convection; however, it exhibits issues similar to the conventional scheme at coarser resolutions. In contrast, the dynamic mixed model closely follows the large‐eddy simulation results across the range of sub‐kilometre simulations. This is achieved by the combined effect of an adaptive length scale and the inclusion of the Leonard terms, which can produce counter‐gradient fluxes through the backscatter of energy from the subgrid to the resolved scales and enable appropriate non‐local contributions. A further sensitivity test on the inclusion of the Leonard terms on all hydrometeor fluxes reveals the strong interaction between turbulent transport and microphysics and the possible need for further optimisation of the dynamic mixed model coefficients together with the microphysical representation.

Pressure Loss Modeling for Multi-Stage Obstacles in Pressurized Ducts

Estimating singular pressure losses for multi-stage obstacles in pressurized hydraulic ducts is a challenging task. An experimental study was conducted in a closed-loop hydrodynamic tunnel to characterize the pressure losses of a system consisting of a porous fibrous foam placed in front of a bar rack. The pressure losses of different foam–rack configurations were measured over a range of inlet velocities in order to highlight the mutual influence of their characteristics on the flow. The interdependence between the two stages has been evidenced by both the experimental results and additional numerical simulations using RANS (Reynolds-Averaged Navier–Stokes Equations) simulations with a k-ω SST turbulent closure model. The pressure losses were first modeled using two approaches based on the assumption of either independence or full dependence between the stages. The respective advantages and limitations of these approaches led to an improved analytical formula that considers the transition of the flow from the porous foam to the bar rack. By taking into account an empirical transition factor, the proposed model improves the head loss prediction for all tested configurations, with an average relative error between the formula and experimental results less than that of the two simpler approaches. This study improves our understanding of global pressure losses in multi-stage systems that include a porous foam or other filtering or clogging media in front of bar racks.

Vertical Energy Fluxes Driven by the Interaction between Wave Groups and Langmuir Turbulence

Abstract Data from an air–sea interaction tower are used to close the turbulent kinetic energy (TKE) budget in the wave-affected surface layer of the upper ocean. Under energetic wind forcing with active wave breaking, the dominant balance is between the dissipation rate of TKE and the downward convergence in vertical energy flux. The downward energy flux is driven by pressure work, and the TKE transport is upward, opposite to the downgradient assumption in most turbulence closure models. The sign and the relative magnitude of these energy fluxes are hypothesized to be driven by an interaction between the vertical velocity of Langmuir circulation (LC) and the kinetic energy and pressure of wave groups, which is the result of small-scale wave–current interaction. Consistent with previous modeling studies, the data suggest that the horizontal velocity anomaly associated with LC refracts wave energy away from downwelling regions and into upwelling regions, resulting in negative covariance between the vertical velocity of LC and the pressure anomaly associated with the wave groups. The asymmetry between downward pressure work and upward TKE flux is explained by the Bernoulli response of the sea surface, which results in groups of waves having a larger pressure anomaly than the corresponding kinetic energy anomaly, consistent with group-bound long-wave theory.

Wind characteristics around a skyway bridge of high-rise buildings

To respond the expansion of urban centers, the proliferation of high-rise buildings demands a better understanding of the aerodynamic phenomena around skyway bridges connecting these structures. This analysis, utilizing the advanced computational fluid dynamics verified by wind tunnel test data, investigates the wind characteristics around such bridges, crucial for structural stability, pedestrian comfort, and aerodynamic efficiency. This study focuses on the interactions between a 2 × 2 building array with a building height-to-street width ratio of 30 and a skyway bridge, investigating those factors such as bridge influence, building structures, building height, and bridge position. Using the three-dimensional steady Reynolds-averaged Navier–Stokes equations along with the Reynolds stress model for turbulence closure, the results show that the presence of skyway bridge significantly modifies local wind patterns. Wind speed and turbulence intensity are impacted differently based on the bridge's upstream or downstream settings. Downstream bridges tend to reduce wind speeds due to the sheltering effects, while upstream placement of bridge can enhance wind flow, affecting both the structural design and pedestrian comfort. Additionally, building height variations adjacent to the bridge influence wind velocity and pressure profiles, with taller buildings intensifying wind speeds at lower levels because of the channeling effects. These insights are pivotal for optimizing the skyway bridge designs to improve airflow distribution, enhance environmental sustainability, and ease wind-caused disturbances, offering a guideline for future architectural and urban planning in high-rise districts.

Recent advances and effectiveness of machine learning models for fluid dynamics in the built environment

ABSTRACT Indoor environmental quality is crucial for human health and comfort, necessitating precise and efficient computational methods to optimise indoor climate parameters. Recent advancements in machine learning (ML) and computational fluid dynamics (CFD) are promising. However, applying ML to complex building airflow presents challenges. This research aims to investigate the integration of ML with CFD in the context of built environment applications using a systematic review approach. It highlights a critical knowledge gap: the need to synthesise innovative approaches that address the limitations of indoor modelling using data-driven ML methods. The review examines contemporary literature, identifying current developments and suggesting potential future directions. It delves into the innovations in combining ML with CFD to predict thermal comfort and indoor air quality, uncovering key limitations such as the lack of high-quality experimental data for training and validation, the computational complexity of detailed CFD simulations, and the interpretability issues of ‘black-box’ ML models. The emergence of data-driven techniques in fluid mechanics offers promising prospects for modelling in the built environment. Future research should focus on incorporating physics-based rules in ML models, adapting turbulence closure models for indoor flows, and enhancing model validation using real-world datasets. The research emphasises the synergistic relationship between ML and CFD; it proposes pathways to overcome current limitations, aiming to enhance the precision and efficiency of indoor environment modelling through their integration.

An Improved Non‐Local Planetary Boundary Layer Parameterization Scheme in Weather Forecasting and Research Model Based on a 1.5‐Order Turbulence Closure Model

AbstractPlanetary boundary layer (PBL) modeling is a primary contributor to uncertainties in a numerical weather prediction (NWP) model due to difficulties in modeling the turbulent transport of surface fluxes. The Weather Research and Forecasting model (WRF) has provided many PBL schemes that may feature a non‐local transport component driven by large eddies or a one‐and‐half order turbulence closure model, but few of them possess the two features at once. In the present study, a turbulent kinetic energy (TKE)‐based eddy diffusivity/viscosity method is integrated into the non‐local Asymmetric Convective Model version 2 (ACM2) PBL scheme and implemented in WRF. The original first‐order eddy‐diffusivity term in ACM2 is discarded and an extra prognostic equation for TKE, which considers the tendency of TKE by buoyancy, wind shear, vertical transport, and dissipative processes, is supplied to calculate the diffusivity/viscosity. Non‐local transport is modeled the same as ACM2 using the transilient matrix method. Idealized tests using prescribed surface heat flux and roughness length are performed. TKE‐ACM2 displays advantages over the PBL scheme developed by Bougeault and Lacarrère (hereinafter referred to as Boulac) and ACM2 in the wind speeds (WS) profile because it better matches large‐eddy simulations results in the surface momentum flux. Real case simulations show that TKE‐ACM2 generally outperforms in the diurnal vertical profiles of WS under stable conditions. TKE‐ACM2 also produces a better alignment under moderately unstable conditions in the early nighttime at the urban LiDAR station. However, the model exhibits discrepancies more apparently under a more unstable condition during the winter daytime.

Numerical and experimental analysis of the formation of nitrogen oxides in a non-premixed industrial gas burner

Nitrogen oxides (NOx) formation in a non-premixed industrial gas burner is analyzed with Reynolds-averaged Navier-Stokes (RANS) computational fluid dynamics, using the Flamelet Generated Manifold (FGM) approach for chemistry and the Reynolds stress model for turbulence closure. Analysis and validation of the numerical methodology use both the RANS and the Large Eddy Simulation methods to study the Sandia flame D experiment. Four regimes of the gas burner, differing by thermal power, are studied experimentally and numerically, comparing measured and computed flue gas temperature and emissions. The flow features leading to NOx formation are analyzed. At low power, most NOx formation occurs in the flame core, where temperatures are highest. At high power, a significant amount of NOx forms where a secondary stream of oxidizer meets the high-temperature burnt gas stream. The NOx formation rate is modeled either by summing the contributions of the main paths relevant to the studied problem or via two variants of the FGM scalar transport approach. One of these variants features a correction to the transported scalar approach that enables accounting for the effects of non-adiabatic phenomena, e.g. heat radiation, on NOx formation without increasing the look-up table dimensionality and the required memory usage. The correction is shown to improve emission predictions over the baseline model, although quantitative differences between predicted and measured emission levels remain significant. The approach computing the NOx formation rate as by summing contributions yields the best agreement with measurements when the thermal path predominates, and the formation of nitrogen dioxide is negligible.

Turbo-RANS: straightforward and efficient Bayesian optimization of turbulence model coefficients

PurposeIndustrial simulations of turbulent flows often rely on Reynolds-averaged Navier-Stokes (RANS) turbulence models, which contain numerous closure coefficients that need to be calibrated. This paper aims to address this issue by proposing a semi-automated calibration of these coefficients using a new framework (referred to as turbo-RANS) based on Bayesian optimization.Design/methodology/approachThe authors introduce the generalized error and default coefficient preference (GEDCP) objective function, which can be used with integral, sparse or dense reference data for the purpose of calibrating RANS turbulence closure model coefficients. Then, the authors describe a Bayesian optimization-based algorithm for conducting the calibration of these model coefficients. An in-depth hyperparameter tuning study is conducted to recommend efficient settings for the turbo-RANS optimization procedure.FindingsThe authors demonstrate that the performance of the k-ω shear stress transport (SST) and generalized k-ω (GEKO) turbulence models can be efficiently improved via turbo-RANS, for three example cases: predicting the lift coefficient of an airfoil; predicting the velocity and turbulent kinetic energy fields for a separated flow; and, predicting the wall pressure coefficient distribution for flow through a converging-diverging channel.Originality/valueTo the best of the authors’ knowledge, this work is the first to propose and provide an open-source black-box calibration procedure for turbulence model coefficients based on Bayesian optimization. The authors propose a data-flexible objective function for the calibration target. The open-source implementation of the turbo-RANS framework includes OpenFOAM, Ansys Fluent, STAR-CCM+ and solver-agnostic templates for user application.

Parameterized Internal Wave Mixing in Three Ocean General Circulation Models

AbstractThe non‐local model of mixing based on internal wave breaking, IDEMIX, is implemented as an enhancement of a turbulent kinetic energy closure model in three non‐eddy resolving general circulation ocean models that differ in the discretization and choice of computational grids. In IDEMIX internal wave energy is generated by an energy flux resulting from near‐inertial waves induced by wind forcing at the surface, and at the bottom, by an energy flux that parameterizes the transfer of energy between baroclinic and barotropic tides. In all model simulations with IDEMIX, the mixing work is increased compared to the reference solutions without IDEMIX, reaching values in better agreement with finestructure observations. Furthermore, the horizontal structure of the mixing work is more realistic as a consequence of the heterogeneous forcing functions. All models with IDEMIX simulate deeper thermocline depths related to stronger shallow overturning cells in the Indo‐Pacific. In the North Atlantic, deeper mixed layers in simulations with IDEMIX are associated with an increased Atlantic overturning circulation and an increase of northward heat transports toward more realistic values. The response of the deep Indo‐Pacific overturning circulation and the weak bottom cell of the Atlantic to the inclusion of IDEMIX is incoherent between the models, suggesting that additional unidentified processes and numerical mixing may confound the analysis. Applying different tidal forcing functions leads to simulation differences that are small compared to differences between the different models or between simulations with IDEMIX and without IDEMIX.

On the Water Jet Quality at Part-Load Operation of Pelton turbines

Abstract Pelton turbines with multiple nozzles are characterised by their excellent performance over a wide operating range of mass flow rates. Thus, Pelton turbine manifolds must also be designed to deliver good flow quality to the nozzles at part-load operating conditions to accomplish the expectations. A reduced number of nozzles are opened at low mass flow rate operating conditions, where the choice of the open nozzles depends on the water jet quality generated at the individual nozzles. We investigate the flow phenomena occurring in the manifold with the operation of a reduced nozzle number employing numerical simulations to analyse beneficial nozzle opening selections. The turbulent two-phase flow is predicted by the Volume-Of-Fluid (VOF) approach and the k–ω SST turbulence closure model. Different nozzle operating scenarios are simulated, and the implications on the water jet quality, i.e. deviation and shape deformation, are reported. The results show that the free water jet quality criteria can change drastically at individual nozzle operating conditions because of the velocity changes at the branch line junctions.

Experimental and numerical investigations of the water surface profile and wave extrema of supercritical flows in a narrow channel bend

Supercritical flows in channel bends, e.g., in steep streams, chute spillways, and flood and sediment bypass tunnels (SBTs), experience cross-waves, which undulate the free surface. The designs of these hydraulic structures and flood protection retaining structures in streams necessitate computing the locations and water depths of the wave extrema. This study numerically and experimentally investigates the water surface profiles along the sidewalls, the wave extrema flow depths, and their angular locations in a narrow channel bend model of the Solis SBT in Switzerland. The 0.2 m wide and 16.75 m long channel has a bend of 6.59 m radius and 46.5° angle of deviation. The tested flow conditions produced Froude numbers ≈ 2 and aspect ratios ranging from 1.14 to 1.83. Two-phase flow simulations were performed in OpenFOAM using the RNG k–ε turbulence closure model and the volume-of-fluid method. The simulated angular locations of the first wave extrema and the corresponding flow depths deviate marginally, within ± 6.3% and ± 2.1%, respectively, from the experimental observations, which signifies good predictions using the numerical model. Larger deviations, especially for the angular locations of the wave extrema, are observed for the existing analytical and empirical approaches. Therefore, the presented numerical approach is a suitable tool in designing the height of the hydraulic structures with bends and conveying supercritical flows. In the future, the model’s application shall be extended to the design of the height and location of retaining walls, embankments, and levees in steep natural streams with bends.

Numerical study of wave-induced flow and wave set-up on a platform reef with steep slope

A 2D numerical model was developed in this study to investigate the wave-induced flow and wave set-up on a platform reef with steep frontal slope. The model solves the Reynolds Averaged Navier-Stokes (RANS) equations, with the classical k−ε turbulence closure model. The inflow boundary was modified by adding fluid replenishment to balance the water level in the numerical flume to ensure long-term simulation in the platform reef system. The numerical model was first validated by systematically comparing the wave breaking, flow structure, and mean water level against laboratory measurements. The validated model was employed to further analyze the detailed flow field and reveal the influence of the reef morphology on the wave-induced flow and wave set-up. The complex distribution pattern of wave-induced flow and wave set-up on the reef flat were found to be related to the Froude number. The results also reveal that reef slope and reef roughness will not only change the wave breaking morphology but also affect the distribution pattern of wave-induced flow and wave set-up.

Evolutionary neural networks for learning turbulence closure models with explicit expressions

Developing physical closure models with explicit expressions based on a given dataset is essential to science and engineering. For such symbolic regression tasks, biology-inspired evolutionary algorithms are most widely used. However, typical evolutionary algorithms do not utilize any structural information inherent in training data, which limits their performance in finding accurate model structures and coefficients. By combining one evolutionary algorithm, gene expression programing (GEP), with an artificial neural network (ANN) for symbolic regression, we propose a novel evolutionary neural network method, in which candidate expressions are specifically designed so that they can be transformed between the GEP and ANN structures during training iterations. By combining the GEP's global searching and the ANN's gradient optimization capabilities, efficient and robust convergence to accurate models can be achieved. In addition, sparsity-enhancing strategies have been introduced to improve the interpretability of the trained models. The present method has been tested for finding different physical laws and then applied to turbulence modeling problems with different configurations, showing advantages compared to the existing GEP and ANN methods.

CFD Simulations and Phenomenological Modelling of Aerodynamic Stall Hysteresis of NACA 0018 Wing

Computational simulations of three-dimensional flow around a NACA 0018 wing with an aspect ratio of AR=5 were carried out by using the Unsteady Reynolds-Averaged Navier–Stokes (URANS) equations with the Shear-Stress Transport turbulence model closure. Simulations were performed to capture aerodynamic stall hysteresis by using the developed pseudo-transient continuation (PTC) method based on a dual-time step approach in CFD OpenFOAM code. The flow was characterized by incompressible Mach number M=0.12 and moderate Reynolds number Re=0.67×106. The results obtained indicate the presence of noticeable aerodynamic hysteresis in the static dependencies of the force and moment coefficients, as well as the manifestation of bi-stable flow separation patterns, accompanied by the development of asymmetry in the stall zone. The URANS simulation results are in good agreement with the experimental data obtained for the NACA 0018 finite-aspect-ratio wing in the low-speed wind tunnel under the same test conditions. A new phenomenological bifurcation model of aerodynamic stall hysteresis under static and dynamic conditions is formulated and is proven to be able to closely match the experimental data.