Pseudospectral Simulations Research Articles

Localness of energy cascade in hydrodynamic turbulence. I. Smooth coarse graining

We introduce a novel approach to scale decomposition of the fluid kinetic energy (or other quadratic integrals) into band-pass contributions from a series of length scales. Our decomposition is based on a multiscale generalization of the “Germano identity” for smooth, graded filter kernels. We employ this method to derive a budget equation that describes the transfers of turbulent kinetic energy both in space and in scale. It is shown that the interscale energy transfer is dominated by local triadic interactions, assuming only the scaling properties expected in a turbulent inertial range. We derive rigorous upper bounds on the contributions of nonlocal triads, extending the work of Eyink [Physica D 207, 91 (2005)] for low-pass filtering. We also propose a physical explanation of the differing exponents for our rigorous upper bounds and for the scaling predictions of Kraichnan [Phys. Fluids 9, 1728 (1966); J. Fluid Mech. 47, 525 (1971)]. The faster decay predicted by Kraichnan is argued to be the consequence of additional cancellations in the signed contributions to transfer from nonlocal triads after averaging over space. This picture is supported by data from a 5123 pseudospectral simulation of Navier–Stokes turbulence with phase-shift dealiasing.

Large Rayleigh number thermal convection: Heat flux predictions and strongly nonlinear solutions

We investigate the structure of strongly nonlinear Rayleigh–Bénard convection cells in the asymptotic limit of large Rayleigh number and fixed, moderate Prandtl number. Unlike the flows analyzed in prior theoretical studies of infinite Prandtl number convection, our cellular solutions exhibit dynamically inviscid constant-vorticity cores. By solving an integral equation for the cell-edge temperature distribution, we are able to predict, as a function of cell aspect ratio, the value of the core vorticity, details of the flow within the thin boundary layers and rising/falling plumes adjacent to the edges of the convection cell, and, in particular, the bulk heat flux through the layer. The results of our asymptotic analysis are corroborated using full pseudospectral numerical simulations and confirm that the heat flux is maximized for convection cells that are roughly square in cross section.

A spurious evolution of turbulence originated from round-off error in pseudo-spectral simulation

We report a slowly-developing, spurious numerical solution in pseudo-spectral direct numerical simulation (DNS) of incompressible fluid turbulence. When the effect of machine round-off on the divergence-free condition is not carefully controlled, a problem can develop slowly (over about 50 large-eddy turnover times) and eventually leads to an unphysical flow field. The problem was found with a previously published, highly-compact algorithm for pseudo-spectral DNS and therefore it is important to document the contamination of this numerical artifact on simulated turbulence structure and statistics. This is a striking example since the problem is not easily noticeable due to its very long development time, and it does not lead to numerical instability but rather a different flow state. A theory is developed to explain the unphysical evolution and predicts the exponential growth of round-off error induced velocity divergence. The theory shows that any correlation of the large-scale forcing with the velocity field at the beginning of the time step could lead to amplification of the velocity divergence. For this reason, the problem is quite reproducible. Several simple remedies are tested and shown to correct the problem. It is shown that all revised algorithms are identical theoretically to the original algorithm, with the only difference in the level of control for the divergence-free condition of the simulated flow field. A general recommendation is that the pressure projection operation should be performed at the end of each time step to ensure that the divergence-free condition is not contaminated by machine round-off.

Turbulent Channel Flows on a Rotating Earth

Abstract This paper deals with flow in a rectilinear channel on a rotating earth. The flow is directed perpendicular to the background planetary vorticity; both an analytical theory and numerical simulations are employed. The analytical approach assumes the existence of an eddy viscosity and employs a perturbation expansion in powers of the reciprocal of the Rossby number (Ro). At lowest order, a cross-channel circulation arises because of the tilting of the planetary vorticity vector by the shear in the along-channel direction. This circulation causes a surface convergence, which achieves its maximum value at a channel aspect ratio (= width/depth) of approximately 10. The location of the maximum surface convergence moves from near the center of the channel to a position very near the sidewalls as the aspect ratio increases from O(1) to O(100). To include the effects of turbulence, direct numerical pseudospectral simulations of the equations of motion are employed. While holding the friction Reynolds number fixed at 230.27, a series of simulations with increasing rotation (Ro = ∞, 10, 1.0, 0.1) are performed. The channelwide circulation cell observed in the analytical theory occurs for the finite Rossby number, but is displaced by lateral self-advection. In addition, turbulence-driven corner circulations appear, which make the along-channel maximum velocity appear at a subsurface location. The most interesting effect is the segregation of the turbulence to one side of the channel, while the turbulence is suppressed on the opposite side.

FORMATION, SURVIVAL, AND DESTRUCTION OF VORTICES IN ACCRETION DISKS

Two dimensional hydrodynamical disks are nonlinearly unstable to the formation of vortices. Once formed, these vortices essentially survive forever. What happens in three dimensions? We show with pseudospectral simulations that in 3D a vortex in a short box forms and survives just as in 2D. But a vortex in a tall box is unstable and is destroyed. In our simulation, the unstable vortex decays into a transient turbulent-like state that transports angular momentum outward at a nearly constant rate for hundreds of orbital times. The 3D instability that destroys vortices is a generalization of the 2D instability that forms them. We derive the conditions for these nonlinear instabilities to act by calculating the coupling between linear modes, and thereby derive the criterion for a vortex to survive in 3D as it does in 2D: the azimuthal extent of the vortex must be larger than the scale height of the accretion disk. When this criterion is violated, the vortex is unstable and decays. Because vortices are longer in azimuthal than in radial extent by a factor that is inversely proportional to their excess vorticity, a vortex with given radial extent will only survive in a 3D disk if it is sufficiently weak. This counterintuitive result explains why previous 3D simulations always yielded decaying vortices: their vortices were too strong. We conclude that in protoplanetary disks weak vortices can trap dust and serve as the nurseries of planet formation. Decaying strong vortices might be responsible for the outwards transport of angular momentum that is required to make accretion disks accrete.

Strongly nonlinear Langmuir circulation and Rayleigh–Bénard convection

Most rational asymptotic studies of non-rotating Rayleigh–Bénard convection and its cousins have been restricted to the linear or weakly nonlinear regime. An important exception occurs for large Rayleigh-number thermal convection at effectively infinite Prandtl number, i.e. fast but very viscous convection. In this scenario, the temperature field exhibits a layer-like structure surrounding an isothermal core and, crucially, the momentum equation linearizes. These features have been exploited by several authors to obtain semi-analytical nonlinear solutions. AtO(1) Prandtl number, the fluid dynamics in the vortex core is dominated by nonlinear inertial rather than linear viscous effects, substantially altering the vortex structure. Here, it is shown that a combination of matched asymptotic analysis and global conservation constraints can be used to obtain a semi-analytic yet strongly nonlinear description of two related flows: (i) Rayleigh–Bénard convection between constant heat-flux boundaries at unit Prandtl number, and (ii) Langmuir circulation (LC), a wind and wave-driven convective flow commonly observed in natural water bodies. A simple analytical prediction is given for the roll-vortex amplitude, which is shown to be independent of the horizontal wavenumber of the convection pattern. In marked contrast to weakly nonlinear convection cells, the fully nonlinear asymptotic solutions exhibit flow features relevant to turbulent convection including the complete vertical redistribution of the basic-state temperature (or, for LC, downwind velocity) field. Comparisons with well-resolved pseudospectral numerical simulations of the full two-dimensional governing equations confirm the accuracy of the asymptotic results.

Tracking Vector Magnetograms with the Magnetic Induction Equation

The differential affine velocity estimator (DAVE) that we developed in 2006 for estimating velocities from line-of-sight magnetograms is modified to directly incorporate horizontal magnetic fields to produce a differential affine velocity estimator for vector magnetograms (DAVE4VM). The DAVE4VM's performance is demonstrated on the synthetic data from the anelastic pseudospectral ANMHD simulations that were used in the recent comparison of velocity inversion techniques by Welsch and coworkers. The DAVE4VM predicts roughly 95% of the helicity rate and 75% of the power transmitted through the simulation slice. Intercomparison between DAVE4VM and DAVE and further analysis of the DAVE method demonstrates that line-of-sight tracking methods capture the shearing motion of magnetic footpoints but are insensitive to flux emergence - the velocities determined from line-of-sight methods are more consistent with horizontal plasma velocities than with flux transport velocities. These results suggest that previous studies that rely on velocities determined from line-of-sight methods such as the DAVE or local correlation tracking may substantially misrepresent the total helicity rates and power through the photosphere.

A Comparison of Spectral Element and Finite Difference Methods Using Statically Refined Nonconforming Grids for the MHD Island Coalescence Instability Problem

A recently developed spectral-element adaptive refinement incompressible magnetohydrodynamic (MHD) code [Rosenberg, Fournier, Fischer, Pouquet, J. Comp. Phys. 215, 59-80 (2006)] is applied to simulate the problem of MHD island coalescence instability (MICI) in two dimensions. MICI is a fundamental MHD process that can produce sharp current layers and subsequent reconnection and heating in a high-Lundquist number plasma such as the solar corona [Ng and Bhattacharjee, Phys. Plasmas, 5, 4028 (1998)]. Due to the formation of thin current layers, it is highly desirable to use adaptively or statically refined grids to resolve them, and to maintain accuracy at the same time. The output of the spectral-element static adaptive refinement simulations are compared with simulations using a finite difference method on the same refinement grids, and both methods are compared to pseudo-spectral simulations with uniform grids as baselines. It is shown that with the statically refined grids roughly scaling linearly with effective resolution, spectral element runs can maintain accuracy significantly higher than that of the finite difference runs, in some cases achieving close to full spectral accuracy.

Atmospheric Circulation of Close‐in Extrasolar Giant Planets. I. Global, Barotropic, Adiabatic Simulations

We present results from a set of over 300 pseudospectral simulations of atmospheric circulation on extrasolar giant planets with circular orbits. The simulations are of high enough resolution (up to 341 total and sectoral modes) to resolve small-scale eddies and waves, required for reasonable physical accuracy. In this work, we focus on the global circulation pattern that emerges in a shallow, equivalent barotropic, turbulent atmosphere on both tidally synchronized and unsynchronized planets. A full exploration of the large physical and numerical parameter space is performed to identify robust features of the circulation. For some validation, the model is first applied to solar system giant planets. For extrasolar giant planets with physical parameters similar to HD 209458b—a presumably synchronized extrasolar giant planet, representative in many dynamical respects—the circulation is characterized by the following features: (1) a coherent polar vortex that revolves around the pole in each hemisphere; (2) a low number (typically two or three) of slowly varying, broad zonal (east-west) jets that form when the maximum jet speed is comparable to, or somewhat stronger than, those observed on the planets in the solar system; and (3) a motion-associated temperature field, whose detectability and variability depend on the strength of the net heating rate and the global rms wind speed in the atmosphere. In many ways, the global circulation is Earth-like, rather than Jupiter-like. However, if extrasolar giant planets rotate faster and are not close-in (therefore not synchronized), their circulations become more Jupiter-like, for Jupiter-like rotation rates.

Pseudo-spectral time-domain simulation of the transmission and the group delay of photonic devices

The pseudo-spectral time-domain (PSTD) method is applied to the simulation of photonic structures. We study transmission and group delay properties of dielectric and dispersive devices. An auxiliary differential equation technique is used to model the dispersive behavior of metallic features. Comparisons with the finite-difference time-domain method show the superiority of the PSTD method in particular when extracting phase information from a photonic device. The validity of the algorithm is demonstrated by investigating a ring resonator, Bragg gratings and metallic nano-particle arrangements.

Baroclinic Vorticity Production in Protoplanetary Disks. I. Vortex Formation

The formation of vortices in protoplanetary disks is explored via pseudo-spectral numerical simulations of an anelastic-gas model. This model is a coupled set of equations for vorticity and temperature in two dimensions which includes baroclinic vorticity production and radiative cooling. Vortex formation is unambiguously shown to be caused by baroclinicity because (1) these simulations have zero initial perturbation vorticity and a nonzero initial temperature distribution; and (2) turning off the baroclinic term halts vortex formation, as shown by an immediate drop in kinetic energy and vorticity. Vortex strength increases with: larger background temperature gradients; warmer background temperatures; larger initial temperature perturbations; higher Reynolds number; and higher resolution. In the simulations presented here vortices form when the background temperatures are $\sim 200K$ and vary radially as $r^{-0.25}$, the initial vorticity perturbations are zero, the initial temperature perturbations are 5% of the background, and the Reynolds number is $10^9$. A sensitivity study consisting of 74 simulations showed that as resolution and Reynolds number increase, vortices can form with smaller initial temperature perturbations, lower background temperatures, and smaller background temperature gradients. For the parameter ranges of these simulations, the disk is shown to be convectively stable by the Solberg-H{\o}iland criteria.

Fourier--Padé approximations and filtering for spectral simulations of an incompressible Boussinesq convection problem

In this paper, we present rational approximations based on Fourier series representation. For periodic piecewise analytic functions, the well-known Gibbs phenomenon hampers the convergence of the standard Fourier method. Here, for a given set of the Fourier coefficients from a periodic piecewise analytic function, we define Fourier-Pade-Galerkin and Fourier-Pade collocation methods by expressing the coefficients for the rational approximations using the Fourier data. We show that those methods converge exponentially in the smooth region and successfully reduce the Gibbs oscillations as the degrees of the denominators and the numerators of the Pade approximants increase. Numerical results are demonstrated in several examples. The collocation method is applied as a postprocessing step to the standard pseudospectral simulations for the one-dimensional inviscid Burgers' equation and the two-dimensional incompressible inviscid Boussinesq convection flow.

In situ acoustic scattering, attenuation, and dispersion measurements in marine sediments

The degree of acoustic velocity dispersion and nonlinear frequency dependency of attenuation in various types of sediments is investigated. A down‐looking chirp sonar and a wideband acoustic probe system are used to measure compressional wave speed and attenuation within 2–150‐kHz frequency band. Measurements indicate slight velocity dispersion and nonlinear frequency dependency of attenuation in silty and sandy sediments that can be effectively modeled by an extended Biot theory (Yamamoto and Turgut, J. Acoust. Soc. Am. 83, 1744–1751 (1988)]. The degree of acoustic scattering due to volume inhomogeneities and surface roughness is also estimated from the incoherent components of the measured signals. Acoustic scattering results are favorably compared with those predicted by pseudospectral FDTD simulations that used colocated vibracore data as input. Robustness of the attenuation and dispersion measurements by using chirp sonar and acoustic probe systems is also verified by the FDTD simulations [Work supported by the ONR.]

Pseudospectral time domain simulations of multiple light scattering in three‐dimensional macroscopic random media

We report a full‐vector, three‐dimensional, numerical solution of Maxwell's equations for optical propagation within, and scattering by, a random medium of macroscopic dimensions. The total scattering cross section is determined using the pseudospectral time domain technique. Specific results reported in this paper indicate that multiply scattered light also contains information that can be extracted by the proposed cross‐correlation analysis. On a broader perspective, our results demonstrate the feasibility of accurately determining the optical characteristics of arbitrary, macroscopic random media, including geometries with continuous variations of refractive index. Specifically, our results point toward the new possibilities of tissue optics; by numerically solving Maxwell's equations, the optical properties of tissue structures can be determined unambiguously.

Pseudo-spectral numerical simulation of miscible fluids with a high density ratio

A computational algorithm is described for direct numerical simulation (DNS) of turbulent mixing of two incompressible miscible fluids having greatly differing densities. The algorithm uses Fourier pseudo-spectral methods to compute spatial derivatives and a fractional step method involving the third-order Adams–Bashforth–Moulton predictor–corrector scheme to advance the solution in time. The pressure projection technique is shown to eliminate stability problems, previously observed, when the ratio of the densities in the two streams is as high as 35. The algorithm is investigated in detail for mixing in isotropic homogeneous turbulence of two fluids with a density ratio of 10. The limit on the density ratio is imposed so that the flow is both everywhere turbulent and spatially resolved. Both fluids have the same molecular viscosities, the nominal Schmidt number is 0.7, and the initial nominal Reynolds number based on the integral length scale and the rms velocity is 158. No body force is considered. It is shown that the pressure projection scheme does not limit the temporal accuracy of the solution when periodic boundary conditions are used, but that it significantly affects the stability of the simulations. It is also shown that the rate at which turbulence kinetic energy dissipates averaged for the whole computational domain is almost unaffected by density ratio.

The Effect of Coherent Structures on Stochastic Acceleration in MHD Turbulence

We investigate the influence of coherent structures on particle acceleration in the strongly turbulent solar corona. By randomizing the Fourier phases of a pseudospectral simulation of isotropic magnetohydrodynamic (MHD) turbulence(Re � 300)andtracingcollisionlesstestprotonsinboththeexact-MHDandphase-randomizedfields,itis foundthatthephasecorrelationsenhancetheaccelerationefficiencyduringthefirstadiabaticstageoftheacceleration process. The underlying physical mechanism is identified as the dynamical MHD alignment of the magnetic field with the electric current, which favors parallel (resistive) electric fields responsible for initial injection. Conversely, the alignment of the magnetic field with the bulk velocity weakens the acceleration by convective electric fields � u< b at a nonadiabatic stage of the acceleration process. We point out that nonphysical parallel electric fields in random-phase turbulence proxies lead to artificial acceleration and that the dynamical MHD alignment can be taken intoaccountonthelevelofthejointtwo-pointfunctionofthemagneticandelectricfieldsandisthereforeamenableto Fokker-Planck descriptions of stochastic acceleration. Subject headings: acceleration of particles — methods: numerical — MHD

Fourier spectral simulations and Gegenbauer reconstructions for electromagnetic waves in the presence of a metal nanoparticle

We describe Fourier pseudospectral time-domain simulations, carried out in order to study light interacting with a metallic nanoscale object. The difficulty of using Fourier methods to accurately predict the electromagnetic scattering in such problems arises from the discontinuity in the dielectric function along the surface of the metallic object. Standard Fourier methods lead to oscillatory behavior in approximating solutions that are nonsmooth or that have steep gradients. By applying the Gegenbauer reconstruction technique as a postprocessing method to the Fourier pseudospectral solution, we successfully reduce the oscillations after postprocessing. Our computational results, including comparison with finite-difference time-domain simulations, demonstrate the efficiency and accuracy of the method.

On the role of the smallest scales of a passive scalar field in a near-wall turbulent flow

Role of the smallest diffusive scales of a passive scalar field in the near-wall turbulent flow was examined with pseudo-spectral numerical simulations. Temperature fields were analyzed at friction Reynolds number Re τ=171 and at Prandtl numbers, Pr=1 and Pr=5.4. Results of direct numerical simulations (DNS) were compared with the under-resolved simulations where the velocity field was still resolved with the DNS accuracy, while a coarser grid was used to describe the temperature fields. Since the smallest temperature scales remained unresolved in these simulations, an appropriate spectral turbulent thermal diffusivity was applied to avoid pile-up at the higher wave numbers. In spite of coarser numerical grids, the temperature fields are still highly correlated with the DNS results, including instantaneous temperature fields. Results point to practically negligible role of the diffusive temperature scales on the macroscopic behavior of the turbulent heat transfer.

Numerical study of the decay of enstrophy in a two-dimensional Navier–Stokes fluid in the limit of very small viscosities

The rate of enstrophy decay in two-dimensional Navier–Stokes turbulence is explored numerically with high resolution pseudospectral simulations. Initial conditions and all other physical parameters are held fixed while viscosity is decreased through about three orders of magnitude, while increasing the numerical resolution as required to resolve the smallest scales that develop. The attempt is to discern whether or not a regime is being approached in which the behavior of the vorticity decay could be approximated as ideal and nondissipative. Although the enstrophy decay shows a (logarithmic) tendency to go to zero as the viscosity is decreased, at fixed times, the extrapolated values of viscosity to approach a zero enstrophy decay regime lie many orders of magnitude out of our currently available numerical resolution.

Large eddy simulation and study of the urban boundary layer

Based on a pseudo-spectral large eddy simulation (LES) model, an LES model with an anisotropy turbulent kinetic energy (TKE) closure model and an explicit multi-stage third-order Runge-Kutta scheme is established. The modeling and analysis show that the LES model can simulate the planetary boundary layer (PBL) with a uniform underlying surface under various stratifications very well. Then, similar to the description of a forest canopy, the drag term on momentum and the production term of TKE by subgrid city buildings are introduced into the LES equations to account for the area-averaged effect of the subgrid urban canopy elements and to simulate the meteorological fields of the urban boundary layer (UBL). Numerical experiments and comparison analysis show that: (1) the result from the LES of the UBL with a proposed formula for the drag coefficient is consistent and comparable with that from wind tunnel experiments and an urban subdomain scale model; (2) due to the effect of urban buildings, the wind velocity near the canopy is decreased, turbulence is intensified, TKE, variance, and momentum flux are increased, the momentum and heat flux at the top of the PBL are increased, and the development of the PBL is quickened; (3) the height of the roughness sublayer (RS) of the actual city buildings is the maximum building height (1.5–3 times the mean building height), and a constant flux layer (CFL) exists in the lower part of the UBL.