Super-Grid Modeling of the Elastic Wave Equation in Semi-Bounded Domains

We develop a super-grid modeling technique for solving the elastic wave equation in semi-bounded two- and three-dimensional spatial domains. In this method, waves are slowed down and dissipated in sponge layers near the far-field boundaries. Mathematically, this is equivalent to a coordinate mapping that transforms a very large physical domain to a significantly smaller computational domain, where the elastic wave equation is solved numerically on a regular grid. To damp out waves that become poorly resolved because of the coordinate mapping, a high order artificial dissipation operator is added in layers near the boundaries of the computational domain. We prove by energy estimates that the super-grid modeling leads to a stable numerical method with decreasing energy, which is valid for heterogeneous material properties and a free surface boundary condition on one side of the domain. Our spatial discretization is based on a fourth order accurate finite difference method, which satisfies the principle of summation by parts. We show that the discrete energy estimate holds also when a centered finite difference stencil is combined with homogeneous Dirichlet conditions at several ghost points outside of the far-field boundaries. Therefore, the coefficients in the finite difference stencils need only be boundary modified near the free surface. This allows for improved computational efficiency and significant simplifications of the implementation of the proposed method in multi-dimensional domains. Numerical experiments in three space dimensions show that the modeling error from truncating the domain can be made very small by choosing a sufficiently wide super-grid damping layer. The numerical accuracy is first evaluated against analytical solutions of Lamb's problem, where fourth order accuracy is observed with a sixth order artificial dissipation. We then use successive grid refinements to study the numerical accuracy in the more complicated motion due to a point moment tensor source in a regularized layered material.     

Two-Grid Method for Miscible Displacement Problem by Mixed Finite Element Methods and Mixed Finite Element Method of Characteristics

The miscible displacement of one incompressible fluid by another in a porous medium is governed by a system of two equations. One is elliptic form equation for the pressure and the other is parabolic form equation for the concentration of one of the fluids. Since only the velocity and not the pressure appears explicitly in the concentration equation, we use a mixed finite element method for the approximation of the pressure equation and mixed finite element method with characteristics for the concentration equation. To linearize the mixed-method equations, we use a two-grid algorithm based on the Newton iteration method for this full discrete scheme problems. First, we solve the original nonlinear equations on the coarse grid, then, we solve the linearized problem on the fine grid used Newton iteration once. It is shown that the coarse grid can be much coarser than the fine grid and achieve asymptotically optimal approximation as long as the mesh sizes satisfy $h=H^2$ in this paper. Finally, numerical experiment indicates that two-grid algorithm is very effective.     

Numerical Modeling of Anisotropic Elastic-Wave Sensitivity Propagation for Optimal Design of Time-Lapse Seismic Surveys

Reliable subsurface time-lapse seismic monitoring is crucial for many geophysical applications, such as enhanced geothermal system characterization, geologic carbon utilization and storage, and conventional and unconventional oil/gas reservoir characterization, etc. We develop an elastic-wave sensitivity propagation method for optimal design of cost-effective time-lapse seismic surveys considering the fact that most of subsurface geologic layers and fractured reservoirs are anisotropic instead of isotropic. For anisotropic media, we define monitoring criteria using qP- and qS-wave sensitivity energies after decomposing qP- and qS-wave components from the total elastic-wave sensitivity wavefield using a hybrid time- and frequency-domain approach. Geophones should therefore be placed at locations with significant qP- and qS-wave sensitivity energies for cost-effective time-lapse seismic monitoring in an anisotropic geology setting. Our numerical modeling results for a modified anisotropic Hess model demonstrate that, compared with the isotropic case, subsurface anisotropy changes the spatial distributions of elastic-wave sensitivity energies. Consequently, it is necessary to consider subsurface anisotropies when designing the spatial distribution of geophones for cost-effective time-lapse seismic monitoring. This finding suggests that it is essential to use our new anisotropic elastic-wave sensitivity modeling method for optimal design of time-lapse seismic surveys to reliably monitor the changes in subsurface reservoirs, fracture zones or target monitoring regions.     

Direct Simulation Methods for Scalar-Wave Envelopes in Two-Dimensional Layered Random Media Based on the Small-Angle Scattering Approximation

This study presents stochastic methods to simulate wave envelopes in layered random media. High-frequency Seismograms of small earthquakes are so complex due to lithospheric inhomogeneity that seismologists often analyze wave envelopes rather than wave traces to quantify the subsurface inhomogeneity. Since the statistical properties of the inhomogeneity vary regionally, it is important to develop and examine direct envelope simulation methods for non-uniform random media. As a simple example, this study supposes plane wave propagation through two-layer random media in 2-D composed of weak and strong inhomogeneity zones. The characteristic spatial-scale of the inhomogeneity is supposed to be larger than the wavelength, where small-angle scattering around the forward direction dominates large-angle scattering. Two envelope simulation methods based on the small-angle scattering approximation are examined. One method is to solve a differential equation for the two-frequency mutual coherence function with the Markov approximation. The other is to solve the stochastic ray bending process by using the Monte Carlo method based on the Markov approximation for the mutual coherence function. The resultant wave envelopes of the two methods showed excellent coincidence both for uniform and for two-layer random media. Furthermore, we confirmed the validity of the two methods comparing with the envelopes made from the finite difference simulations of waves. The two direct envelope simulation methods presented in this study can be a mathematical base for the study of high-frequency wave propagation through randomly inhomogeneous lithosphere in seismology.     

Multiscale Finite Element Methods for Flows on Rough Surfaces

In this paper, we present the Multiscale Finite Element Method (MsFEM) for problems on rough heterogeneous surfaces. We consider the diffusion equation on oscillatory surfaces. Our objective is to represent small-scale features of the solution via multiscale basis functions described on a coarse grid. This problem arises in many applications where processes occur on surfaces or thin layers. We present a unified multiscale finite element framework that entails the use of transformations that map the reference surface to the deformed surface. The main ingredients of MsFEM are (1) the construction of multiscale basis functions and (2) a global coupling of these basis functions. For the construction of multiscale basis functions, our approach uses the transformation of the reference surface to a deformed surface. On the deformed surface, multiscale basis functions are defined where reduced (1D) problems are solved along the edges of coarse-grid blocks to calculate nodal multiscale basis functions. Furthermore, these basis functions are transformed back to the reference configuration. We discuss the use of appropriate transformation operators that improve the accuracy of the method. The method has an optimal convergence if the transformed surface is smooth and the image of the coarse partition in the reference configuration forms a quasiuniform partition. In this paper, we consider such transformations based on harmonic coordinates (following H. Owhadi and L. Zhang [Comm. Pure and Applied Math., LX(2007), pp. 675–723]) and discuss gridding issues in the reference configuration. Numerical results are presented where we compare the MsFEM when two types of deformations are used for multiscale basis construction. The first deformation employs local information and the second deformation employs a global information. Our numerical results show that one can improve the accuracy of the simulations when a global information is used.     

Sequential Multiscale Modeling Using Sparse Representation

The main obstacle in sequential multiscale modeling is the pre-computation of the constitutive relation which often involves many independent variables. The constitutive relation of a polymeric fluid is a function of six variables, even after making the simplifying assumption that stress depends only on the rate of strain. Precomputing such a function is usually considered too expensive. Consequently the value of sequential multiscale modeling is often limited to "parameter passing". Here we demonstrate that sparse representations can be used to drastically reduce the computational cost for precomputing functions of many variables. This strategy dramatically increases the efficiency of sequential multiscale modeling, making it very competitive in many situations.     

A Wavelet-Adaptive Method for Multiscale Simulation of Turbulent Flows in Flying Insects

We present a wavelet-based adaptive method for computing 3D multiscale flows in complex, time-dependent geometries, implemented on massively parallel computers. While our focus is on simulations of flapping insects, it can be used for other flow problems. We model the incompressible fluid with an artificial compressibility approach in order to avoid solving elliptical problems. No-slip and in/outflow boundary conditions are imposed using volume penalization. The governing equations are discretized on a locally uniform Cartesian grid with centered finite differences, and integrated in time with a Runge–Kutta scheme, both of 4th order. The domain is partitioned into cubic blocks with different resolution and, for each block, biorthogonal interpolating wavelets are used as refinement indicators and prediction operators. Thresholding the wavelet coefficients allows to generate dynamically evolving grids, and an adaption strategy tracks the solution in both space and scale. Blocks are distributed among MPI processes and the grid topology is encoded using a tree-like data structure. Analyzing the different physical and numerical parameters allows us to balance their errors and thus ensures optimal convergence while minimizing computational effort. Different validation tests score accuracy and performance of our new open source code, WABBIT. Flow simulations of flapping insects demonstrate its applicability to complex, bio-inspired problems.     

A Fourth-Order Upwinding Embedded Boundary Method (UEBM) for Maxwell's Equations in Media with Material Interfaces: Part I

In this paper, we present a new fourth-order upwinding embedded boundary method (UEBM) over Cartesian grids, originally proposed in the Journal of Computational Physics [190 (2003), pp. 159-183.] as a second-order method for treating material interfaces for Maxwell's equations. In addition to the idea of the UEBM to evolve solutions at interfaces, we utilize the ghost fluid method to construct finite difference approximation of spatial derivatives at Cartesian grid points near the material interfaces. As a result, Runge-Kutta type time discretization can be used for the semidiscretized system to yield an overall fourth-order method, in contrast to the original second-order UEBM based on a Lax-Wendroff type difference. The final scheme allows time step sizes independent of the interface locations. Numerical examples are given to demonstrate the fourth-order accuracy as well as the stability of the method. We tested the scheme for several wave problems with various material interface locations, including electromagnetic scattering of a plane wave incident on a planar boundary and a two-dimensional electromagnetic application with an interface parallel to the y-axis.     

Multilevel Markov Chain Monte Carlo Method for High-Contrast Single-Phase Flow Problems

In this paper we propose a general framework for the uncertainty quantification of quantities of interest for high-contrast single-phase flow problems. It is based on the generalized multiscale finite element method (GMsFEM) and multilevel Monte Carlo (MLMC) methods. The former provides a hierarchy of approximations of different resolution, whereas the latter gives an efficient way to estimate quantities of interest using samples on different levels. The number of basis functions in the online GMsFEM stage can be varied to determine the solution resolution and the computational cost, and to efficiently generate samples at different levels. In particular, it is cheap to generate samples on coarse grids but with low resolution, and it is expensive to generate samples on fine grids with high accuracy. By suitably choosing the number of samples at different levels, one can leverage the expensive computation in larger fine-grid spaces toward smaller coarse-grid spaces, while retaining the accuracy of the final Monte Carlo estimate. Further, we describe a multilevel Markov chain Monte Carlo method, which sequentially screens the proposal with different levels of approximations and reduces the number of evaluations required on fine grids, while combining the samples at different levels to arrive at an accurate estimate. The framework seamlessly integrates the multiscale features of the GMsFEM with the multilevel feature of the MLMC methods following the work in [26], and our numerical experiments illustrate its efficiency and accuracy in comparison with standard Monte Carlo estimates.     

Mixed Multiscale Finite Volume Methods for Elliptic Problems in Two-Phase Flow Simulations

We develop a framework for constructing mixed multiscale finite volume methods for elliptic equations with multiple scales arising from flows in porous media. Some of the methods developed using the framework are already known [20]; others are new. New insight is gained for the known methods and extra flexibility is provided by the new methods. We give as an example a mixed MsFV on uniform mesh in 2-D. This method uses novel multiscale velocity basis functions that are suited for using global information, which is often needed to improve the accuracy of the multiscale simulations in the case of continuum scales with strong non-local features. The method efficiently captures the small effects on a coarse grid. We analyze the new mixed MsFV and apply it to solve two-phase flow equations in heterogeneous porous media. Numerical examples demonstrate the accuracy and efficiency of the proposed method for modeling the flows in porous media with non-separable and separable scales.     

3D hierarchical geometric modeling and multiscale FE analysis as a base for individualized medical diagnosis of bone structure

This paper describes a new alternative for individualized mechanical analysis of bone trabecular structure. This new method closes the gap between the classic homogenization approach that is applied to macro-scale models and the modern micro-finite element method that is applied directly to micro-scale high-resolution models. The method is based on multiresolution geometrical modeling that generates intermediate structural levels. A new method for estimating multiscale material properties has also been developed to facilitate reliable and efficient mechanical analysis. What makes this method unique is that it enables direct and interactive analysis of the model at every intermediate level. Such flexibility is of principal importance in the analysis of trabecular porous structure. The method enables physicians to zoom-in dynamically and focus on the volume of interest (VOI), thus paving the way for a large class of investigations into the mechanical behavior of bone structure. This is one of the very few methods in the field of computational bio-mechanics that applies mechanical analysis adaptively on large-scale high resolution models. The proposed computational multiscale FE method can serve as an infrastructure for a future comprehensive computerized system for diagnosis of bone structures. The aim of such a system is to assist physicians in diagnosis, prognosis, drug treatment simulation and monitoring. Such a system can provide a better understanding of the disease, and hence benefit patients by providing better and more individualized treatment and high quality healthcare. In this paper, we demonstrate the feasibility of our method on a high-resolution model of vertebra L3.     

An Efficient Two-Grid Scheme for the Cahn-Hilliard Equation

A two-grid method for solving the Cahn-Hilliard equation is proposed in this paper. This two-grid method consists of two steps. First, solve the Cahn-Hilliard equation with an implicit mixed finite element method on a coarse grid. Second, solve two Poisson equations using multigrid methods on a fine grid. This two-grid method can also be combined with local mesh refinement to further improve the efficiency. Numerical results including two and three dimensional cases with linear or quadratic elements show that this two-grid method can speed up the existing mixed finite method while keeping the same convergence rate.     

Stable and Efficient Modeling of Anelastic Attenuation in Seismic Wave Propagation

We develop a stable finite difference approximation of the three-dimensional viscoelastic wave equation. The material model is a super-imposition of N standard linear solid mechanisms, which commonly is used in seismology to model a material with constant quality factor Q. The proposed scheme discretizes the governing equations in second order displacement formulation using 3N memory variables, making it significantly more memory efficient than the commonly used first order velocity-stress formulation. The new scheme is a generalization of our energy conserving finite difference scheme for the elastic wave equation in second order formulation [SIAM J. Numer. Anal., 45 (2007), pp. 1902–1936]. Our main result is a proof that the proposed discretization is energy stable, even in the case of variable material properties. The proof relies on the summation-by-parts property of the discretization. The new scheme is implemented with grid refinement with hanging nodes on the interface. Numerical experiments verify the accuracy and stability of the new scheme. Semi-analytical solutions for a half-space problem and the LOH.3 layer over half-space problem are used to demonstrate how the number of viscoelastic mechanisms and the grid resolution influence the accuracy. We find that three standard linear solid mechanisms usually are sufficient to make the modeling error smaller than the discretization error.     

Numerical Investigation on the Boundary Conditions for the Multiscale Base Functions

We study the multiscale finite element method for solving multiscale elliptic problems with highly oscillating coefficients, which is designed to accurately capture the large scale behaviors of the solution without resolving the small scale characters. The key idea is to construct the multiscale base functions in the local partial differential equation with proper boundary conditions. The boundary conditions are chosen to extract more accurate boundary information in the local problem. We consider periodic and non-periodic coefficients with linear and oscillatory boundary conditions for the base functions. Numerical examples will be provided to demonstrate the effectiveness of the proposed multiscale finite element method.     

A Three-Level Multi-Continua Upscaling Method for Flow Problems in Fractured Porous Media

Traditional two level upscaling techniques suffer from a high offline cost when the coarse grid size is much larger than the fine grid size one. Thus, multilevel methods are desirable for problems with complex heterogeneities and high contrast. In this paper, we propose a novel three-level upscaling method for flow problems in fractured porous media. Our method starts with a fine grid discretization for the system involving fractured porous media. In the next step, based on the fine grid model, we construct a nonlocal multi-continua upscaling (NLMC) method using an intermediate grid. The system resulting from NLMC gives solutions that have physical meaning. In order to enhance locality, the grid size of the intermediate grid needs to be relatively small, and this motivates using such an intermediate grid. However, the resulting NLMC upscaled system has a relatively large dimension. This motivates a further step of dimension reduction. In particular, we will apply the idea of the Generalized Multiscale Finite Element Method (GMsFEM) to the NLMC system to obtain a final reduced model. We present simulation results for a two-dimensional model problem with a large number of fractures using the proposed three-level method.     

A Tailored Finite Point Method for Solving Steady MHD Duct Flow Problems with Boundary Layers

In this paper we propose a development of the finite difference method, called the tailored finite point method, for solving steady magnetohydrodynamic (MHD) duct flow problems with a high Hartmann number. When the Hartmann number is large, the MHD duct flow is convection-dominated and thus its solution may exhibit localized phenomena such as the boundary layer. Most conventional numerical methods can not efficiently solve the layer problem because they are lacking in either stability or accuracy. However, the proposed tailored finite point method is capable of resolving high gradients near the layer regions without refining the mesh. Firstly, we devise the tailored finite point method for the scalar inhomogeneous convection-diffusion problem, and then extend it to the MHD duct flow which consists of a coupled system of convection-diffusion equations. For each interior grid point of a given rectangular mesh, we construct a finite-point difference operator at that point with some nearby grid points, where the coefficients of the difference operator are tailored to some particular properties of the problem. Numerical examples are provided to show the high performance of the proposed method.     

Numerical Methods for Two-Fluid Dispersive Fast MHD Phenomena

The finite volume wave propagation method and the finite element RungeKutta discontinuous Galerkin (RKDG) method are studied for applications to balance laws describing plasma fluids. The plasma fluid equations explored are dispersive and not dissipative. The physical dispersion introduced through the source terms leads to the wide variety of plasma waves. The dispersive nature of the plasma fluid equations explored separates the work in this paper from previous publications. The linearized Euler equations with dispersive source terms are used as a model equation system to compare the wave propagation and RKDG methods. The numerical methods are then studied for applications of the full two-fluid plasma equations. The two-fluid equations describe the self-consistent evolution of electron and ion fluids in the presence of electromagnetic fields. It is found that the wave propagation method, when run at a CFL number of 1, is more accurate for equation systems that do not have disparate characteristic speeds. However, if the oscillation frequency is large compared to the frequency of information propagation, source splitting in the wave propagation method may cause phase errors. The Runge-Kutta discontinuous Galerkin method provides more accurate results for problems near steady-state as well as problems with disparate characteristic speeds when using higher spatial orders.     

A Simulation Approach Including Underresolved Scales for Two-Component Fluid Flows in Multiscale Porous Structures

In this study, we develop computational models and a methodology for accurate multicomponent flow simulation in underresolved multiscale porous structures [1]. It is generally impractical to fully resolve the flow in porous structures with large length-scale differences due to the tremendously high computational expense. The flow contributions from underresolved scales should be taken into account with proper physics modeling and simulation processes. Using precomputed physical properties such as the absolute permeability, $K_0,$ the capillary pressure-saturation curve, and the relative permeability, $K_r,$ in typical resolved porous structures, the local fluid force is determined and applied to simulations in the underresolved regions, which are represented by porous media. In this way, accurate flow simulations in multiscale porous structures become feasible.To evaluate the accuracy and robustness of this method, a set of benchmark test cases are simulated for both single-component and two-component flows in artificially constructed multiscale porous structures. Using comparisons with analytic solutions and results with much finer resolution resolving the porous structures, the simulated results are examined. Indeed, in all cases, the results successfully show high accuracy with proper input of $K_0,$ capillary pressure, and $K_r.$ Specifically, imbibition patterns, entry pressure, residual component patterns, and absolute/relative permeability are accurately captured with this approach.     

Multiscale Modeling and Simulations of Flows in Naturally Fractured Karst Reservoirs

Modeling and numerical simulations of fractured, vuggy, porus media is a challenging problem which occurs frequently in reservoir engineering. The problem is especially relevant in flow simulations of karst reservoirs where vugs and caves are embedded in a porous rock and are connected via fracture networks at multiple scales. In this paper we propose a unified approach to this problem by using the Stokes-Brinkman equations at the fine scale. These equations are capable of representing porous media such as rock as well as free flow regions (fractures, vugs, caves) in a single system of equations. We then consider upscaling these equations to a coarser scale. The cell problems, needed to compute coarse-scale permeability of Representative Element of Volume (REV) are discussed. A mixed finite element method is then used to solve the Stokes-Brinkman equation at the fine scale for a number of flow problems, representative for different types of vuggy reservoirs. Upscaling is also performed by numerical solutions of Stokes-Brinkman cell problems in selected REVs. Both isolated vugs in porous matrix as well as vugs connected by fracture networks are analyzed by comparing fine-scale and coarse-scale flow fields. Several different types of fracture networks, representative of short- and long-range fractures are studied numerically. It is also shown that the Stokes-Brinkman equations can naturally be used to model additional physical effects pertaining to vugular media such as partial fracture with fill-in by some material and/or fluids with suspended solid particles.     

Effective Two-Level Domain Decomposition Preconditioners for Elastic Crack Problems Modeled by Extended Finite Element Method

In this paper, we propose some effective one- and two-level domain decomposition preconditioners for elastic crack problems modeled by extended finite element method. To construct the preconditioners, the physical domain is decomposed into the "crack tip" subdomain, which contains all the degrees of freedom (dofs) of the branch enrichment functions, and the "regular" subdomains, which contain the standard dofs and the dofs of the Heaviside enrichment function. In the one-level additive Schwarz and restricted additive Schwarz preconditioners, the "crack tip" subproblem is solved directly and the "regular" subproblems are solved by some inexact solvers, such as ILU. In the two-level domain decomposition preconditioners, traditional interpolations between the coarse and the fine meshes destroy the good convergence property. Therefore, we propose an unconventional approach in which the coarse mesh is exactly the same as the fine mesh along the crack line, and adopt the technique of a non-matching grid interpolation between the fine and the coarse meshes. Numerical experiments demonstrate the effectiveness of the two-level domain decomposition preconditioners applied to elastic crack problems.