Local approximate radiation boundary conditions of optimal efficiency for the convective wave equation and the linearized Euler equations in waveguide geometry are formulated, analyzed, and tested. The results extend and improve for the convective case the general formulation of high-order local radiation boundary condition sequences for anisotropic scalar equations developed in [4].

We apply the immersed boundary (or IB) method to simulate deformation and detachment of a periodic array of wall-bounded biofilm colonies in response to a linear shear flow. The biofilm material is represented as a network of Hookean springs that are placed along the edges of a triangulation of the biofilm region. The interfacial shear stress, lift and drag forces acting on the biofilm colony are computed by using fluid stress jump method developed by Williams, Fauci and Gaver [Disc. Con-tin. Dyn. Sys. B 11(2):519–540, 2009], with a modified version of their exclusion filter. Our detachment criterion is based on the novel concept of an averaged equivalent continuum stress tensor defined at each IB point in the biofilm which is then used to determine a corresponding von Mises yield stress; wherever this yield stress exceeds a given critical threshold the connections to that node are severed, thereby signalling the onset of a detachment event. In order to capture the deformation and detachment behaviour of a biofilm colony at different stages of growth, we consider a family of four biofilm shapes with varying aspect ratio. For each aspect ratio, we varied the spacing between colonies to investigate role of spatial clustering in offering protection against detachment. Our numerical simulations focus on the behaviour of weak biofilms (with relatively low yield stress threshold) and investigate features of the fluid-structure interaction such as locations of maximum shear and increased drag. The most important conclusions of this work are: (a) reducing the spacing between colonies reduces drag by from 50 to 100% and alters the interfacial shear stress profile, suggesting that even weak biofilms may be able to grow into tall structures because of the protection they gain from spatial proximity with other colonies; (b) the commonly employed detachment strategy in biofilm models based only on interfacial shear stress can lead to incorrect or inaccurate results when applied to the study of shear induced detachment of weak biofilms. Our detachment strategy based on equivalent continuum stresses provides a unified and consistent IB framework that handles both sloughing and erosion modes of biofilm detachment, and is consistent with strategies employed in many other continuum based biofilm models.

In this contribution, a new class of lattice Boltzmann schemes is introduced and studied. These schemes are presented in a framework that generalizes the multiple relaxation times method of d’Humières. They extend also the Geier’s cascaded method. The relaxation phase takes place in a moving frame involving a set of moments depending on a given relative velocity field. We establish with the Taylor expansion method that the equivalent partial differential equations are identical to the ones obtained with the multiple relaxation times method up to the second order accuracy. The method is then performed to derive the equivalent equations up to third order accuracy.

Numerical methods of a 3D multiphysics, two-phase transport model of proton exchange membrane fuel cell (PEMFC) is studied in this paper. Due to the coexistence of multiphase regions, the standard finite element/finite volume method may fail to obtain a convergent nonlinear iteration for a two-phase transport model of PEMFC [49,50]. By introducing Kirchhoff transformation technique and a combined finite element-upwind finite volume approach, we efficiently achieve a fast convergence and reasonable solutions for this multiphase, multiphysics PEMFC model. Numerical implementation is done by using a novel automated finite element/finite volume program generator (FEPG). By virtue of a high-level algorithm description language (script), component programming and human intelligence technologies, FEPG can quickly generate finite element/finite volume source code for PEMFC simulation. Thus, one can focus on the efficient algorithm research without being distracted by the tedious computer programming on finite element/finite volume methods. Numerical success confirms that FEPG is an efficient tool for both algorithm research and software development of a 3D, multiphysics PEMFC model with multicomponent and multiphase mechanism.

Motivated by the magneto hydrodynamic (MHD) simulation for Tokamaks with Isogeometric analysis, we present splines defined over a rectangular mesh with a complex topological structure, i.e., with extraordinary vertices. These splines are piecewise polynomial functions of bi-degree (d,d) and parameter continuity. And we compute their dimension and exhibit basis functions called Hermite bases for bicubic spline spaces. We investigate their potential applications for solving partial differential equations (PDEs) over a physical domain in the framework of Isogeometric analysis. For instance, we analyze the property of approximation of these spline spaces for the L 2-norm; we show that the optimal approximation order and numerical convergence rates are reached by setting a proper parameterization, although the fact that the basis functions are singular at extraordinary vertices.

Moment equations provide a flexible framework for the approximation of the Boltzmann equation in kinetic gas theory. While moments up to second order are sufficient for the description of equilibrium processes, the inclusion of higher order moments, such as the heat flux vector, extends the validity of the Euler equations to non-equilibrium gas flows in a natural way.

Unfortunately, the classical closure theory proposed by Grad leads to moment equations, which suffer not only from a restricted hyperbolicity region but are also affected by non-physical sub-shocks in the continuous shock-structure problem if the shock velocity exceeds a critical value. Amore recently suggested closure theory based on the maximum entropy principle yields symmetric hyperbolic moment equations. However, if moments higher than second order are included, the computational demand of this closure can be overwhelming. Additionally, it was shown for the 5-moment system that the closing flux becomes singular on a subset of moments including the equilibrium state.

Motivated by recent promising results of closed-form, singular closures based on the maximum entropy approach, we study regularized singular closures that become singular on a subset of moments when the regularizing terms are removed. In order to study some implications of singular closures, we use a recently proposed explicit closure for the 5-moment equations. We show that this closure theory results in a hyperbolic system that can mitigate the problem of sub-shocks independent of the shock wave velocity and handle strongly non-equilibrium gas flows.

In this work, we concern with the numerical comparison between different kinds of design points in least square (LS) approach on polynomial spaces. Such a topic is motivated by uncertainty quantification (UQ). Three kinds of design points are considered, which are the Sparse Grid (SG) points, the Monte Carlo (MC) points and the Quasi Monte Carlo (QMC) points. We focus on three aspects during the comparison: (i) the convergence properties; (ii) the stability, i.e. the properties of the resulting condition number of the design matrix; (iii) the robustness when numerical noises are present in function values. Several classical high dimensional functions together with a random ODE model are tested. It is shown numerically that (i) neither the MC sampling nor the QMC sampling introduce the low convergence rate, namely, the approach achieves high order convergence rate for all cases provided that the underlying functions admit certain regularity and enough design points are used; (ii)The use of SG points admits better convergence properties only for very low dimensional problems (say d ≤ 2); (iii)The QMC points, being deterministic, seem to be a good choice for higher dimensional problems not only for better convergence properties but also in the stability point of view.

A new class of supersonic nonequilibrium flows is studied on the basis of solving the Boltzmann and model kinetic equations with the aim to consider new non-linear structures in open systems and to study anomalous transfer properties in relaxation zones. The Unified Flow Solver is applied for numerical simulations. Simple gases and gases with inner degrees of freedom are considered. The experimental data related to the influence of the so-called optical lattices on the supersonic molecular beams are considered and numerical analysis of the nonequilibrium states obtained on this basis is made. The nonuniform relaxation problem with these distributions is simulated numerically and anomalous transport is confirmed. The conditions for strong changes of the temperature in the anomalous transfer zones are discussed and are realized in computations.

In this paper we present a 2D/3D high order accurate finite volume scheme in the context of direct Arbitrary-Lagrangian-Eulerian algorithms for general hyperbolic systems of partial differential equations with non-conservative products and stiff source terms. This scheme is constructed with a single stencil polynomial reconstruction operator, a one-step space-time ADER integration which is suitably designed for dealing even with stiff sources, a nodal solver with relaxation to determine the mesh motion, a path-conservative integration technique for the treatment of non-conservative products and an a posteriori stabilization procedure derived from the so-called Multidimensional Optimal Order Detection (MOOD) paradigm. In this work we consider the seven equation Baer-Nunziato model of compressible multi-phase flows as a representative model involving non-conservative products as well as relaxation source terms which are allowed to become stiff. The new scheme is validated against a set of test cases on 2D/3D unstructured moving meshes on parallel machines and the high order of accuracy achieved by the method is demonstrated by performing a numerical convergence study. Classical Riemann problems and explosion problems with exact solutions are simulated in 2D and 3D. The overall numerical code is also profiled to provide an estimate of the computational cost required by each component of the whole algorithm.

In this short note we present a derivation of the Spectral Difference Scheme from a Discontinuous Galerkin (DG) discretization of a nonlinear conservation law. This allows interpretation of the Spectral Difference Scheme as a particular discretization under the quadrature-free nodal DG paradigm. Moreover, it enables identification of the key differences between the Spectral Difference Scheme and standard nodal DG schemes.

To finite-difference model elastic wave propagation in a combined structure with solid, fluid and porous subregions, a set of modified Biot’s equations are used, which can be reduced to the governing equations in solids, fluids as well as fluid-saturated porous media. Based on the modified Biot’s equations, the field quantities are finite-difference discretized into unified forms in the whole structure, including those on any interface between the solid, fluid and porous subregions. For the discrete equations on interfaces, however, the harmonic mean of shear modulus and the arithmetic mean of the other parameters on both sides of the interfaces are used. These parameter averaging equations are validated by deriving from the continuity conditions on the interfaces. As an example of using the parameter averaging technique, a 2-D finite-difference scheme with a velocity-stress staggered grid in cylindrical coordinates is implemented to simulate the acoustic logs in porous formations. The finite-difference simulations of the acoustic logging in a homogeneous formation agree well with those obtained by the analytical method. The acoustic logs with mud cakes clinging to the borehole well are simulated for investigating the effect of mud cake on the acoustic logs. The acoustic logs with a varying radius borehole embedded in a horizontally stratified formation are also simulated by using the proposed finite-difference scheme.

We study the static and dynamical behavior of the contact line between two fluids and a solid plate by means of the Lattice Boltzmann method (LBM). The different fluid phases and their contact with the plate are simulated by means of standard Shan-Chen models. We investigate different regimes and compare the multicomponent vs. the multiphase LBM models near the contact line. A static interface profile is attained with the multiphase model just by balancing the hydrostatic pressure (due to gravity) with a pressure jump at the bottom. In order to study the same problem with the multicomponent case we propose and validate an idea of a body force acting only on one of the two fluid components. In order to reproduce results matching an infinite bath, boundary conditions at the bath side play a key role. We quantitatively compare open and wall boundary conditions and study their influence on the shape of the meniscus against static and lubrication theory solution.

In this paper, we develop, analyze and test local discontinuous Galerkin (LDG) methods for solving the Degasperis-Procesi equation which contains nonlinear high order derivatives, and possibly discontinuous or sharp transition solutions. The LDG method has the flexibility for arbitrary h and p adaptivity. We prove the L2 stability for general solutions. The proof of the total variation stability of the schemes for the piecewise constant P0 case is also given. The numerical simulation results for different types of solutions of the nonlinear Degasperis-Procesi equation are provided to illustrate the accuracy and capability of the LDG method.

In this paper, we numerically investigate the effects of surfactant on drop-drop interactions in a 2D shear flow using a coupled level-set and immersed interface approach proposed in (Xu et al., J. Comput. Phys., 212 (2006), 590-616). We find that surfactant plays a critical and nontrivial role in drop-drop interactions. In particular, we find that the minimum distance between the drops is a non-monotone function of the surfactant coverage and Capillary number. This non-monotonic behavior, which does not occur for clean drops, is found to be due to the presence of Marangoni forces along the drop interfaces. This suggests that there are non-monotonic conditions for coalescence of surfactant-laden drops, as observed in recent experiments of Leal and co-workers. Although our study is two-dimensional, we believe that drop-drop interactions in three-dimensional flows should be qualitatively similar as the Maragoni forces in the near contact region in 3D should have a similar effect.

The computation of the radiative transfer equation is expensive mainly due to two stiff terms: the transport term and the collision operator. The stiffness in the former comes from the fact that particles (such as photons) travel at the speed of light, while that in the latter is due to the strong scattering in the optically thick region. We study the fully implicit scheme for this equation to account for the stiffness. The main challenge in the implicit treatment is the coupling between the spacial and angular coordinates that requires the large size of the to-be-inverted matrix, which is also ill-conditioned and not necessarily symmetric. Our main idea is to utilize the spectral structure of the ill-conditioned matrix to construct a pre-conditioner, which, along with an exquisite split of the spatial and angular dependence, significantly improve the condition number and allows a matrix-free treatment. We also design a fast solver to compute this pre-conditioner explicitly in advance. Our method is shown to be efficient in both diffusive and free streaming limit, and the computational cost is comparable to the state-of-the-art method. Various examples including anisotropic scattering and two-dimensional problems are provided to validate the effectiveness of our method.

The development of high-order schemes has been mostly concentrated on the limiters and high-order reconstruction techniques. In this paper, the effect of the flux functions on the performance of high-order schemes will be studied. Based on the same WENO reconstruction, two schemes with different flux functions, i.e., the fifth-order WENO method and the WENO-Gas-kinetic scheme (WENO-GKS), will be compared. The fifth-order finite difference WENO-SW scheme is a characteristic variable reconstruction based method which uses the Steger-Warming flux splitting for inviscid terms, the sixth-order central difference for viscous terms, and three stages Runge-Kutta time stepping for the time integration. On the other hand, the finite volume WENO-GKS is a conservative variable reconstruction based method with the same WENO reconstruction. But, it evaluates a time dependent gas distribution function along a cell interface, and updates the flow variables inside each control volume by integrating the flux function along the boundary of the control volume in both space and time. In order to validate the robustness and accuracy of the schemes, both methods are tested under a wide range of flow conditions: vortex propagation, Mach 3 step problem, and the cavity flow at Reynolds number 3200. Our study shows that both WENO-SW and WENO-GKS yield quantitatively similar results and agree with each other very well provided a sufficient grid resolution is used. With the reduction of mesh points, the WENO-GKS behaves to have less numerical dissipation and present more accurate solutions than those from the WENO-SW in all test cases. For the Navier-Stokes equations, since the WENO-GKS couples inviscid and viscous terms in a single flux evaluation, and the WENO-SW uses an operator splitting technique, it appears that the WENO-SW is more sensitive to the WENO reconstruction and boundary treatment. In terms of efficiency, the finite volume WENO-GKS is about 4 times slower than the finite differenceWENO-SW in two dimensional simulations. The current study clearly shows that besides high-order reconstruction, an accurate gas evolution model or flux function in a high-order scheme is also important in the capturing of physical solutions. In a physical flow, the transport, stress deformation, heat conduction, and viscous heating are all coupled in a single gas evolution process. Therefore, it is preferred to develop such a scheme with multi-dimensionality, and unified treatment of inviscid and dissipative terms. A high-order scheme does prefer a high-order gas evolution model. Even with the rapid advances of high-order reconstruction techniques, the first-order dynamics of the Riemann solution becomes the bottleneck for the further development of high-order schemes. In order to avoid the weakness of the low order flux function, the development of high-order schemes relies heavily on the weak solution of the original governing equations for the update of additional degree of freedom, such as the non-conservative gradients of flow variables, which cannot be physically valid in discontinuous regions.

The Jiles-Atherton (J-A) model is a commonly used physics-based model in describing the hysteresis characteristics of ferromagnetic materials. However, citations and interpretation of this model in literature have been non-uniform. Solution methods for solving numerically this model has not been studied adequately. In this paper, through analyzing the mathematical properties of equations and the physical mechanism of energy conservation, we point out some unreasonable descriptions of this model and develop a relatively more accurate, modified J-A model together with its numerical solution method. Our method employs a fixed point method to compute anhysteretic magnetization. We obtain the susceptibility value of the anhysteretic magnetization analytically and apply the 4th order Runge-Kutta method to the solution of total magnetization. Computational errors are estimated and then precisions of the solving method in describing various materials are verified. At last, through analyzing the effects of the accelerating method, iterative error and step size on the computational errors, we optimize the numerical method to achieve the effects of high precision and short computing time. From analysis, we determine the range of best values of some key parameters for fast and accurate computation.

In multi physics computations where a compressible fluid is coupled with a linearly elastic solid, it is standard to enforce continuity of the normal velocities and of the normal stresses at the interface between the fluid and the solid. In a numerical scheme, there are many ways that velocity- and stress-continuity can be enforced in the discrete approximation. This paper performs a normal mode stability analysis of the linearized problem to investigate the stability of different numerical interface conditions for a model problem approximated by upwind type finite difference schemes. The analysis shows that depending on the ratio of densities between the solid and the fluid, some numerical interface conditions are stable up to the maximal CFL-limit, while other numerical interface conditions suffer from a severe reduction of the stable CFL-limit. The paper also presents a new interface condition, obtained as a simplified characteristic boundary condition, that is proved to not suffer from any reduction of the stable CFL-limit. Numerical experiments in one space dimension show that the new interface condition is stable also for computations with the non-linear Euler equations of compressible fluid flow coupled with a linearly elastic solid.

A second order Ghost Fluid method is proposed for the treatment of interface problems of elliptic equations with discontinuous coefficients. By appropriate use of auxiliary virtual points, physical jump conditions are enforced at the interface. The signed distance function is used for the implicit description of irregular domain. With the additional unknowns, high order approximation considering the discontinuity can be built. To avoid the ill-conditioned matrix, the interpolation stencils are selected adaptively to balance the accuracy and the numerical stability. Additional equations containing the jump restrictions are assembled with the original discretized algebraic equations to form a new sparse linear system. Several Krylov iterative solvers are tested for the newly derived linear system. The results of a series of 1-D, 2-D tests show that the proposed method possesses second order accuracy in L∞ norm. Besides, the method can be extended to the 3-D problems straightforwardly. Numerical results reveal the present method is highly efficient and robust in dealing with the interface problems of elliptic equations.

High intensity focused ultrasound is a rapidly developing technology for the ablation of tumors. Liver cancer is one of the most common malignancies worldwide. Since liver has a large number of blood vessels, blood flow cooling can reduce the necrosed volume and may cause regeneration of the tumor to occur. All cancer cells should be ablated without damaging of the critical tissues. Today, treatment planning tools consider liver as a homogeneous organ. This paper is a step towards the development of surgical planning platform for a non-invasive HIFU tumor ablative therapy in a real liver geometry based on CT/MRI image. This task requires coupling of different physical fields: acoustic, thermal and hydrodynamic. These physical fields can influence each other. In this paper we illustrate how a computational model can be used to improve the treatment efficiency. In large blood vessel both convective cooling and acoustic streaming can change the temperature considerably near blood vessel. The whole tumor ablation took only 30 seconds in the considered simulation case, which is very small comparing with the current treatment time of several hours. Through this study we are convinced that high ultrasound power and nonlinear propagation effects with appropriate treatment planning can sufficiently reduce the treatment time.