The Three Dimensional Non-conforming Finite Element Solution of the Chapman–Ferraro Problem

We demonstrate the feasibility of using a non-conforming, piecewise harmonic finite element method on an unstructured grid in solving a magnetospheric physics problem. We use this approach to construct a global discrete model of the magnetic field of the magnetosphere that includes the effects of shielding currents at the outer boundary (the magnetopause). As in the approach of F. R. Toffolettoet al.(1994,Geophys. Res. Lett.21, 7) the internal magnetospheric field model is that of R. V. Hilmer and G.-H. Voigt (1995,J. Geophys. Res.) while the magnetopause shape is based on an empirically determined approximation (1997, J. Shueet al.,J. Geophys. Res.102, 9497). The results is a magnetic field model whose field lines are completely confined within the magnetosphere. The presented numerical results indicate that the discrete non-conforming finite element model is well-suited for magnetospheric field modeling.     

A Finite Volume Method for the Approximation of Diffusion Operators on Distorted Meshes

A new finite volume method is presented for discretizing general linear or nonlinear elliptic second-order partial-differential equations with mixed boundary conditions. The advantage of this method is that arbitrary distorted meshes can be used without the numerical results being altered. The resulting algorithm has more unknowns than standard methods like finite difference or finite element methods. However, the matrices that need to be inverted are positive definite, so the most powerful linear solvers can be applied. The method has been tested on a few elliptic and parabolic equations, either linear, as in the case of the standard heat diffusion equation, or nonlinear, as in the case of the radiation diffusion equation and the resistive diffusion equation with Hall term.     

A High-Order Discontinuous Galerkin Method for 2D Incompressible Flows

In this paper we introduce a high-order discontinuous Galerkin method for two-dimensional incompressible flow in the vorticity stream-function formulation. The momentum equation is treated explicitly, utilizing the efficiency of the discontinuous Galerkin method. The stream function is obtained by a standard Poisson solver using continuous finite elements. There is a natural matching between these two finite element spaces, since the normal component of the velocity field is continuous across element boundaries. This allows for a correct upwinding gluing in the discontinuous Galerkin framework, while still maintaining total energy conservation with no numerical dissipation and total enstrophy stability. The method is efficient for inviscid or high Reynolds number flows. Optimal error estimates are proved and verified by numerical experiments.     

Physics-Based Preconditioning and the Newton–Krylov Method for Non-equilibrium Radiation Diffusion

An algorithm is presented for the solution of the time dependent reaction-diffusion systems which arise in non-equilibrium radiation diffusion applications. This system of nonlinear equations is solved by coupling three numerical methods, Jacobian-free Newton–Krylov, operator splitting, and multigrid linear solvers. An inexact Newton's method is used to solve the system of nonlinear equations. Since building the Jacobian matrix for problems of interest can be challenging, we employ a Jacobian–free implementation of Newton's method, where the action of the Jacobian matrix on a vector is approximated by a first order Taylor series expansion. Preconditioned generalized minimal residual (PGMRES) is the Krylov method used to solve the linear systems that come from the iterations of Newton's method. The preconditioner in this solution method is constructed using a physics-based divide and conquer approach, often referred to as operator splitting. This solution procedure inverts the scalar elliptic systems that make up the preconditioner using simple multigrid methods. The preconditioner also addresses the strong coupling between equations with local 2×2 block solves. The intra-cell coupling is applied after the inter-cell coupling has already been addressed by the elliptic solves. Results are presented using this solution procedure that demonstrate its efficiency while incurring minimal memory requirements.     

Using the seventh-order numerical method to solve first-order nonlinear coupled-wave equations for degenerate two-wave and four-wave mixing

Using a new seventh-order numerical method [theO(h ⁷) method] for solving two-point boundary value problems, numerical solutions of the first-order nonlinear coupledwave equations for degenerate two-wave and four-wave mixing in a reflection geometry have been obtained. A computer program employing the Gauss-Jordan elimination technique has also been adopted to effectively solve the resultant large, sparse and unsymmetric matrix, obtained from theO(h ⁷) method and the Newton-Raphson iteration method. Numerical results from the computer calculations are presented graphically. A comparison between thisO(h ⁷) method and the shooting method, mainly from the viewpoint of computational efficiency, is also made.     

The Adjoint Method for an Inverse Design Problem in the Directional Solidification of Binary Alloys

This paper presents a systematic procedure based on the adjoint method for solving a class of inverse directional alloy solidification design problems in which a desired growth velocityv_fis achieved under stable growth conditions. To the best of our knowledge, this is the first time that a continuum adjoint formulation is proposed for the solution of an inverse problem with simultaneous heat and mass transfer, thermo-solutal convection, and phase change. In this paper, the interfacial stability is considered to imply a sharp solid–liquid freezing interface. This condition is enforced using the constitutional undercooling criterion in the form of an inequality constraint between the thermal and solute concentration gradients,GandG_c, respectively, at the freezing front. The main unknowns of the design problem are the heating and/or cooling boundary conditions on the mold walls. The inverse design problem is formulated as a functional optimization problem. The cost functional is defined by the square of theL₂norm of the deviation of the freezing interface temperature from the temperature corresponding to thermodynamic equilibrium. A continuum adjoint system is derived to calculate the adjoint temperature, concentration, and velocity fields such that the gradient of the cost functional can be expressed analytically. The cost functional minimization process is realized by the conjugate gradient method via the finite element method solutions of the continuum direct, sensitivity, and adjoint problems. The developed formulation is demonstrated with an example of designing the directional solidification of a binary aqueous solution in a rectangular mold such that a stable vertical interface advances from left to right with a desired growth velocity.     

Approximation and Reconstruction of the Electrostatic Field in Wire–Plate Precipitators by a Low-Order Model

The numerical computation of the ionic space charge and electric field produced by corona discharge in a wire–plate electrostatic precipitator (ESP) is considered. The electrostatic problem is defined by a reduced set of the Maxwell equations. Since self-consistent conditions at the wire and at the plate cannot be specified a priori, a time-consuming iterative numerical procedure is required. The efficiency of all numerical solvers of the reduced Maxwell equations depends in particular on the accuracy of the initial guess solution. The objectives of this work are two: first, we propose a semianalytical technique based on the Karhunen–Loève (KL) decomposition of the current density field J, which can significantly improve the performance of a numerical solver; second, we devise a procedure to reconstruct the complete electric field from a given J. The approximate solution of the current density field is based on the derivation of an analytical approximation , which, added to a linear combination of few KL basis functions, constitutes an accurate approximation of J. In the first place, this result is useful for optimization procedures of the current density field, which involve the computation of many different configurations. Second, we show that from the current density field we can obtain an accurate estimate for the complete electrostatic field which can be used to speed up the convergence of the iterative procedure of standard numerical solvers.     

非线性波动方程的新数值迭代方法

提出了一种新的求解非线性波动方程的数值迭代法,它是一种半解析的方法.与完全的数值计算方法扰法相比,它能够考虑各阶谐波的相互作用,且能够满足能量守恒定律.用它研究了非线性声波在液体中的传播性质,结果表明,在微扰法适用的声强范围内迭代法也适用,在微扰法不适用的一个较宽的声强范围内迭代法依然适用.     

Convection-Diffusion Lattice Boltzmann Scheme for Irregular Lattices

In this paper, a lattice Boltzmann (LB) scheme for convection diffusion on irregular lattices is presented, which is free of any interpolation or coarse graining step. The scheme is derived using the axioma that the velocity moments of the equilibrium distribution equal those of the Maxwell–Boltzmann distribution. The axioma holds for both Bravais and irregular lattices, implying a single framework for LB schemes for all lattice types. By solving benchmark problems we have shown that the scheme is indeed consistent with convection diffusion. Furthermore, we have compared the performance of the LB schemes with that of finite difference and finite element schemes. The comparison shows that the LB scheme has a similar performance as the one-step second-order Lax–Wendroff scheme: it has little numerical diffusion, but has a slight dispersion error. By changing the relaxation parameter ω the dispersion error can be balanced by a small increase of the numerical diffusion.     

A nonconforming combination for solving elliptic problems with interfaces

Nonconforming combinations are provided for solving interface problems of elliptic equations. In these approaches, the Ritz-Galerkin method with particular solutions is used for the part of a solution domain where there are interface singular points; and the conventional finite element method is used for the rest of the solution domain. In addition, admissible functions chosen are constrained to be continuous only at the element nodes on the common boundary of the subdomains. Error bounds are derived in the Sobolev norms, and numerical experiments are given for solving a model interface problem of the equation, −Δu + U = 0. Moreover, a significant coupling relation, L + 1 = O(|ln h|), is found for interface problems by using the nonconforming combinations, where (L + 1) is the total number of particular solutions used in the Ritz-Galerkin method, and h is the maximal boundary length of triangular elements in the finite element method.     

Unconditionally Stable Methods for Hamilton–Jacobi Equations

We present new numerical methods for constructing approximate solutions to the Cauchy problem for Hamilton–Jacobi equations of the form u_t+H(D_xu)=0. The methods are based on dimensional splitting and front tracking for solving the associated (non-strictly hyperbolic) system of conservation laws p_t+D_xH(p)=0, where p=D_xu. In particular, our methods depend heavily on a front tracking method for one-dimensional scalar conservation laws with discontinuous coefficients. The proposed methods are unconditionally stable in the sense that the time step is not limited by the space discretization and they can be viewed as “large-time-step” Godunov-type (or front tracking) methods. We present several numerical examples illustrating the main features of the proposed methods. We also compare our methods with several methods from the literature.     

Implementation Aspects of 3D Lattice-BGK: Boundaries, Accuracy, and a New Fast Relaxation Method

In many realistic fluid-dynamical simulations the specification of the boundary conditions, the error sources, and the number of time steps to reach a steady state are important practical considerations. In this paper we study these issues in the case of the lattice-BGK model. The objective is to present a comprehensive overview of some pitfalls and shortcomings of the lattice-BGK method and to introduce some new ideas useful in practical simulations. We begin with an evaluation of the widely used bounce-back boundary condition in staircase geometries by simulating flow in an inclined tube. It is shown that the bounce-back scheme is first-order accurate in space when the location of the non-slip wall is assumed to be at the boundary nodes. Moreover, for a specific inclination angle of 45 degrees, the scheme is found to be second-order accurate when the location of the non-slip velocity is fitted halfway between the last fluid nodes and the first solid nodes. The error as a function of the relaxation parameter is in that case qualitatively similar to that of flat walls. Next, a comparison of simulations of fluid flow by means of pressure boundaries and by means of body force is presented. A good agreement between these two boundary conditions has been found in the creeping-flow regime. For higher Reynolds numbers differences have been found that are probably caused by problems associated with the pressure boundaries. Furthermore, two widely used 3D models, namelyD₃Q₁₅andD₃Q₁₉, are analysed. It is shown that theD₃Q₁₅model may induce artificial checkerboard invariants due to the connectivity of the lattice. Finally, a new iterative method, which significantly reduces the saturation time, is presented and validated on different benchmark problems.     

A moving grid finite element method applied to a model biological pattern generator

Many problems in biology involve growth. In numerical simulations it can therefore be very convenient to employ a moving computational grid on a continuously deforming domain. In this paper we present a novel application of the moving grid finite element method to compute solutions of reaction–diffusion systems in two-dimensional continuously deforming Euclidean domains. A numerical software package has been developed as a result of this research that is capable of solving generalised Turing models for morphogenesis.     

球几何中子输运保正线性间断有限元格式

针对球几何中子输运方程线性间断有限元方法计算的负中子通量问题,构造了保正线性间断有限元格式,该格式保持中子角通量0阶矩和1阶矩。现有方法计算中子角通量非负时,采用传统的线性间断有限元方法,求解线性方程组;原方法计算出现负通量,则采用构造的保正格式,求解非线性方程组。编制了球几何中子输运问题保正格式程序模块,并集成到应用程序。数值算例表明构造的保正格式计算的中子通量非负,有效降低数值误差,提高数值计算的精度。     

Solution of Time-Dependent Diffusion Equations with Variable Coefficients Using Multiwavelets

A new numerical algorithm is developed for the solution of time-dependent differential equations of diffusion type. It allows for an accurate and efficient treatment of multidimensional problems with variable coefficients, nonlinearities, and general boundary conditions. For space discretization we use the multiwavelet bases introduced by Alpert (1993,SIAM J. Math. Anal.24, 246–262), and then applied to the representation of differential operators and functions of operators presented by Alpert, Beylkin, and Vozovoi (Representation of operators in the multiwavelet basis, in preparation). An important advantage of multiwavelet basis functions is the fact that they are supported only on non-overlapping subdomains. Thus multiwavelet bases are attractive for solving problems in finite (non periodic) domains. Boundary conditions are imposed with a penalty technique of Hesthaven and Gottlieb (1996,SIAM J. Sci. Comput., 579–612) which can be used to impose rather general boundary conditions. The penalty approach was extended to a procedure for ensuring the continuity of the solution and its first derivative across interior boundaries between neighboring subdomains while time stepping the solution of a time dependent problem. This penalty procedure on the interfaces allows for a simplification and sparsification of the representation of differential operators by discarding the elements responsible for interactions between neighboring subdomains. Consequently the matrices representing the differential operators (on the finest scale) have block-diagonal structure. For a fixed order of multiwavelets (i.e., a fixed number of vanishing moments) the computational complexity of the present algorithm is proportional to the number of subdomains. The time discretization method of Beylkin, Keiser, and Vozovoi (1998, PAM Report 347) is used in view of its favorable stability properties. Numerical results are presented for evolution equations with variable coefficients in one and two dimensions.     

An energy law preserving C finite element scheme for simulating the kinematic effects in liquid crystal dynamics

In this paper, we use finite element methods to simulate the hydrodynamical systems governing the motions of nematic liquid crystals in a bounded domain Ω. We reformulate the original model in the weak form which is consistent with the continuous dissipative energy law for the flow and director fields in W^1,2+σ(Ω) (σ > 0 is an arbitrarily small number). This enables us to use convenient conformal C⁰ finite elements in solving the problem. Moreover, a discrete energy law is derived for a modified midpoint time discretization scheme. A fixed iterative method is used to solve the resulted nonlinear system so that a matrix free time evolution may be achieved and velocity and director variables may be solved separately. A number of hydrodynamical liquid crystal examples are computed to demonstrate the effects of the parameters and the performance of the method.     

An Iterative Two-Grid Method of a Finite Element PML Approximation for the Two Dimensional Maxwell Problem

In this paper, we propose an iterative two-grid method for the edge finite element discretizations (a saddle-point system) of Perfectly Matched Layer (PML) equations to the Maxwell scattering problem in two dimensions. Firstly, we use a fine space to solve a discrete saddle-point system of $H(grad)$ variational problems, denoted by auxiliary system 1. Secondly, we use a coarse space to solve the original saddle-point system. Then, we use a fine space again to solve a discrete$\boldsymbol{H}(curl)$-elliptic variational problems, denoted by auxiliary system 2. Furthermore, we develop a regularization diagonal block preconditioner for auxiliary system 1 and use $H$-$X$ preconditioner for auxiliary system 2. Hence we essentially transform the original problem in a fine space to a corresponding (but much smaller) problem on a coarse space, due to the fact that the above two preconditioners are efficient and stable. Compared with some existing iterative methods for solving saddle-point systems, such as PMinres, numerical experiments show the competitive performance of our iterative two-grid method.     

Numerical comparisons of time–space iterative method and spatial iterative methods for the stationary Navier–Stokes equations

This paper makes some numerical comparisons of time–space iterative method and spatial iterative methods for solving the stationary Navier–Stokes equations. The time–space iterative method consists in solving the nonstationary Stokes equations based on the time–space discretization by the Euler implicit/explicit scheme under a weak uniqueness condition (A2). The spatial iterative methods consist in solving the stationary Stokes scheme, Newton scheme, Oseen scheme based on the spatial discretization under some strong uniqueness assumptions. We compare the stability and convergence conditions of the time–space iterative method and the spatial iterative methods. Moreover, the numerical tests show that the time–space iterative method is the more simple than the spatial iterative methods for solving the stationary Navier–Stokes problem. Furthermore, the time–space iterative method can solve the stationary Navier–Stokes equations with some small viscosity and the spatial iterative methods can only solve the stationary Navier–Stokes equations with some large viscosities.     

Moving Mesh Methods in Multiple Dimensions Based on Harmonic Maps

In practice, there are three types of adaptive methods using the finite element approach, namely the h-method, p-method, and r-method. In the h-method, the overall method contains two parts, a solution algorithm and a mesh selection algorithm. These two parts are independent of each other in the sense that the change of the PDEs will affect the first part only. However, in some of the existing versions of the r-method (also known as the moving mesh method), these two parts are strongly associated with each other and as a result any change of the PDEs will result in the rewriting of the whole code. In this work, we will propose a moving mesh method which also contains two parts, a solution algorithm and a mesh-redistribution algorithm. Our efforts are to keep the advantages of the r-method (e.g., keep the number of nodes unchanged) and of the h-method (e.g., the two parts in the code are independent). A framework for adaptive meshes based on the Hamilton–Schoen–Yau theory was proposed by Dvinsky. In this work, we will extend Dvinsky's method to provide an efficient solver for the mesh-redistribution algorithm. The key idea is to construct the harmonic map between the physical space and a parameter space by an iteration procedure. Each iteration step is to move the mesh closer to the harmonic map. This procedure is simple and easy to program and also enables us to keep the map harmonic even after long times of numerical integration. The numerical schemes are applied to a number of test problems in two dimensions. It is observed that the mesh-redistribution strategy based on the harmonic maps adapts the mesh extremely well to the solution without producing skew elements for multi-dimensional computations.     

Accurate ω-ψ Spectral Solution of the Singular Driven Cavity Problem

This article provides accurate spectral solutions of the driven cavity problem, calculated in the vorticity–stream function representation without smoothing the corner singularities—a prima facie impossible task. As in a recent benchmark spectral calculation by primitive variables of Botella and Peyret, closed-form contributions of the singular solution for both zero and finite Reynolds numbers are subtracted from the unknown of the problem tackled here numerically in biharmonic form. The method employed is based on a split approach to the vorticity and stream function equations, a Galerkin–Legendre approximation of the problem for the perturbation, and an evaluation of the nonlinear terms by Gauss–Legendre numerical integration. Results computed for Re=0, 100, and 1000 compare well with the benchmark steady solutions provided by the aforementioned collocation–Chebyshev projection method. The validity of the proposed singularity subtraction scheme for computing time-dependent solutions is also established.