In this work, a decoupled, parallel, iterative finite element method for solving the steady Boussinesq equations is proposed and analyzed. Starting from an initial guess, an iterative algorithm is designed to decouple the Naiver-Stokes equations and the heat equation based on certain explicit treatment with the solution from the previous iteration step. At each step of the iteration, the two equations can be solved in parallel by using finite element discretization. The existence and uniqueness of the solution to each step of the algorithm is proved. The stability analysis and error estimation are also carried out. Numerical tests are presented to verify the analysis results and illustrate the applicability of the proposed method.

In [15], the computational performance of various weak Galerkin finite element methods in terms of stability, convergence, and supercloseness is explored and numerical results are listed in 31 tables. Some of the phenomena can be explained by the existing theoretical results and the others are to be explained. The main purpose of this paper is to provide a unified theoretical foundation to a class of WG schemes, where $(P_k(T), P_{k+1}(e), [P_{k+1}(T)]^2)$ elements are used for solving the second order elliptic equations (1)-(2) on a triangle grid in 2D. With this unified treatment, all of the existing results become special cases. The theoretical conclusions are corroborated by a number of numerical examples.

In this paper, we investigate numerically the long-time dynamics of a two-dimensional dendritic precipitate. We focus our study on the self-similar scaling behavior of the primary dendritic arm with profile $x∼t^{α_1}$ and $y∼t^{α_2},$ and explore the dependence of parameters $α_1$ and $α_2$ on applied driving forces of the system (e.g. applied far-field flux or strain). We consider two dendrite forming mechanisms: the dendritic growth driven by (i) an anisotropic surface tension and (ii) an applied strain at the far-field of the elastic matrix. We perform simulations using a spectrally accurate boundary integral method, together with a rescaling scheme to speed up the intrinsically slow evolution of the precipitate. The method enables us to accurately compute the dynamics far longer times than could previously be accomplished. Comparing with the original work on the scaling behavior $α_1 = 0.6$ and $α_2 = 0.4$ [Phys. Rev. Lett. 71(21) (1993) 3461–3464], where a constant flux was used in a diffusion only problem, we found at long times this scaling still serves a good estimation of the dynamics though it deviates from the asymptotic predictions due to slow retreats of the dendrite tip at later times. In particular, we find numerically that the tip grows self-similarly with $α_1 = 1/3$ and $α_2 = 1/3$ if the driving flux $J ∼ 1/R(t)$ where $R(t)$ is the equivalent size of the evolving precipitate. In the diffusive growth of precipitates in an elastic media, we examine the tip of the precipitate under elastic stress, under both isotropic and anisotropic surface tension, and find that the tip also follows a scaling law.

In the phase-field modeling of brittle fracture, anisotropic constitutive assumptions for the degradation of stored elastic energy due to fracture are crucial to preventing cracking in compression and obtaining physically sound numerical solutions. Three energy decomposition models, the spectral decomposition, the volumetric-deviatoric split, and a modified volumetric-deviatoric split, and their effects on the performance of the phase-field modeling are studied. Meanwhile, anisotropic degradation of stiffness may lead to a small amount of energy remaining on crack surfaces, which violates crack boundary conditions and can cause unphysical crack openings and propagation. A simple yet effective treatment for this is proposed: define a critically damaged zone with a threshold parameter and then degrade both the active and passive energies in the zone. A dynamic mesh adaptation finite element method is employed for the numerical solution of the corresponding elasticity system. Four examples, including two benchmark ones, one with complex crack systems, and one based on an experimental setting, are considered. Numerical results show that the spectral decomposition and modified volumetric-deviatoric split models, together with the improvement treatment of crack boundary conditions, can lead to crack propagation results that are comparable with the existing computational and experimental results. It is also shown that the numerical results are not sensitive to the parameter defining the critically damaged zone.

Turbulent mixing due to hydrodynamic instabilities occurs in a wide range of science and engineering applications such as supernova explosions and inertial confinement fusion. The experimental, theoretical and numerical studies help us to understand the dynamics of hydrodynamically unstable interfaces between fluids in these important problems. In this paper, we present an increasingly accurate and robust front tracking method for the numerical simulations of Richtmyer-Meshkov Instability (RMI) to estimate the growth rate. The single-mode shock tube experiments of Collins and Jacobs 2002 [1] for two incident shock strengths $(M = 1.11$ and $M = 1.21)$ are used to validate the RMI simulations. The simulations based on the classical fifth order weighted essentially non-oscillatory (WENO) scheme of Jiang and Shu [2] with Yang’s artificial compression [3] are compared with Collins and Jacobs 2002 shock tube experiments. We investigate the resolution effects using front tracking with WENO schemes on the two-dimensional RMI of an ${\rm air}/SF_6$ interface. We achieve very good agreement on the early time interface displacement and amplitude growth rate between simulations and experiments for Mach number $M = 1.11.$ A 4% discrepancy on early-time amplitude is observed between the fine grid simulation and the $M = 1.21$ experiments of Collins and Jacobs 2002.

In this paper we provide a detailed analysis of the preconditioned steepest descent (PSD) iteration solver for a convex splitting numerical scheme to the Cahn-Hilliard equation with variable mobility function. In more details, the convex-concave decomposition is applied to the energy functional, which in turn leads to an implicit treatment for the nonlinear term and the surface diffusion term, combined with an explicit update for the expansive concave term. In addition, the mobility function, which is solution-dependent, is explicitly computed, which ensures the elliptic property of the operator associated with the temporal derivative. The unique solvability of the numerical scheme is derived following the standard convexity analysis, and the energy stability analysis could also be carefully established. On the other hand, an efficient implementation of the numerical scheme turns out to be challenging, due to the coupling of the nonlinear term, the surface diffusion part, and a variable-dependent mobility elliptic operator. Since the implicit parts of the numerical scheme are associated with a strictly convex energy, we propose a preconditioned steepest descent iteration solver for the numerical implementation. Such an iteration solver consists of a computation of the search direction (involved with a Poisson-like equation), and a one-parameter optimization over the search direction, in which the Newton’s iteration becomes very powerful. In addition, a theoretical analysis is applied to the PSD iteration solver, and a geometric convergence rate is proved for the iteration. A few numerical examples are presented to demonstrate the robustness and efficiency of the PSD solver.

In this paper, we present and analyze a new ultra-weak discontinuous Galerkin (UWDG) finite element method for two-dimensional semilinear second-order elliptic problems on Cartesian grids. Unlike the traditional local discontinuous Galerkin (LDG) method, the proposed UWDG method can be applied without introducing any auxiliary variables or rewriting the original equation into a system of equations. The UWDG scheme is presented in details, including the definition of the numerical fluxes, which are necessary to obtain optimal error estimates. The proposed scheme can be made arbitrarily high-order accurate in two-dimensional space. The error estimates of the presented scheme are analyzed. The order of convergence is proved to be $p + 1$ in the $L^2$-norm, when tensor product polynomials of degree at most $p$ and grid size $h$ are used. Several numerical examples are provided to confirm the theoretical results.

A finite element solution of an ion channel dielectric continuum model such as Poisson-Boltzmann equation (PBE) and a system of Poisson-Nernst-Planck equations (PNP) requires tetrahedral meshes for an ion channel protein region, a membrane region, and an ionic solvent region as well as an interface fitted irregular tetrahedral mesh of a simulation box domain. However, generating these meshes is very difficult and highly technical due to the related three regions having very complex geometrical shapes. Currently, an ion channel mesh generation software package developed in Lu’s research group is the only one available in the public domain. To significantly improve its mesh quality and computer performance, in this paper, new numerical schemes for generating membrane and solvent meshes are presented and implemented in Python, resulting in a new ion channel mesh generation software package. Numerical results are then reported to demonstrate the efficiency of the new numerical schemes and the quality of meshes generated by the new package for ion channel proteins with ion channel pores having different geometric complexities.