This paper presents a new approach to verify the accuracy of computational simulations. We develop mathematical theorems which can serve as robust a posteriori error estimation techniques to identify numerical pollution, check the performance of adaptive meshes, and verify numerical solutions. We demonstrate performance of this methodology on problems from flow thorough porous media. However, one can extend it to other models. We construct mathematical properties such that the solutions to Darcy and Darcy-Brinkman equations satisfy them. The mathematical properties include the total minimum mechanical power, minimum dissipation theorem, reciprocal relation, and maximum principle for the vorticity. All the developed theorems have firm mechanical bases and are independent of numerical methods. So, these can be utilized for solution verification of finite element, finite volume, finite difference, lattice Boltzmann methods and so forth. In particular, we show that, for a given set of boundary conditions, Darcy velocity has the minimum total mechanical power of all the kinematically admissible vector fields. We also show that a similar result holds for Darcy-Brinkman velocity. We then show for a conservative body force, the Darcy and Darcy-Brinkman velocities have the minimum total dissipation among their respective kinematically admissible vector fields. Using numerical examples, we show that the minimum dissipation and total mechanical power theorems can be utilized to identify pollution errors in numerical solutions. The solutions to Darcy and Darcy-Brinkman equations are shown to satisfy a reciprocal relation, which has the potential to identify errors in the numerical implementation of boundary conditions. It is also shown that the vorticity under both steady and transient Darcy-Brinkman equations satisfy maximum principles if the body force is conservative and the permeability is homogeneous and isotropic. A discussion on the nature of vorticity under steady and transient Darcy equations is also presented. Using several numerical examples, we will demonstrate the predictive capabilities of the proposed a posteriori techniques in assessing the accuracy of numerical solutions for a general class of problems, which could involve complex domains and general computational grids.

The level set method is one of the most successful methods for the simulation of multi-phase flows. To keep the level set function close the signed distance function, the level set function is constantly reinitialized by solving a Hamilton-Jacobi type of equation during the simulation. When the fluid interface intersects with a solid wall, a moving contact line forms and the reinitialization of the level set function requires a boundary condition in certain regions on the wall. In this work, we propose to use the dynamic contact angle, which is extended from the contact line, as the boundary condition for the reinitialization of the level set function. The reinitialization equation and the equation for the normal extension of the dynamic contact angle form a coupled system and are solved simultaneously. The extension equation is solved on the wall and it provides the boundary condition for the reinitialization equation; the level set function provides the directions along which the contact angle is extended from the contact line. The coupled system is solved using the 3rd order TVD Runge-Kutta method and the Godunov scheme. The Godunov scheme automatically identifies the regions where the angle condition needs to be imposed. The numerical method is illustrated by examples in three dimensions.

The statistical error associated with sampling in the DSMC method can be categorized as type I and II, which are caused by the incorrect rejection and acceptance of the null hypothesis, respectively. In this study, robust global and local automatic steady state detection methods were developed based on an ingenious method based purely on the statistics and kinetics of particles. The key concept is built upon probabilistic automatic reset sampling (PARS) to minimize the type II error caused by incorrect acceptance of the samples that do not belong to the steady state. The global steady state method is based on a relative standard variation of collisional invariants, while the local steady state method is based on local variations in the distribution function of particles at each cell. In order to verify the capability of the new methods, two benchmark cases — the one-dimensional shear-driven Couette flow and the two-dimensional high speed flow past a vertical wall—were extensively investigated. Owing to the combined effects of the automatic detection and local reset sampling, the local steady state detection method yielded a substantial gain of 30-36% in computational cost for the problem studied. Moreover, the local reset feature outperformed the automatic detection feature in overall computational savings.

In this paper, a diffuse-interface immersed boundary method (IBM) is proposed to treat three different thermal boundary conditions (Dirichlet, Neumann, Robin) in thermal flow problems. The novel IBM is implemented combining with the lattice Boltzmann method (LBM). The present algorithm enforces the three types of thermal boundary conditions at the boundary points. Concretely speaking, the IBM for the Dirichlet boundary condition is implemented using an iterative method, and its main feature is to accurately satisfy the given temperature on the boundary. The Neumann and Robin boundary conditions are implemented in IBM by distributing the jump of the heat flux on the boundary to surrounding Eulerian points, and the jump is obtained by applying the jump interface conditions in the normal and tangential directions. A simple analysis of the computational accuracy of IBM is developed. The analysis indicates that the Taylor-Green vortices problem which was used in many previous studies is not an appropriate accuracy test example. The capacity of the present thermal immersed boundary method is validated using four numerical experiments: (1) Natural convection in a cavity with a circular cylinder in the center; (2) Flows over a heated cylinder; (3) Natural convection in a concentric horizontal cylindrical annulus; (4) Sedimentation of a single isothermal cold particle in a vertical channel. The numerical results show good agreements with the data in the previous literatures.

It is well-known that the traditional full integral quadrilateral element fails to provide accurate results to the Helmholtz equation with large wave numbers due to the “pollution error” caused by the numerical dispersion. To overcome this deficiency, this paper proposed an element decomposition method (EDM) for analyzing 2D acoustic problems by using quadrilateral element. In the present EDM, the quadrilateral element is first subdivided into four sub-triangles, and the local acoustic gradient in each sub-triangle is obtained using linear interpolation function. The acoustic gradient field of the whole quadrilateral is then formulated through a weighted averaging operation, which means only one integration point is adopted to construct the system matrix. To cure the numerical instability of one-point integration, a variation gradient item is complemented by variance of the local gradients. The discretized system equations are derived using the generalized Galerkin weakform. Numerical examples demonstrate that the EDM can achieves better accuracy and higher computational efficiency. Besides, as no mapping or coordinate transformation is involved, restrictions on the shape elements can be easily removed, which makes the EDM works well even for severely distorted meshes.

In this paper, we present accurate and economic integration quadratures for hypersingular functions over three simple geometric shapes in ℝ3 (spheres, cubes, and cylinders). The quadrature nodes are made of the tensor-product of 1-D Gauss nodes on [–1,1] for non-periodic variables or uniform nodes on [0,2π] or [0,π] for periodic ones. The quadrature weights are converted from a brute-force integration of the hypersingular function through interpolating the smooth component of the integrand. Numerical results are presented to validate the accuracy and efficiency of computing hypersingular integrals, as in the computations of Cauchy principal values, with a minimum number of quadrature nodes. The pre-calculated quadrature tables can be then readily used to implement Nyström collocation methods of hypersingular volume integral equations such as the one for Maxwell equations.

In this paper, we develop a novel energy-preserving wavelet collocation method for solving general multi-symplectic formulations of Hamiltonian PDEs. Based on the autocorrelation functions of Daubechies compactly supported scaling functions, the wavelet collocation method is conducted for spatial discretization. The obtained semi-discrete system is shown to be a finite-dimensional Hamiltonian system, which has an energy conservation law. Then, the average vector field method is used for time integration, which leads to an energy-preserving method for multi-symplectic Hamiltonian PDEs. The proposed method is illustrated by the nonlinear Schrödinger equation and the Camassa-Holm equation. Since differentiation matrix obtained by the wavelet collocation method is a cyclic matrix, we can apply Fast Fourier transform to solve equations in numerical calculation. Numerical experiments show the high accuracy, effectiveness and conservation properties of the proposed method.

In this paper, we study the numerical solution of singularly perturbed time-dependent convection-diffusion problems. To solve these problems, the backward Euler method is first applied to discretize the time derivative on a uniform mesh, and the classical upwind finite difference scheme is used to approximate the spatial derivative on an arbitrary nonuniform grid. Then, in order to obtain an adaptive grid for all temporal levels, we construct a positive monitor function, which is similar to the arc-length monitor function. Furthermore, the ε-uniform convergence of the fully discrete scheme is derived for the numerical solution. Finally, some numerical results are given to support our theoretical results.

We propose an implicit finite-difference method to study the time evolution of the director field of a nematic liquid crystal under the influence of an electric field with weak anchoring at the boundary. The scheme allows us to study the dynamics of transitions between different director equilibrium states under varying electric field and anchoring strength. In particular, we are able to simulate the transition to excited states of odd parity, which have previously been observed in experiments, but so far only analyzed in the static case.

Even though there are various fast methods and preconditioning techniques available for the simulation of Poisson problems, little work has been done for solving Poisson's equation by using the Helmholtz decomposition scheme. To bridge this issue, we propose a novel efficient algorithm to solve Poisson's equation in irregular two dimensional domains for electrostatics through a quasi-Helmholtz decomposition technique—the loop-tree basis decomposition. It can handle Dirichlet, Neumann or mixed boundary problems in which the filling media can be homogeneous or inhomogeneous. A novel point of this method is to first find the electric flux efficiently by applying the loop-tree basis functions. Subsequently, the potential is obtained by finding the inverse of the gradient operator. Furthermore, treatments for both Dirichlet and Neumann boundary conditions are addressed. Finally, the validation and efficiency are illustrated by several numerical examples. Through these simulations, it is observed that the computational complexity of our proposed method almost scales as , where N is the triangle patch number of meshes. Consequently, this new algorithm is a feasible fast Poisson solver.

In this paper, a conservative parallel iteration scheme is constructed to solve nonlinear diffusion equations on unstructured polygonal meshes. The design is based on two main ingredients: the first is that the parallelized domain decomposition is embedded into the nonlinear iteration; the second is that prediction and correction steps are applied at subdomain interfaces in the parallelized domain decomposition method. A new prediction approach is proposed to obtain an efficient conservative parallel finite volume scheme. The numerical experiments show that our parallel scheme is second-order accurate, unconditionally stable, conservative and has linear parallel speed-up.

Nanosecond (ns) pulsed dielectric barrier discharge (DBD) actuator in a laminar flat plate boundary layer is investigated numerically in an attempt to gain some new insights into the understanding of ns DBD actuation mechanism. Special emphasis is put on the examination, separation and comparison of behaviors of discharge induced micro shock wave and residual heat as well as on the investigation of response of external flow to the two effects. The shock wave is found to introduce highly transient, localized perturbation to the flow and be able to significantly alter the flow pattern shortly after its initiation. The main flow tends to quickly recover to close to its undisturbed state due to the transient nature of perturbation. However, with the shock decay and final disappearance, another perturbation source in the vicinity of discharge region, which contains contribution from both residual heat and shock, becomes increasingly pronounced and eventually develops into a perturbation wave train in the boundary layer. The perturbation is relatively weak and may not be a Tollmien-Schlichting (TS) wave and not trigger the laminar-turbulent transition of boundary layer. Instead, it is more likely to manipulate the flow stability to achieve the strong control authority of this kind of actuation in the case of flow separation control. In addition, a parametric study over the different electrical and hydrodynamic parameters is also conducted.

The radiative transfer equation (RTE) arises in many different areas of science and engineering. In this paper, we propose and investigate a discrete-ordinate discontinuous-streamline diffusion (DODSD) method for solving the RTE, which is a combination of the discrete-ordinate technique and the discontinuous-streamline diffusion method. Different from the discrete-ordinate discontinuous Galerkin (DODG) method for the RTE, an artificial diffusion parameter is added to the test functions in the spatial discretization. Stability and error estimates in certain norms are proved. Numerical results show that the proposed method can lead to a more accurate approximation in comparison with the DODG method.