K. B. Nakshatrala

Numerical Method for Darcy Flow Derived Using Discrete Exterior Calculus

We derive a numerical method for Darcy flow, and also for Poisson’s equation in mixed (first order) form, based on discrete exterior calculus (DEC). Exterior calculus is a generalization of vector calculus to smooth manifolds and DEC is one of its discretizations on simplicial complexes such as triangle and tetrahedral meshes. DEC is a coordinate invariant discretization, in that it does not depend on the embedding of the simplices or the whole mesh. We start by rewriting the governing equations of Darcy flow using the language of exterior calculus. This yields a formulation in terms of flux differential form and pressure. The numerical method is then derived by using the framework provided by DEC for discretizing differential forms and operators that act on forms. We also develop a discretization for a spatially dependent Hodge star that varies with the permeability of the medium. This also allows us to address discontinuous permeability. The matrix representation for our discrete non-homogeneous Hodge star is diagonal, with positive diagonal entries. The resulting linear system of equations for flux and pressure are saddle type, with a diagonal matrix as the top left block. The performance of the proposed numerical method is illustrated on many standard test problems. These include patch tests in two and three dimensions, comparison with analytically known solutions in two dimensions, layered medium with alternating permeability values, and a test with a change in permeability along the flow direction. We also show numerical evidence of convergence of the flux and the pressure. A convergence experiment is included for Darcy flow on a surface. A short introduction to the relevant parts of smooth and discrete exterior calculus is included in this article. We also include a discussion of the boundary condition in terms of exterior calculus.

A numerical framework for diffusion-controlled bimolecular-reactive systems to enforce maximum principles and the non-negative constraint

We present a novel computational framework for diffusive-reactive systems that satisfies the non-negative constraint and maximum principles on general computational grids. The governing equations for the concentration of reactants and product are written in terms of tensorial diffusion-reaction equations. We restrict our studies to fast irreversible bimolecular reactions. If one assumes that the reaction is diffusion-limited and all chemical species have the same diffusion coefficient, one can employ a linear transformation to rewrite the governing equations in terms of invariants, which are unaffected by the reaction. This results in two uncoupled tensorial diffusion equations in terms of these invariants, which are solved using a novel non-negative solver for tensorial diffusion-type equations. The concentrations of the reactants and the product are then calculated from invariants using algebraic manipulations. The novel aspect of the proposed computational framework is that it will always produce physically meaningful non-negative values for the concentrations of all chemical species. Several representative numerical examples are presented to illustrate the robustness, convergence, and the numerical performance of the proposed computational framework. We will also compare the proposed framework with other popular formulations. In particular, we will show that the Galerkin formulation (which is the standard single-field formulation) does not produce reliable solutions, and the reason can be attributed to the fact that the single-field formulation does not guarantee non-negative solutions. We will also show that the clipping procedure (which produces non-negative solutions but is considered as a variational crime) does not give accurate results when compared with the proposed computational framework.

On multi-time-step monolithic coupling algorithms for elastodynamics

We present a way of constructing multi-time-step monolithic coupling methods for elastodynamics. The governing equations for constrained multiple subdomains are written in dual Schur form and the continuity of velocities is enforced at system time levels. The resulting equations will be in the form of differential-algebraic equations. To crystallize the ideas we shall employ Newmark family of time-stepping schemes. The proposed method can handle multiple subdomains, and allows different time-steps as well as different time stepping schemes from the Newmark family in different subdomains. We shall use the energy method to assess the numerical stability, and quantify the influence of perturbations under the proposed coupling method. Two different notions of energy preservation are introduced and employed to assess the performance of the proposed method. Several numerical examples are presented to illustrate the accuracy and stability properties of the proposed method. We shall also compare the proposed multi-time-step coupling method with some other methods available in the literature.

On mesh restrictions to satisfy comparison principles, maximum principles, and the non-negative constraint: Recent developments and new results

ABSTRACT This article concerns mesh restrictions that are needed to satisfy several important mathematical properties—maximum principles, comparison principles, and the nonnegative constraint—for a general linear second-order elliptic partial differential equation. We critically review some recent developments in the field of discrete maximum principles, derive new results, and discuss some possible future research directions in this area. In particular, we derive restrictions for a three-node triangular (T3) element and a four-node quadrilateral (Q4) element to satisfy comparison principles, maximum principles, and the nonnegative constraint under the standard single-field Galerkin formulation. Analysis is restricted to uniformly elliptic linear differential operators in divergence form with Dirichlet boundary conditions specified on the entire boundary of the domain. Various versions of maximum principles and comparison principles are discussed in both continuous and discrete settings. In the literature, it is well-known that an acute-angled triangle is sufficient to satisfy the discrete weak maximum principle for pure isotropic diffusion. Herein, we show that this condition can be either too restrictive or not sufficient to satisfy various discrete principles when one considers anisotropic diffusivity, advection velocity field, or linear reaction coefficient. Subsequently, we derive appropriate restrictions on the mesh for simplicial (e.g., T3 element) and nonsimplicial (e.g., Q4 element) elements. Based on these conditions, an iterative algorithm is developed to construct simplicial meshes that preserve discrete maximum principles using existing open source mesh generators. Various numerical examples based on different types of triangulations are presented to show the pros and cons of placing restrictions on a computational mesh. We also quantify local and global mass conservation errors using representative numerical examples and illustrate the performance of metric.

A monolithic multi-time-step computational framework for first-order transient systems with disparate scales

Abstract Developing robust simulation tools for problems involving multiple mathematical scales has been a subject of great interest in computational mathematics and engineering. A desirable feature to have in a numerical formulation for multiscale transient problems is to be able to employ different time-steps (multi-time-step coupling), and different time integrators and different numerical formulations (mixed methods) in different regions of the computational domain. To this end, we present two new monolithic multi-time-step mixed coupling methods for first-order transient systems . We shall employ unsteady advection–diffusion–reaction equation with linear decay as the model problem, which offers several unique challenges in terms of non-self-adjoint spatial operator and rich features in the solutions. We shall employ the dual Schur domain decomposition technique to split the computational domain into an arbitrary number of subdomains. It will be shown that the governing equations of the decomposed problem, after spatial discretization, will be differential/algebraic equations. This is a crucial observation to obtain stable numerical results. Two different methods of enforcing compatibility along the subdomain interface will be used in the time discrete setting. A systematic theoretical analysis (which includes numerical stability, influence of perturbations, bounds on drift along the subdomain interface) will be performed. The first coupling method ensures that there is no drift along the subdomain interface, but does not facilitate explicit/implicit coupling. The second coupling method allows explicit/implicit coupling with controlled (but non-zero) drift in the solution along the subdomain interface. Several canonical problems will be solved to numerically verify the theoretical predictions, and to illustrate the overall performance of the proposed coupling methods. Finally, we shall illustrate the robustness of the proposed coupling methods using a multi-time-step transient simulation of a fast bimolecular advective–diffusive–reactive system.

Large-Scale Optimization-Based Non-negative Computational Framework for Diffusion Equations: Parallel Implementation and Performance Studies

It is well-known that the standard Galerkin formulation, which is often the formulation of choice under the finite element method for solving self-adjoint diffusion equations, does not meet maximum principles and the non-negative constraint for anisotropic diffusion equations. Recently, optimization-based methodologies that satisfy maximum principles and the non-negative constraint for steady-state and transient diffusion-type equations have been proposed. To date, these methodologies have been tested only on small-scale academic problems. The purpose of this paper is to systematically study the performance of the non-negative methodology in the context of high performance computing (HPC). PETSc and TAO libraries are, respectively, used for the parallel environment and optimization solvers. For large-scale problems, it is important for computational scientists to understand the computational performance of current algorithms available in these scientific libraries. The numerical experiments are conducted on the state-of-the-art HPC systems, and a single-core performance model is used to better characterize the efficiency of the solvers. Our studies indicate that the proposed non-negative computational framework for diffusion-type equations exhibits excellent strong scaling for real-world large-scale problems.Graphical AbstractThis figure shows the fate of chromium after 180 days using the single-field Galerkin formulation. The white regions indicate the violation of the non-negative constraint.

Material degradation due to moisture and temperature. Part 1: mathematical model, analysis, and analytical solutions

The mechanical response, serviceability, and load-bearing capacity of materials and structural components can be adversely affected due to external stimuli, which include exposure to a corrosive chemical species, high temperatures, temperature fluctuations (i.e., freezing–thawing), cyclic mechanical loading, just to name a few. It is, therefore, of paramount importance in several branches of engineering—ranging from aerospace engineering, civil engineering to biomedical engineering—to have a fundamental understanding of degradation of materials, as the materials in these applications are often subjected to adverse environments. As a result of recent advancements in material science, new materials such as fiber-reinforced polymers and multi-functional materials that exhibit high ductility have been developed and widely used, for example, as infrastructural materials or in medical devices (e.g., stents). The traditional small-strain approaches of modeling these materials will not be adequate. In this paper, we study degradation of materials due to an exposure to chemical species and temperature under large strain and large deformations. In the first part of our research work, we present a consistent mathematical model with firm thermodynamic underpinning. We then obtain semi-analytical solutions of several canonical problems to illustrate the nature of the quasi-static and unsteady behaviors of degrading hyperelastic solids.

Variational inequality approach to enforcing the non-negative constraint for advection–diffusion equations

Abstract Predictive simulations are crucial for the success of many subsurface applications, and it is highly desirable to obtain accurate non-negative solutions for transport equations in these numerical simulations. To this end, optimization-based methodologies based on quadratic programming (QP) have been shown to be a viable approach to ensuring discrete maximum principles and the non-negative constraint for anisotropic diffusion equations. In this paper, we propose a computational framework based on the variational inequality (VI) which can also be used to enforce important mathematical properties (e.g., maximum principles) and physical constraints (e.g., the non-negative constraint). We demonstrate that this framework is not only applicable to diffusion equations but also to non-symmetric advection–diffusion equations. An attractive feature of the proposed framework is that it works with any weak formulation for the advection–diffusion equations, including single-field formulations, which are computationally attractive. A particular emphasis is placed on the parallel and algorithmic performance of the VI approach across large-scale and heterogeneous problems. It is also shown that QP and VI are equivalent under certain conditions. State-of-the-art QP and VI solvers available from the PETSc library are used on a variety of steady-state 2D and 3D benchmarks, and a comparative study on the scalability between the QP and VI solvers is presented. We then extend the proposed framework to transient problems by simulating the miscible displacement of fluids in a heterogeneous porous medium and illustrate the importance of enforcing maximum principles for these types of coupled problems. Our numerical experiments indicate that VIs are indeed a viable approach for enforcing the maximum principles and the non-negative constraint in a large-scale computing environment. Also provided are Firedrake project files as well as a discussion on the computer implementation to help facilitate readers in understanding the proposed framework.

A Mixed Formulation for a Modification to Darcy Equation Based on Picard Linearization and Numerical Solutions to Large-Scale Realistic Problems

In this paper, we consider a modification to the Darcy equation by taking into account the dependence of viscosity on the pressure. We present a stabilized mixed formulation for the resulting governing equations. Equal-order interpolation for the velocity and pressure is considered, and shown to be stable (which is not the case under the classical mixed formulation). The proposed mixed formulation is tested using a wide variety of numerical examples. The proposed formulation is also implemented in a parallel setting, and the performance of the formulation for large-scale problems is illustrated using a representative problem. Two practical and technologically important problems, one each on enhanced oil recovery and geological carbon-dioxide sequestration, are solved using the proposed formulation. The numerical examples show that the predictions based on the Darcy model are qualitatively and quantitatively different from the predictions based on the modified Darcy model, which takes into account the dependence of the viscosity on the pressure. In particular, the numerical example on the geological carbon-dioxide sequestration shows that the Darcy model over-predicts the leakage into an abandoned well when compared to that of the modified Darcy model. On the other hand, the modified Darcy model predicts higher pressures and higher pressure gradients near the injection well. These predictions have dire consequences in predicting damage and fracture zones, and in designing the seal, whose integrity is crucial to the safety of a geological carbon-dioxide sequestration geosystem.

Modeling Flow in Porous Media with Double Porosity/Permeability: Mathematical Model, Properties, and Analytical Solutions

Geo-materials such as vuggy carbonates are known to exhibit multiple spatial scales. A common manifestation of spatial scales is the presence of (at least) two different scales of pores, which is commonly referred to as double porosity. To complicate things, the pore-network at each scale exhibits different permeability, and these networks are connected through fissure and conduits. Although some models are available in the literature, they lack a strong theoretical basis. This paper aims to fill this lacuna by providing the much needed theoretical foundations of the flow in porous media which exhibit double porosity/permeability. We first obtain a mathematical model for double porosity/permeability using the maximization of rate of dissipation hypothesis, and thereby providing a firm thermodynamic underpinning. We then present, along with mathematical proofs, several important mathematical properties that the solutions to the double porosity/permeability model satisfy. These properties are important in their own right as well as serve as good (mechanics-based) a posteriori measures to assess the accuracy of numerical solutions. We also present several canonical problems and obtain the corresponding analytical solutions, which are used to gain insights into the velocity and pressure profiles, and the mass transfer across the two pore-networks. In particular, we highlight how the solutions under the double porosity/permeability differ from the corresponding solutions under Darcy equations.