Perfectly matched layers for time-harmonic elastodynamics of unbounded domains: theory and finite-element implementation

One approach to the numerical solution of a wave equation on an unbounded domain uses a bounded domain surrounded by an absorbing boundary or layer that absorbs waves propagating outwards from the bounded domain. A perfectly matched layer (PML) is an unphysical absorbing layer model for linear wave equations that absorbs, almost perfectly, outgoing waves of all non-tangential angles-of-incidence and of all non-zero frequencies. This paper develops the PML concept for time-harmonic elastodynamics in Cartesian coordinates, utilising insights obtained with electromagnetics PMLs, and presents a novel displacement-based, symmetric finite-element implementation of the PML for time-harmonic plane-strain or three-dimensional motion. The PML concept is illustrated through the example of a one-dimensional rod on elastic foundation and through the anti-plane motion of a two-dimensional continuum. The concept is explored in detail through analytical and numerical results from a PML model of the semi-infinite rod on elastic foundation, and through numerical results for the anti-plane motion of a semi-infinite layer on a rigid base. Numerical results are presented for the classical soil–structure interaction problems of a rigid strip-footing on a (i) half-plane, (ii) layer on a half-plane, and (iii) layer on a rigid base. The analytical and numerical results obtained for these canonical problems demonstrate the high accuracy achievable by PML models even with small bounded domains.     

Long Time Behavior of the Perfectly Matched Layer Equations in Computational Electromagnetics

We investigate the long time behavior of two unsplit PML methods for the absorption of electromagnetic waves. Computations indicate that both methods suffer from a temporal instability after the fields reach a quiescent state. The analysis reveals that the source of the instability is the undifferentiated terms of the PML equations and that it is associated with a degeneracy of the quiescent systems of equations. This highlights why the instability occurs in special cases only and suggests a remedy to stabilize the PML by removing the degeneracy. Computational results confirm the stability of the modified equations and is used to address the efficacy of the modified schemes for absorbing waves.     

An Energy-Stable High-Order Central Difference Scheme for the Two-Dimensional Shallow Water Equations

An energy-stable high-order central finite difference scheme is derived for the two-dimensional shallow water equations. The scheme is mathematically formulated using the semi-discrete energy method for initial boundary value problems described in Olsson (1995, Math. Comput. 64, 1035–1065): after symmetrizing the equations via a change to entropy variables, the flux derivatives are entropy-split enabling the formulation of a semi-discrete energy estimate. We show experimentally that the entropy-splitting improves the stability properties of the fully discretized equations. Thus, the dependence on numerical dissipation to keep the scheme stable for long term time integrations is reduced relative to the original unsplit form, thereby decreasing non-physical damping of solutions. The numerical dissipation term used with the entropy-split equations is in a form which preserves the semi-discrete energy estimate. A random one-dimensional dam break calculation is performed showing that the shock speed is computed correctly for this particular case, however it is an open question whether the correct shock speed will be computed in generalMSC: 35Q35; 65M12; 65M06Supported in part by the New Zealand Marsden Fund, grant UOA827     

Mixed analytical/numerical method applied to the high Reynolds number k-ε turbulence model

In [Comput. Fluids 32 (2003) 659], a mixed analytical/numerical method for partial differential equation with an oscillating source term was proposed. The inhomogeneous partial differential equation is split into a homogeneous one plus an ODE for the source term using the time splitting method. The homogeneous part is then integrated numerically while the source term ODE is integrated analytically. This method was demonstrated to be efficient when the source term has a time scale much smaller than the mean flow time scale. In this paper, this approach is extended to the k-ε turbulence model for high Reynolds number flows. It is found that the mixed method based on operator splitting does not converge to a stable steady state solution. We thus propose an unsplit method. Numerical experiments for homogeneous turbulent flow and for a plane jet show that the unsplit mixed method substantially improves the accuracy of a time dependent problem while it has better convergence properties for a steady flow problem.     

Asymptotic Burnett model of the gas flow in a thin shock layer near a cylinder

Using the two-layer ideology of the theory of a viscous shock layer assuming the presence in the structure of the flow between the streamlined surface and the overrunning unperturbed stream of two (each with its own specific) characteristic areas-layers (smeared-out shock + shock layer itself), the complete Burnett equations are simplified in accordance with the problem of the cross-sectional flow-around of a round cylinder by the hypersound flow of rarified gas. The unitary asymptotic composite system of equations is formulated describing the flow in the entire thin shock layer including both of its above-mentioned structural areas (of the sublayer). Unlike the full Burnett equations, the obtained equation of the Burnett thin shock layer has an order not exceeding the order of the corresponding Navier-Stokes problem of the hypersonic thin shock layer, for which the research tool is sufficiently well developed.     

A 3D unsplit-advection volume tracking algorithm with planarity-preserving interface reconstruction

A new volume tracking method is introduced for tracking interfaces in three-dimensional (3D) geometries partitioned with orthogonal hexahedra. The method approximates interface geometries as piecewise planar, and advects volumes in a single unsplit step using fully multidimensional fluxes that have their definition based in backward-trajectory remapping. By using multidimensional unsplit advection, the expense of high-order interface reconstruction is incurred only once per timestep. Simple departures from strict backward-trajectory remapping remove any need for consideration of volume computations involving shapes consisting of non-planar ruled surfaces. Second-order accuracy of the method is demonstrated even for vigorous 3D deformations.     

Constrained multibody system dynamics an automated approach

The governing equations for constrained multibody systems are formulated in a manner suitable for their automated, numerical development and solution. Specifically, the “closed loop” problem of multibody chain systems is addressed.

The governing equations are developed by modifying dynamical equations obtained from Lagrange's form of d'Alembert's principle. This modification, which is based upon a solution of the constraint equations obtained through a “zero eigenvalues theorem,” is, in effect, a contraction of the dynamical equations.

It is observed that, for a system with n generalized coordinates and m constraint equations, the coefficients in the constraint equations may be viewed as “constraint vectors” in n-dimensional space. Then, in this setting the system itself is free to move in the n−m directions which are “orthogonal” to the constraint vectors.     

Magnetohydrodynamic three-dimensional flow of nanofluids with slip and thermal radiation over a nonlinear stretching sheet: a numerical study

A numerical simulation for mixed convective three-dimensional slip flow of water-based nanofluids with temperature jump boundary condition is presented. The flow is caused by nonlinear stretching surface. Conservation of energy equation involves the radiation heat flux term. Applied transverse magnetic effect of variable kind is also incorporated. Suitable nonlinear similarity transformations are used to reduce the governing equations into a set of self-similar equations. The subsequent equations are solved numerically by using shooting method. The solutions for the velocity and temperature distributions are computed for several values of flow pertinent parameters. Further, the numerical values for skin-friction coefficients and Nusselt number in respect of different nanoparticles are tabulated. A comparison between our numerical and already existing results has also been made. It is found that the velocity and thermal slip boundary condition showed a significant effect on momentum and thermal boundary layer thickness at the wall. The presence of nanoparticles stabilizes the thermal boundary layer growth.

     

Dynamical study of metallic clusters using the statistical method of time series clustering

We perform common neighbor analysis on the long-time series data generated by isothermal Brownian-type molecular dynamics simulations to study the thermal and dynamical properties of metallic clusters. In our common neighbor analysis, we introduce the common neighbor label (CNL) which is a group of atoms of a smaller size (than the cluster) designated by four numeric digits. The CNL thus describes topologically smaller size atomic configurations and is associated an abundance value which is the number of “degenerate” four digits all of which characterize the same CNL. When the cluster is in its lowest energy state, it has a fixed number of CNLs and hence abundances. At nonzero temperatures, the cluster undergoes different kinds of atomic activities such as vibrations, migrational relocation, permutational and topological isomer transitions, etc. depending on its lowest energy structure. As a result, the abundances of CNLs at zero temperature will change and new CNLs with their respective new abundances are created. To understand the temperature dependence of the CNL dynamics, and hence shed light on the cluster dynamics itself, we employ a novel method of statistical time series analysis. In this method, we perform statistical clustering at two time scales. First, we examine, at given temperature, the signs of abundance changes at a short-time scale, and assign CNLs to two short-time clusters. Quasi-periodic features can be seen in the time evolution of these short-time clusters, based on which we choose a long-time scale to compute the long-time correlations between CNL pairs. We then exploit the separation of correlation levels seen in these long-time correlations to extract strongly-correlated collections of CNLs, which we will identify as effective variables for the long-time cluster dynamics. It is found that certain effective variables show subtleties in their temperature dependences and these thermal traits bear a delicate relation to prepeaks and main peaks seen in clusters Ag₁₄, Cu₁₄ and Cu₁₃Au₁. We therefore infer from the temperature changes of effective variables and locate the temperatures at which these prepeaks and principal peaks appear, and they are evaluated by comparing with those deduced from the specific heat data.     

Unified framework for numerical methods to solve the time-dependent Maxwell equations

We present a comparative study of numerical algorithms to solve the time-dependent Maxwell equations for systems with spatially varying permittivity and permeability. We show that the Lie-Trotter-Suzuki product-formula approach can be used to construct a family of unconditionally stable algorithms, the conventional Yee algorithm, and two new variants of the Yee algorithm that do not require the use of the staggered-in-time grid. We also consider a one-step algorithm, based on the Chebyshev polynomial expansion, and compare the computational efficiency of the one-step, the Yee-type, the alternating-direction-implicit, and the unconditionally stable algorithms. For applications where the long-time behavior is of main interest, we find that the one-step algorithm may be orders of magnitude more efficient than present multiple time-step, finite-difference time-domain algorithms.     

Symplectic and multi-symplectic methods for coupled nonlinear Schrödinger equations with periodic solutions

We consider for the integration of coupled nonlinear Schrödinger equations with periodic plane wave solutions a splitting method from the class of symplectic integrators and the multi-symplectic six-point scheme which is equivalent to the Preissman scheme. The numerical experiments show that both methods preserve very well the mass, energy and momentum in long-time evolution. The local errors in the energy are computed according to the discretizations in time and space for both methods. Due to its local nature, the multi-symplectic six-point scheme preserves the local invariants more accurately than the symplectic splitting method, but the global errors for conservation laws are almost the same.     

Influence of centrifugal forces on the development of wave packets in boundary layers with a uniformly valid model

The behavior of three-dimensional wave packets in the boundary layer on curved surfaces is analyzed in this study based on a modification of the triple-deck theory referred to as the “criss-cross” interaction model. The equations of the criss-cross interaction describe a particular type of boundary layer instability mode that possesses underlying properties of both the Tollmien-Schlichting waves and Taylor-Görtler vortices. Previous analysis of the criss-cross interaction regime suggests a possibility for upstream propagation of perturbations in the boundary layer and possible absolute instability of the flow. However, these results cannot be considered as conclusive because the initial-value problem for the criss-cross interaction equations is ill-posed. In a recent work [Turkyilmazoglu M, Ruban AI. A uniformly valid well-posed asymptotic approach to the inviscid-viscous interaction theory. Stud Appl Math 2002;108:161-85] a regularized non-asymptotic model to describe criss-cross interaction has been proposed. Whereas in the original version of the theory, perturbations have an unbounded growth rate as the longitudinal wave number ∣k∣ → ∞, in the new model of [Turkyilmazoglu M, Ruban AI. A uniformly valid well-posed asymptotic approach to the inviscid-viscous interaction theory. Stud Appl Math 2002;108:161-85], as physically expected the amplification rate remains bounded for both spatially growing and temporally growing waves. A Fourier transform method is used in the present study to solve the linearized equations for the flow over concave roughness and humps and it is found that disturbances develop and are convected downstream as wave packets. The behavior of the wave packets is consistent with convective instability, and the upstream influence is no longer present at large times.     

Chaotic vibration of the one-dimensional linear wave equation with a van der Pol nonlinear boundary condition

The one-dimensional linear wave equation with a van der Pol nonlinear boundary condition is one of the simplest models that may cause isotropic or nonisotropic chaotic vibrations. It characterizes the nonisotropic chaotic vibration by means of the total variation theory. Some results are derived on the exponential growth of total variation of the snapshots on the spatial interval in the long-time horizon when the map and the initial condition satisfy some conditions.     

On the immersed interface method for solving time-domain Maxwell's equations in materials with curved dielectric interfaces

This paper deals with accurate numerical simulation of two-dimensional time-domain Maxwell's equations in materials with curved dielectric interfaces. The proposed fully second-order scheme is a hybridization between the immersed interface method (IIM), introduced to take into account curved geometries in structured schemes, and the Lax-Wendroff scheme, usually used to improve order of approximations in time for partial differential equations. In particular, the IIM proposed for two-dimensional acoustic wave equations with piecewise constant coefficients [C. Zhang, R.J. LeVeque, The immersed interface method for acoustic wave equations with discontinuous coefficients, Wave Motion 25 (1997) 237-263] is extended through a simple least squares procedure to such Maxwell's equations. Numerical results from the simulation of electromagnetic scattering of a plane incident wave by a dielectric circular cylinder appear to indicate that, compared to the original IIM for the acoustic wave equations, the augmented IIM with the proposed least squares fitting greatly improves the long-time stability of the time-domain solution. Semi-discrete finite difference schemes using the IIM for spatial discretization are also discussed and numerically tested in the paper.     

A Chebyshev pseudospectral multidomain method for the soliton solution of coupled nonlinear Schrödinger equations

The collision of solitary waves is an important problem in both physics and applied mathematics. In this paper, we study the solution of coupled nonlinear Schrödinger equations based on pseudospectral collocation method with domain decomposition algorithm for approximating the spatial variable. The problem is converted to a system of nonlinear ordinary differential equations which will be integrated in time by explicit Runge–Kutta method of order four. The multidomain scheme has much better stability properties than the single domain. Thus this permits using much larger step size for the time integration which fulfills stability restrictions. The proposed scheme reduces the effects of round-of-error for the Chebyshev collocation and also uses less memory without sacrificing the accuracy. The numerical experiments are presented which show the multidomain pseudospectral method has excellent long-time numerical behavior and preserves energy conservation property.     

Numerical simulation of growth of a vapor bubble during flow boiling of water in a microchannel

The present study is performed to numerically analyze the growth of a vapor bubble during flow of water through a microchannel. The complete Navier–Stokes equations, along with continuity and energy equations, are solved using the SIMPLER (semi-implicit method for pressure-linked equations revised) method. The liquid–vapor interface is captured using the level set technique. The microchannel is 200-m square in cross-section and the bubble is placed at the center of the channel with superheated liquid around it. The results show steady initial bubble growth followed by a rapid axial expansion after the bubble fills the channel cross-section. A trapped liquid layer is observed between the bubble and the channel as it elongates. The bubble growth rate increased with the incoming liquid superheat, but decreased with Reynolds number. The formation of a vapor patch at the walls is found to be dependent on the time the bubble takes to fill up the channel. The upstream interface of the bubble is found to exhibit both forward and reverse movement during bubble growth. The results show little effect of gravity on the bubble growth rate under the specified conditions. The bubble growth features obtained from the numerical results are found to be qualitatively similar to experimental observations.     

Numerical study of a nonlinear heat equation for plasma physics

This paper is devoted to the numerical approximation of a nonlinear parabolic balance equation, which describes the heat evolution of a magnetically confined plasma in the edge region of a tokamak. The nonlinearity implies some numerical difficulties, in particular for the long-time behaviour approximation, when solved with standard methods. An efficient numerical scheme is presented in this paper, based on a combination of a directional splitting scheme and the implicit–explicit scheme introduced in Filbet and Jin [A class of asymptotic preserving schemes for kinetic equations and related problems with stiff sources, J. Comput. Phys. 229 (2010), pp. 7625–7648].     

Solution of a hyperbolic system of turbulence-model equations by the “box” scheme

The momentum and continuity equations for a two-dimensional boundary layer, together with an empirical first-order partial differential equation for shear-stress transport, form a hyperbolic system of equations for u, v and τ. Here the Bradshaw-Ferriss-Atwell version of this turbulence model is solved by the Keller-Cebeci “box” scheme, which is particularly suited to systems of equations that are individually of first order. Computing time is about equal to that taken by a method-of-characteristics program if the same number of grid points are used across the layer and in the streamwise direction. However, the box scheme allows larger x-steps to be taken in the streamwise direction leading to smaller computing times.     

A multi-director formulation for nonlinear elastic-viscoelastic layered shells

A C⁰ finite element formulation for nonlinear analysis of multi-layered shells comprised of elastic and viscoelastic layers is presented for applications involving small strains but finite rotations. The elastic and viscoelastic layers may occupy arbitrary layer locations and the formulation is applicable to thick and thin shells. The formulation utilizes a three-dimensional variational approach in which the layered shell is represented as a multi-director field. The incorporated kinematic theory describes, within individual layers, the effects of transverse shear and transverse normal strain to arbitrary orders in the layer thickness coordinate. Stresses are computed through the three-dimensional constitutive equations and the usual “zero normal stress” shell hypothesis is not employed. Sufficiently general constitutive equations for the viscoelastic layers are proposed in objective rate form and a product algorithm, based on an operator split in the complete set of constitutive equations, is used for the temporal integration of the rate equations. The definition of the tangent operator, used in Newton's method for the solution of the nonlinear equations, is derived consistently from the product algorithm. Observations on the use of reduced/selective integration in the presence of high order kinematics are made and a number of numerical examples are presented to illustrate the capability of the formulation.     

Normal mode analysis for resistive cylindrical plasmas

The compressible, resistive MHD equations are linearized around an equilibrium with cylindrical symmetry and solved numerically as a complex eigenvalue problem. This normal mode code allows one to solve for very small resistivity γ≈10^-10. The scaling of growth rates and layer width agrees very well with analytical theory. Especially, the influence of both current and pressure on the instabilities is studied in detail; the effect of resistivity on the ideally unstable internal kink is analyzed.