We provide an overview of current techniques and typical applications of numerical bifurcation analysis in fluid dynamical problems. Many of these problems are characterized by high-dimensional dynamical systems which undergo transitions as parameters are changed. The computation of the critical conditions associated with these transitions, popularly referred to as 'tipping points', is important for understanding the transition mechanisms. We describe the two basic classes of methods of numerical bifurcation analysis, which differ in the explicit or implicit use of the Jacobian matrix of the dynamical system. The numerical challenges involved in both methods are mentioned and possible solutions to current bottlenecks are given. To demonstrate that numerical bifurcation techniques are not restricted to relatively low-dimensional dynamical systems, we provide several examples of the application of the modern techniques to a diverse set of fluid mechanical problems.

We present a new numerical method to approximate the solutions of an Euler-Poisson model, which is inherent to astrophysical flows where gravity plays an important role. We propose a discretization of gravity which ensures adequate coupling of the Poisson and Euler equations, paying particular attention to the gravity source term involved in the latter equations. In order to approximate this source term, its discretization is introduced into the approximate Riemann solver used for the Euler equations. A relaxation scheme is involved and its robustness is established. The method has been implemented in the software HERACLES [29] and several numerical experiments involving gravitational flows for astrophysics highlight the scheme.

In this study, we present a new numerical model for crystal growth in a vertical solidification system. This model takes into account the buoyancy induced convective flow and its effect on the crystal growth process. The evolution of the crystal growth interface is simulated using the phase-field method. A semi-implicit lattice kinetics solver based on the Boltzmann equation is employed to model the unsteady incompressible flow. This model is used to investigate the effect of furnace operational conditions on crystal growth interface profiles and growth velocities. For a simple case of macroscopic radial growth, the phase-field model is validated against an analytical solution. The numerical simulations reveal that for a certain set of temperature boundary conditions, the heat transport in the melt near the phase interface is diffusion dominant and advection is suppressed.

The conventional Poisson-Nernst-Planck equations do not account for the finite size of ions explicitly. This leads to solutions featuring unrealistically high ionic concentrations in the regions subject to external potentials, in particular, near highly charged surfaces. A modified form of the Poisson-Nernst-Planck equations accounts for steric effects and results in solutions with finite ion concentrations. Here, we evaluate numerical methods for solving the modified Poisson-Nernst-Planck equations by modeling electric field-driven transport of ions through a nanopore. We describe a novel, robust finite element solver that combines the applications of the Newton's method to the nonlinear Galerkin form of the equations, augmented with stabilization terms to appropriately handle the drift-diffusion processes.
To make direct comparison with particle-based simulations possible, our method is specifically designed to produce solutions under periodic boundary conditions and to conserve the number of ions in the solution domain. We test our finite element solver on a set of challenging numerical experiments that include calculations of the ion distribution in a volume confined between two charged plates, calculations of the ionic current though a nanopore subject to an external electric field, and modeling the effect of a DNA molecule on the ion concentration and nanopore current.

We focus on the problem of evaluating the velocity field outside a solid object moving in an incompressible Stokes flow using the boundary integral formulation. For points near the boundary, the integral is nearly singular, and accurate computation of the velocity is not routine. One way to overcome this problem is to regularize the integral kernel. The method of regularized Stokeslet (MRS) is a systematic way to regularize the kernel in this situation. For a specific blob function which is widely used, the error of the MRS is only of first order with respect to the blob parameter. We prove that this is the case for radial blob functions with decay property $ϕ(r)$=$\mathcal{O}$($r^{−3−α}$) when $r→∞$ for some constant $α>1$. We then find a class of blob functions for which the leading local error term can be removed to get second and third order errors with respect to blob parameter. Since the addition of these terms might give a flow field that is not divergence free, we introduce a modification of these terms to make the divergence of the corrected flow field close to zero while keeping the desired accuracy. Furthermore, these dominant terms are explicitly expressed in terms of blob function and so the computation time is negligible.

Assumed having axial symmetry, the streamer discharge is often described by a fluid model in cylindrical coordinate system, which consists of convection dominated (diffusion) equations with source terms, coupled with a Poisson's equation. Without additional care for a stricter CFL condition or special treatment to the negative source term, popular methods used in streamer discharge simulations, e.g., FEM-FCT, FVM, cannot ensure the positivity of the particle densities for the cases in attaching gases. By introducing the positivity-preserving limiter proposed by Zhang and Shu [15] and Strang operator splitting, this paper proposes a finite difference scheme with a provable positivity-preserving property in cylindrical coordinate system, for the numerical simulation of streamer discharges in non-attaching and attaching gases. Numerical examples in non-attaching gas (N₂) and attaching gas (SF₆) are given to illustrate the effectiveness of the scheme.

We investigate connections between nonlocal continuum models and molecular dynamics. A continuous upscaling of molecular dynamics models of the form of the embedded-atom model is presented, providing means for simulating molecular dynamics systems at greatly reduced cost. Results are presented for structured and structureless material models, supported by computational experiments. The nonlocal continuum models are shown to be instances of the state-based peridynamics theory. Connections relating multibody peridynamic models and upscaled nonlocal continuum models are derived.

Orographic drag formation is investigated using a numerical wave model (NWM), based on the pressure-coordinate dynamics of non-hydrostatic HIRLAM. The surface drag, wave stress (vertical flux of horizontal momentum), and wave drag are split to the longitudinal and transverse components and presented as Fourier sums of their spectral amplitudes weighted with the power spectrum of relative orographic height. The NWM is accomplished, enabling a spectral investigation of the buoyancy wave stress, and drag generation by orography and is then applied to a cold front, characterised by low static stability of the upper troposphere, large vertical and directional wind variations, and intensive trapped wave generation downstream of obstacles. Resonances are discovered in the stress and drag spectra in the form of high narrow peaks. The stress conservation problem is revisited. Longitudinal stress conserves in unidirectional flow, 2D orography conditions, but becomes convergent for rotating wind or 3D orography. Even in the convergent case the vertical momentum flux from the troposphere to stratosphere remains substantial. The transverse stress never conserves. Disappearing at the surface and on the top, it realises the main momentum exchange between lower an upper parts of the troposphere. Existence of stationary stratospheric quasi-turbulence (SQT) is established above wind minimum in the stratosphere.

Aiming at simulating micro gas flows with accurate boundary conditions, an efficient hybrid algorithm is developed by combining the molecular dynamics (MD) method with the direct simulation Monte Carlo (DSMC) method. The efficiency comes from the fact that the MD method is applied only within the gas-wall interaction layer, characterized by the cut-off distance of the gas-solid interaction potential, to resolve accurately the gas-wall interaction process, while the DSMC method is employed in the remaining portion of the flow field to efficiently simulate rarefied gas transport outside the gas-wall interaction layer. A unique feature about the present scheme is that the coupling between the two methods is realized by matching the molecular velocity distribution function at the DSMC/MD interface, hence there is no need for one-to-one mapping between a MD gas molecule and a DSMC simulation particle. Further improvement in efficiency is achieved by taking advantage of gas rarefaction inside the gas-wall interaction layer and by employing the "smart-wall model" proposed by Barisik et al. The developed hybrid algorithm is validated on two classical benchmarks namely 1-D Fourier thermal problem and Couette shear flow problem. Both the accuracy and efficiency of the hybrid algorithm are discussed. As an application, the hybrid algorithm is employed to simulate thermal transpiration coefficient in the free-molecule regime for a system with atomically smooth surface. Result is utilized to validate the coefficients calculated from the pure DSMC simulation with Maxwell and Cercignani-Lampis gas-wall interaction models.

In this paper, the power-law fluid flows in a two-dimensional square cavity are investigated in detail with multi-relaxation-time lattice Boltzmann method (MRT-LBM). The influence of the Reynolds number (Re) and the power-law index (n) on the vortex strength, vortex position and velocity distribution are extensively studied. In our numerical simulations, Re is varied from 100 to 10000, and n is ranged from 0.25 to 1.75, covering both cases of shear-thinning and shear-thickening. Compared with the Newtonian fluid, numerical results show that the flow structure and number of vortex of power-law fluid are not only dependent on the Reynolds number, but also related to power-law index.