Implementation of a Locally Conservative Numerical Subgrid Upscaling Scheme for Two-Phase Darcy Flow

A method for performing nested grid calculations with a large-eddy simulation code is described. A common numerical method is used for all meshes, and the grid architecture consists of a single outer or coarse grid, and nested or fine grids, which overlap in some common region. Inter-grid communication matches the velocity, pressure and potential temperature fields in the overlap region. Resolved and sub-grid scale (SGS) turbulent fluxes and kinetic energy on the fine grid are averaged to the coarse grid using a conservation rule equivalent to Germano's identity used to develop dynamic SGS models. Simulations of a slightly convective, strong shear planetary boundary layer were carried out with varying surface-layer resolutions. Grid refinements in the (x, y, z) directions of up to (5, 5, 2) times were employed. Two-way interaction solutions on the coarse and fine meshes are successfully matched in the overlap region on an instantaneous basis, and the turbulent motions on the fine grid blend smoothly into the coarse grid across the grid interface. With surface-layer grid nesting, significant increases in resolved eddy fluxes and variances are found. The energy-scale content of the vertical velocity, and hence vertical turbulent fluxes, appear to be most influenced by increased grid resolution. Vertical velocity spectra show that the dominant scale shifts towards higher wavenumbers (smaller scales) and the magnitude of the peak energy is increased by more than a factor of 3 with finer resolution. Outside of the nested region the average heat and momentum fluxes and spectra are slightly influenced by the fine resolution in the surface layer. From these results we conclude that fine resolution is required to resolve the details of the turbulent motions in the surface layer. At the same time, however, increased resolution in the surface layer does not appreciably alter the ensemble statistics of the resolved and SGS motions outside of the nested region.

An error arises in large eddy simulations when the order of spatial differentiation and filtering operations is commuted on a non-uniform grid. The magnitude of the commutation error is controlled by spatial variations in the width of the filter function, which is conventionally defined in terms of the local grid size. We have previously found large commutation errors in the sub-grid scale (SGS) terms of a numerical dynamo model, particularly near boundaries where the grid size changes most rapidly. In the present study, we propose a correction for the commutation error which is specifically designed for use with the nonlinear gradient SGS model. The commutation error correction for each SGS term can be expressed as the product of the spatial derivative of the second-order moment of the filter function and the second spatial derivative of the SGS term. We test the correction using output from a fully resolved convection-driven dynamo simulation in a rotating plane layer. As an example, we evaluate the SGS heat flux on a coarser grid with and without the commutation error correction. The result is tested using the resolved dynamo solution on a finer grid. We find that the commutation error correction is comparable in magnitude to the SGS heat flux near the boundary. Addition of the correction term significantly improves the modeled SGS heat flux near the boundaries.

Flow fields in a centrifugal fan for an indoor unit of an air-conditioner were calculated with finite element method-based large eddy simulation (LES) with the aim of predicting fan performance and aerodynamic noise in this study. The numerical simulation code employed throughout the LES was called FrontFlow/blue (FFB). We compared 10M grid [coarse grid] and 60M grid [fine grid] calculation results for investigation of influence of grid resolution. In the fine grid, the number of grid elements in blade-to-blade direction, and of region between the shroud and the bell mouth increased in particular. By calculating with the fine grid, calculated distributions of absolute velocities at blade exit reasonably agreed with experimental results. Because of this, maximum absolute velocity by fine grid near hub decreased as compared to those by coarse grid. Calculated sound pressure level by fine grid was therefore smaller than that by coarse grid, and the overestimation of sound pressure was suppressed by calculating with fine grid. This decrease of the absolute velocity was a first factor for the improvement of calculation accuracy. Moreover, number of captured streaks on the blade, hub, and shroud surfaces by fine grid increased as compared to those by coarse grid. As a result, size of streak by fine grid became smaller than that by coarse grid. Static pressure fluctuations by fine grid on the blade, hub, and shroud surfaces therefore reduced as compared to those by coarse grid. Aerodynamic noise was related to static pressure fluctuations according to Curle’s equation. This reduction of static pressure fluctuations was therefore a second factor for improvement of calculation accuracy.

Adaptive gridding methods are of fundamental importance both for industry and academia. As one of the computing methods, the Boundary Element Method (BEM) is used to simulate problems whose fundamental solutions are available. The method is usually characterised as constant elements BEM or linear elements BEM depending on the type of interpolation used at the elements. Its popularity is steadily growing because of its advantages over other numerical methods most important of which being that it only involves obtaining data at the boundary and computation of the solution in the domain is merely a case of post processing through the use of an identity. This results in the reduction of the problem dimension by unity. Although there is a reduction in dimension when we use BEM, the method usually results in full matrices which can be expensive to solve. This makes the method costly for problems that require very fine grids like those with hot spots as in the example of impressed current cathodic protection systems (ICCP). The BEM is a global method in nature in that the solution in one node depends on the solutions in all the other nodes of the grid. Hence an error in one node can pollute the solution in all the other nodes. In this thesis, we first focus our attention on defining and studying both the local and global errors for the BEM. This is not a completely new study as the literature suggests, however, our approach is different. We use the basic foundations of the the method to define the errors. Since the method is a global method, first we use the interpolation error on each element to define what we have called a sublocal error. Then using the sublocal error we have defined the local error. Understanding the local errors enabled us study the global error. Theoretical and numerical results show that these errors are second order in grid size for both the constant and linear element cases. Then, having explored errors, we study a method for adaptive grid refinement for the BEM. Rather than using a truly nonuniform grid, we present a method called local defect correction (LDC) that is based on local uniform grid refinement. This method is already developed and documented for other numerical methods such as finite difference and finite volume methods but not for BEM. In the LDC method, the discretization on a composite grid is based on a combination of standard discretizations on several uniform grids with different grid sizes that cover different parts of the domain. At least one grid, the uniform global coarse grid, should cover the entire boundary. The size of the global coarse grid is chosen in agreement with the relatively smooth behaviour of the solution outside the hot spots. Then several uniform local fine grids each of which covers only a (small) part of the boundary are used in the hot spots. The grid sizes of the local grids are chosen in agreement with the behaviour of the continuous solution in that part of the boundary. The LDC method is an iterative process whereby a basic global discretisation is improved by local discretisations defined in subdomains. The update of the coarse grid solution is achieved by adding a defect correction term to the right hand side of the coarse grid problem. At each iteration step, the process yields a discrete approximation of the continuous solution on the corresponding composite grid. We have shown how this discretisation can be achieved for the BEM. We apply the discretisation to an academic example to demonstrate its implementation and later show how to use it for an application as the ICCP system. The results show that it is a cheaper method than solving on a truly composite grid and converges in a single step.

Running a fine grid model with 107 - 109 of cells is possible using a supercomputer with 103 - 106 of CPUs but may not be always cost-effective. The most cost-effective way is to use a coarse grid model that is much smaller but with static/dynamic profiles very close to the fine grid model. This paper proposes a new layer optimization and upscaling method with the aim for creating a consistent coarse grid model. Unlike the industry's existing layer optimization and upscaling methods, the proposed method performs layer optimization and upscaling fully integrated with the Lorenz coefficient and curves (LCC). Coarse grid layers and their permeabilities are created by minimizing the difference between fine and coarse grid LCCs. The process consists of static and dynamic optimizations. The former is measured by LCC while the latter by pressure, GOR, and water-cut. A new LCC-based permeability upscaling method is developed to preserve the fine grid multiphase flow behaviors. A satisfactory coarse grid model is achieved when both static and dynamic criteria are met. The proposed method has been successfully applied to a giant carbonate oil field in the Caspian Sea that consists of a matrix dominated platform and a fracture/karst dominated rim. Due to the field's complex geology and high H2S content (15%), a dual porosity, dual permeability compositional model has been created to model compositional sour crude flow within and between the matrix and fracture/karst features. The reservoir drive mechanisms are fluid expansion, miscible gas injection and aquifer drive. The reservoir is undersaturated and has an abnormally high initial reservoir pressure. The fine-grid static model contains 104 million cells (370×225×625×2) and the optimized upscaled coarse-grid dynamic model has 8.3 million cells (370×225×50×2). The upscaled model can be run efficiently on the company's existing HPC infrastructure with a maximum of 64 CPUs. Excellent matches of the Lorenz coefficient maps for reservoir total/zones and Lorenz curves at all wells between the fine and coarse grid models have been achieved. Matches on the dynamic variables, e.g., pressure, gas breakthrough time, and GOR growth, in all producers are within the defined acceptable tolerances. The high quality of the static and dynamic matches between the coarse- and fine-grid models confirms that the reservoir properties of the coarse-grid model is very close to the fine-grid model and can be used a base model for history matching and uncertainty analysis.

Seismic data are routinely used to infer in situ properties of earth materials on many scales, ranging from global studies to investigations of surficial geological formations. While inversion and imaging algorithms utilizing these data have improved steadily, there are remaining challenges that make detailed measurements of the properties of some geologic materials very difficult. For example, the determination of the concentration and orientation of fracture systems is prohibitively expensive to simulate on the fine grid and, thus, some type of coarse-grid simulations are needed. In this paper, we describe a new multiscale finite element algorithm for simulating seismic wave propagation in heterogeneous media. This method solves the wave equation on a coarse grid using multiscale basis functions and a global coupling mechanism to relate information between fine and coarse grids. Using a mixed formulation of the wave equation and staggered discontinuous basis functions, the proposed multiscale methods have the following properties. • The total wave energy is conserved. • Mass matrix is diagonal on a coarse grid and explicit energy-preserving time discretization does not require solving a linear system at each time step. • Multiscale basis functions can accurately capture the subgrid variations of the solution and the time stepping is performed on a coarse grid. We discuss various subgrid capturing mechanisms and present some preliminary numerical results.

Permeability and porosity parameters of a porous medium are known only in a statistical sense. For risk assessment, one must perform multiple flow simulations of a single site, varying these input parameters. Because multiple simulations of large sites are computationally prohibitive, upscaling from fine to coarse scales is necessary. Traditional upscaling techniques determine a new effective or upscaled permeability field defined on a coarser scale, which is then used in a standard coarse grid discretization operator. We develop here a method of determining a new coarse grid discretization operator that provides an upscaled solution but bypasses the determination of effective permeability and porosity fields. The method has two steps. We first solve for fine scale flow information internal to each coarse grid cell. Because the problems are small, this step is relatively fast. Then we determine a modified coarse grid operator for solving the upscaled problem that includes the fine scale flow information from the first step. The method is developed for single-phase flow in the context of the mixed finite element method; therefore, the method is locally mass conservative. Unlike traditional upscaling methods (such as homogenization) we do not impose arbitrary boundary conditions on the coarse grid. We present comparisons of our method with the harmonic average permeability upscaling technique.

This paper presents a multidomain (MD) computational approach to multiscale (MS) compositional reservoir simulation amenable to parallel processing. The main objective is to minimize the upscaling process and preserve the pertinent flow characteristics of the fine-scale geological features of the reservoir. The preservation of the fine-scale reservoir characteristics yields higher-accuracy solution compared to an upscaled model. The computing algorithm involves breaking down the fine-grid (FG) computing mesh into an underlying coarse-grid (CG) subdomain. Unlike other MS techniques, the CG subdomain only serves as an assemblage of the FG cells and not as a computational CG overlaying the FG cells. The flow equations are then restructured into a global pressure equation which is solved on the FG using a large time step without solving for the compositions, saturations, etc. The convergence tolerance for the global pressure solution is soft (no iterations) to reduce computing time. The pressure solution is then used to calculate the molar flux of components only at the FG boundaries of the CG subdomain cells. Finally, each CG subdomain subject to its boundary conditions is solved independently using a smaller time step to compute cell pressures, compositions, and phase saturations. This procedure is amenable to parallel processing. This MS scheme was applied to time domain, in which the global pressure solution is obtained with a larger time step than in the local FG solution. This means that the molar flux boundary condition at the CG subdomain boundaries is not updated every time step, thus, buying additional computational advantage. This new technique was tested on several volatile oil and gas-condensate systems and excellent agreement was obtained with FG simulation with about 60% reduction in computing time using a single processor. Much greater time reduction is expected in a multi processing environment.

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 = H2 in this paper. Finally, numerical experiment indicates that two-grid algorithm is very effective.

A multiscale algorithm for a multiphase filtration problem is proposed. Filtration fluxes on a fine grid are determined from the solution of the pressure equation on a coarse grid. Further, the domain is decomposed into subdomains with an acceptable number of cells and the full second boundary condition filtration problem is solved using the fluxes. The support operator method has been improved for a complex structure cell for solution of the pressure equation on a coarse grid. This method is a high resolution one: a divergence operator has the approximation of the second order, while fluxes have the approximation of the first order. At the same time, the method allows revealing the solution’s properties related to a fine grid structure.

Two-grid mixed finite element method for nonlinear hyperbolic equations

Convergence analysis of two-grid fixed stress split iterative scheme for coupled flow and deformation in heterogeneous poroelastic media

This paper studies and analyzes a model describing the flow of contaminated brines through the porous media under severe thermal conditions caused by the radioactive contaminants. The problem is approximated based on combining the mixed finite element method with the modified method of characteristics. In order to solve the resulting algebraic nonlinear equations efficiently, a two-grid method is presented and discussed in this paper. This approach includes a small nonlinear system on a coarse grid with size H and a linear system on a fine grid with size h . It follows from error estimates that asymptotically optimal accuracy can be obtained as long as the mesh sizes satisfy H = O h 1 / 3 .

: To illustrate Brandt's multigrid method this report considers the Helmholtz equation in cases where the frequency is close to an eigenfrequency. In a second part the eigenvalue problem is discussed. The basic idea for the treatment of the inhomogeneous problem is contained in Brandt's publications. Difficulties are encountered because the coarse grid eigenvalues differ from those of the finer grid. As a consequence, approximations to a solution in the coarse grid may differ widely from the corresponding approximation in a finer grid, if the frequency defined by the inhomogeneous problem is close to one of the eigenfrequencies. This may be detrimental to the aim of the coarse grid step; namely to remove from an approximation in a fine grid long waves which cause divergence of fine grid iterations. This is remedied by eliminating long wave contributions to the residual directly in the next finer grid. The success of this procedure still depends upon the choice of the subspace from which the long wave corrections in the finer grid are taken. Even this modified method fails if the frequency pertaining to the inhomogeneous term is too close to an eigenvalue, unless the subspace contains the pertinent eigenfunction.

A three dimensional adaptive multiscale method for crack growth in Silicon