Implicit Discretization Research Articles

A numerical simulation model for multi-scale flow in tight oil reservoirs

Abstract A discrete fracture model for multi-scale flow in large-scale fractured tight oil reservoirs is proposed considering the compressibility of reservoir rock and fluid, and the non-linear flow in the tight matrix. Validation of the model is performed, followed by the field application of the model. The two-point flux-approximation scheme is adopted in the model to calculate conductivity, and small grids at the fracture intersections are eliminated by the “star-delta” transformation method to improve the computational stability. The fully implicit discretization scheme is performed on the temporal domain. Automatic differentiation technique which can improve model establishment efficiency and computational accuracy is applied in the model to solve the numerical model. The model is validated with the simulation results of Eclipse and the historical production data of a long fractured horizontal well in a tight oil reservoir in Xinjiang oilfield. Simulation results of a field-scale reservoir show that the model proposed can simulate reservoirs with large-scale complex fracture systems; well productivity is positively correlated with the scale of the stimulated reservoir volume, and the difference in planar fracture density and fracture connectivity are proved to be the key factors that lead to the heterogeneous distribution of remaining oil in tight oil reservoirs.

An energy stable monolithic Eulerian fluid‐structure finite element method

SummaryWhen written in an Eulerian frame, the conservation laws of continuum mechanics are similar for fluids and solids leading to a single set of variables for a monolithic formulation. Such formulations are well adapted to large displacement fluid‐structure configurations, but stability is a challenging problem because of moving geometries. In this article, the method is presented; time implicit discretizations are proposed with iterative algorithms well posed at each step, at least for small displacements; stability is discussed for an implicit in time finite element method in space by showing that energy decreases with time. The key numerical ingredient is the Characterics‐Galerkin method coupled with a powerful mesh generator. A numerical section discusses implementation issues and presents a few simple tests. It is also shown that contacts are easily handled by extending the method to variational inequalities. This paper deals only with incompressible neo‐Hookean Mooney‐Rivlin hyperelastic material in 2 dimensions in a Newtonian fluid, but the method is not limited to these; compressible and 3D cases will be presented later.

An efficient and stable hydrodynamic model with novel source term discretization schemes for overland flow and flood simulations

AbstractNumerical models solving the full 2‐D shallow water equations (SWEs) have been increasingly used to simulate overland flows and better understand the transient flow dynamics of flash floods in a catchment. However, there still exist key challenges that have not yet been resolved for the development of fully dynamic overland flow models, related to (1) the difficulty of maintaining numerical stability and accuracy in the limit of disappearing water depth and (2) inaccurate estimation of velocities and discharges on slopes as a result of strong nonlinearity of friction terms. This paper aims to tackle these key research challenges and present a new numerical scheme for accurately and efficiently modeling large‐scale transient overland flows over complex terrains. The proposed scheme features a novel surface reconstruction method (SRM) to correctly compute slope source terms and maintain numerical stability at small water depth, and a new implicit discretization method to handle the highly nonlinear friction terms. The resulting shallow water overland flow model is first validated against analytical and experimental test cases and then applied to simulate a hypothetic rainfall event in the 42 km2 Haltwhistle Burn, UK.

Wave Equation-Based Implicit Subdomain DGTD Method for Modeling of Electrically Small Problems

A second-order wave equation-based implicit discontinuous Galerkin time-domain (DGTD) method is proposed to efficiently model electrically small problems. The proposed method employs the second-order wave equation for electric field (or magnetic field) as the governing equation of the DG formulation, instead of the first-order Maxwell’s curl equations. A modified version of the Riemann solver (upwind flux) is introduced to evaluate the numerical flux resulting from the weak form of the wave equation. Compared with previous first-order Maxwell’s curl equation-based implicit DGTD methods, which typically solve all electric and magnetic field unknowns for each subdomain, the proposed method only needs to solve for the electric field unknowns plus the surface magnetic field unknowns at subdomain interfaces. This reduces the dimensions of the resultant linear system and thus allows for modeling larger problems. Furthermore, unlike element-based DGTD methods, the proposed method is subdomain-based. The computational region is divided into multiple subdomains based on the domain-decomposition method, and each subdomain may contain multiple elements. Different element types and orders of basis functions can be employed in different subdomains to exploit the geometry property of the model. A nonconformal mesh is allowed between different subdomains to increase meshing flexibility. The Newmark-beta time-integration scheme is used for implicit temporal discretization, and fast direct linear solvers, such as the lower-diagonal-upper decomposition algorithm, are employed to accelerate time integration when all the subdomains are in a sequential order. Numerical results show that the proposed method is more efficient in terms of CPU time, and also saves memory with respect to the previous implicit DGTD method when modeling electrically small problems.

Numerical modeling of non-Newtonian fluid flow in fractures and porous media

Non-Newtonian fluids having Bingham or power-law rheology are common in many applications within drilling and reservoir engineering. Examples of such fluids are drilling muds, foams, heavy oil, hydraulic-fracturing and other stimulation fluids, and cement slurries. Despite the importance of non-Newtonian rheology, it is rarely used in reservoir simulators and fracture flow simulations. We study two types of non-Newtonian rheology: the truncated power-law (Ostwald-de Waele) fluid and the Bingham fluid. For either of the two types of non-Newtonian rheology, we construct relationships between the superficial fluid velocity and the pressure gradient in fractures and porous media. The Bingham fluid is regularized by means of Papanastasiou-type regularization for porous media and by means of a simple hyperbolic function for fracture flow. Approximation by Taylor expansion is used to evaluate the fluid velocity for small pressure gradients to reduce rounding errors. We report simulations of flow in rough-walled fractures for different rheologies and study the effect of fluid parameters on the flow channelization in rough-walled fractures. This effect is known from previous studies. We demonstrate how rheologies on different domains can be included in a fully-unstructured reservoir simulation that incorporates discrete fracture modeling (DFM). The above formulation was implemented in the open-source MATLAB Reservoir Simulation Toolbox (MRST), which uses fully implicit discretization on general polyhedral grids, including industry standard grids with DFM. This robust implementation is an important step towards hydro-mechanically coupled simulation of hydraulic fracturing with realistic non-Newtonian fluid rheology since most hydraulic fracturing models implemented so far make use of oversimplified rheological models (e.g., Newtonian or pure power-law).

A matrix-free, implicit finite volume lattice Boltzmann method for steady flows

Abstract In the present paper, a matrix-free, implicit finite volume lattice Boltzmann method for steady flow on unstructured mesh is proposed. The approximate linear system arising from the implicit finite volume discretization of lattice Boltzmann equation (LBE) is solved by a novel algorithm, which combines the lower-upper symmetric Gauss–Seidel (LU-SGS) and the Jacobi iteration schemes. A remarkable feature of the present implicit method is that the storage of the Jacobian matrixes of the convection and collision terms can be completely eliminated by approximating the Jacobian matrix-solution incremental vector product with appropriate numerical flux incremental vector and the numerical increment of collision term, resulting a matrix-free implicit scheme. The present method is validated by several two-dimensional testing cases.

A Parallel Finite Element Method for 3D Two-Phase Moving Contact Line Problems in Complex Domains

Moving contact line problem plays an important role in fluid-fluid interface motion on solid surfaces. The problem can be described by a phase-field model consisting of the coupled Cahn–Hilliard and Navier–Stokes equations with the generalized Navier boundary condition (GNBC). Accurate simulation of the interface and contact line motion requires very fine meshes, and the computation in 3D is even more challenging. Thus, the use of high performance computers and scalable parallel algorithms are indispensable. In this paper, we generalize the GNBC to surfaces with complex geometry and introduce a finite element method on unstructured 3D meshes with a semi-implicit time integration scheme. A highly parallel solution strategy using different solvers for different components of the discretization is presented. More precisely, we apply a restricted additive Schwarz preconditioned GMRES method to solve the systems arising from implicit discretization of the Cahn–Hilliard equation and the velocity equation, and an algebraic multigrid preconditioned CG method to solve the pressure Poisson system. Numerical experiments show that the strategy is efficient and scalable for 3D problems with complex geometry and on a supercomputer with a large number of processors.

Finite-difference modeling à la Mickens of the distribution of the stopping time in a stochastic differential equation

Departing from a general stochastic differential equation with Brownian diffusion, we establish that the distribution of probability of the stopping time is governed by a parabolic partial differential equation. A particular form of the problem under investigation may be associated to a stochastic generalization of the well-known Paris’ law from structural mechanics, in which case, the solution of the boundary-value problem represents the probability distribution of the hitting time. An implicit, convergent and probability-based discretization to approximate the solution of the boundary-value problem is proposed in this work. Using a convenient vector representation of our scheme, we prove that the method preserves the most relevant properties of a probability distribution function, namely, the non-negativity, the boundedness from above by 1, and the monotonicity. In addition, we establish that our method is a convergent technique, and provide some illustrative comparisons against known exact solutions.

Fatigue life estimates under non-proportional loading through continuum damage evolution law

This paper seeks to use Lemaitre’s Continuum Damage Model to obtain fatigue life estimates under multiaxial proportional and non-proportional loadings. Initially, the Lemaitre’s mathematical model is presented, by formulating the necessary constitutive relations. Chaboche’s model was chosen to describe the kinematic hardening law. Then, the numerical model necessary to solve the constitutive relations is developed, utilizing the Chaboche’s law with 3 terms and Euler’s implicit discretization. Then, the material parameters are identified for 304 and S460N steels and 6061-T6 aluminum alloy. The model is implemented in a FORTRAN routine, which is submitted to uniaxial and proportional and non-proportional multiaxial loading histories. The fatigue life data obtained from Lemaitre’s damage model is compared to experimental data. Then, a stress amplitude analysis is conducted and the numerical stresses are compared to experimental data. Damage evolution curves are also obtained for each material and loading. The results show that Lemaitre’s damage model describes adequately the behavior of the analyzed materials under low cycle fatigue, when low strain amplitudes are being applied to the specimen.

Numerical studies of a class of linear solvers for fine-scale petroleum reservoir simulation

Numerical simulation based on fine-scale reservoir models helps petroleum engineers in understanding fluid flow in porous media and achieving higher recovery ratio. Fine-scale models give rise to large-scale linear systems, and thus require effective solvers for solving these linear systems to finish simulation in reasonable turn-around time. In this paper, we study convergence, robustness, and efficiency of a class of multi-stage preconditioners accelerated by Krylov subspace methods for solving Jacobian systems from a fully implicit discretization. We compare components of these preconditioners, including decoupling and sub-problem solvers, for fine-scale reservoir simulation. Several benchmark and real-world problems, including a ten-million-cell reservoir problem, were simulated on a desktop computer. Numerical tests show that the combination of the alternating block factorization method and multi-stage subspace correction preconditioner gives a robust and memory-efficient solver for fine-scale reservoir simulation.

Improving the convergence behaviour of a fixed‐point‐iteration solver for multiphase flow in porous media

SummaryA new method to admit large Courant numbers in the numerical simulation of multiphase flow is presented. The governing equations are discretized in time using an adaptive θ‐method. However, the use of implicit discretizations does not guarantee convergence of the nonlinear solver for large Courant numbers. In this work, a double‐fixed point iteration method with backtracking is presented, which improves both convergence and convergence rate. Moreover, acceleration techniques are presented to yield a more robust nonlinear solver with increased effective convergence rate. The new method reduces the computational effort by strengthening the coupling between saturation and velocity, obtaining an efficient backtracking parameter, using a modified version of Anderson's acceleration and adding vanishing artificial diffusion. © 2016 The Authors. International Journal for Numerical Methods in Fluids Published by John Wiley & Sons Ltd.

A local level set method based on a finite element method for unstructured meshes

A local level set method for unstructured meshes has been implemented by using a finite element method. A least-square weighted residual method was employed for implicit discretization to solve the level set advection equation. By contrast, a direct re-initialization method, which is directly applicable to the local level set method for unstructured meshes, was adopted to re-correct the level set function to become a signed distance function after advection. The proposed algorithm was constructed such that the advection and direct reinitialization steps were conducted only for nodes inside the narrow band around the interface. Therefore, in the advection step, the Gauss–Seidel method was used to update the level set function using a node-by-node solution method. Some benchmark problems were solved by using the present local level set method. Numerical results have shown that the proposed algorithm is accurate and efficient in terms of computational time.

A deterministic model for the distribution of the stopping time in a stochastic equation and its numerical solution

In this work, we consider a stochastic differential equation that generalizes the well known Paris’ equation from fracture of materials. The model describes the propagation of cracks on solids, and it includes a deterministic summand and a stochastic component in terms of a Brownian motion. The use of Itô’s stochastic integral gives an equivalent stochastic integral equation that is further generalized here. We note that the probability distribution of the stopping time of the general model satisfies a deterministic diffusion–advection partial differential equation for which the solution is known only in a reduced number of particular cases. Motivated by these analytical results, we develop a fast finite-difference method to approximate the distribution of the stopping time. The method is an explicit exponential-like technique that preserves the main features of a probability distribution, namely, the non-negativity, the boundedness from above by 1 as well as the spatial monotonicity. Moreover, the method is a monotone technique that is also capable of preserving the temporal monotonicity of the approximations. These properties of the proposed methodology are thoroughly established in the present manuscript. A continuity condition of the numerical solutions in terms of the initial conditions and the temporal computational parameter is established also, together with a limiting property of the methodology when the free parameter tends to infinity. For comparison purposes, we are providing an implicit and stable discretization of the mathematical model which has a second order of convergence but for which conditions that guarantee the positivity, the boundedness and the monotonicity of approximations are not available. The numerical simulations obtained with implementations of our techniques show that the explicit method is an efficient scheme that preserves the characteristics of interest (non-negativity, boundedness from above by 1 and monotonicity), and that the numerical approximations are in good agreement with the known exact solutions.

A dynamic adaptation method based on unstructured mesh for solving sloshing problems

A dynamic mesh-adaptation algorithm based on an unstructured mesh was used to solve sloshing problems by a finite element method. The free surface evolution of sloshing with a low filling ratio is so violent that a local mesh coarsening/refinement (LCR) technique can be effectively used to resolve the flow field near the interface of two immiscible fluids. Since an implicit discretization method was employed to solve the incompressible Navier-Stokes equations, an assembled global matrix was generated using a dynamic compressed sparse row (CSR) method. The proposed algorithm was validated by comparing the simulation results with those of a non-adaptive method with respect to wall impact pressures and mass conservation. Numerical results showed that relative mass error strongly depends on mesh resolution near the interface. Results of adaptive simulations were found to be comparable with those of non-adaptive simulations only if a similar mesh resolution was used near the interface. Moreover, adaptive simulations were about two times faster than the non-adaptive ones. The effect of the adaptive zone and smoothing zone on the impact pressure was also examined for the proposed algorithm.

Sensitivity of the ECMWF Model to Semi-Lagrangian Departure Point Iterations

Accurate estimation of the position of the departure points (d.p.) is crucial for the accuracy of a semi-Lagrangian NWP model. This calculation is often performed applying an implicit discretization to a kinematic equation solved by a fixed-point iteration scheme. A small number of iterations is typically used, assuming that this is sufficient for convergence. This assumption, derived from a past theoretical analysis, is revisited here. Analyzing the convergence of a generic d.p. iteration scheme and testing the ECMWF Integrated Forecast System (IFS) model, it is demonstrated that 2–3 iterations may not be sufficient for convergence to satisfactory accuracy in a modern high-resolution global model. Large forecast improvements can be seen by increasing the number of iterations. The extratropical geopotential error decreases and the simulated vertical structure of tropical cyclones improves. These findings prompted the implementation of an algorithm in which stopping criteria based on estimated convergence rates are used to “dynamically” stop d.p. iterations when an error tolerance criterion is satisfied. This is applied consistently to the nonlinear forecast, tangent linear, and adjoint models used by the ECMWF data assimilation system (4DVAR). Although the additional benefit of dynamic iteration is small, its testing reinforces the conclusion that a larger number of iterations is needed in regions of strong winds and shear. Furthermore, experiments suggest that dynamic iteration may prevent occasional 4DVAR failures in cases of strong stratospheric cross-polar flow in which the tangent linear model becomes unstable.

Two-Level Space–Time Domain Decomposition Methods for Flow Control Problems

For time-dependent control problems, the class of sub-optimal algorithms is popular and the parallelization is usually applied in the spatial dimension only. In the paper, we develop a class of fully-optimal methods based on space–time domain decomposition methods for some boundary and distributed control of fluid flow and heat transfer problems. In the fully-optimal approach, we focus on the use of an inexact Newton solver for the necessary optimality condition arising from the implicit discretization of the optimization problem and the use of one-level and two-level space–time overlapping Schwarz preconditioners for the Jacobian system. We show that the numerical solution from the fully-optimal approach is generally better than the solution from the sub-optimal approach in terms of meeting the objective of the optimization problem. To demonstrate the robustness and parallel scalability and efficiency of the proposed algorithm, we present some numerical results obtained on a parallel computer with a few thousand processors.

A block-coupled Finite Volume methodology for linear elasticity and unstructured meshes

The current article presents a fully coupled cell-centred Finite Volume solution methodology for linear elasticity and unstructured meshes. Details are given of the novel implicit discretisation and block coupled solution procedure, including use of a Finite Area method for face tangential gradient calculations and coupled treatment of non-orthogonal corrections. A number of 2-D and 3-D linear-elastic benchmark test cases are examined using hexahedral, tetrahedral and general polyhedral meshes; solution accuracy and efficiency are compared with that of a segregated procedure and a commercial Finite Element software, where the new method is shown to be faster in all cases.

A computationally efficient spectral method for modeling core dynamics

AbstractAn efficient, spectral numerical method is presented for solving problems in a spherical shell geometry that employs spherical harmonics in the angular dimensions and Chebyshev polynomials in the radial direction. We exploit the three‐term recurrence relation for Chebyshev polynomials that renders all matrices sparse in spectral space. This approach is significantly more efficient than the collocation approach and is generalizable to both the Galerkin and tau methodologies for enforcing boundary conditions. The sparsity of the matrices reduces the computational complexity of the linear solution of implicit‐explicit time stepping schemes to O(N) operations, compared to operations for a collocation method. The method is illustrated by considering several example problems of important dynamical processes in the Earth's liquid outer core. Results are presented from both fully nonlinear, time‐dependent numerical simulations and eigenvalue problems arising from the investigation of the onset of convection and the inertial wave spectrum. We compare the explicit and implicit temporal discretization of the Coriolis force; the latter becomes computationally feasible given the sparsity of the differential operators. We find that implicit treatment of the Coriolis force allows for significantly larger time step sizes compared to explicit algorithms; for hydrodynamic and dynamo problems at an Ekman number of , time step sizes can be increased by a factor of 3 to 16 times that of the explicit algorithm, depending on the order of the time stepping scheme. The implementation with explicit Coriolis force scales well to at least 2048 cores, while the implicit implementation scales to 512 cores.

A new family of fourth-order locally one-dimensional schemes for the three-dimensional wave equation

In this paper, we present a new family of locally one-dimensional (LOD) schemes with fourth-order accuracy in both space and time for the three-dimensional (3D) acoustic wave equation. It is well-known that high order explicit schemes offer efficient time steppings but with restrictive CFL conditions; implicit discretization gives unconditional stable schemes but with inefficient time steppings. Our LOD schemes can be seen as a compromise between explicit and implicit schemes, in the way that our schemes have more relax CFL conditions compared with traditional explicit schemes and only require solutions of tridiagonal systems of linear algebraic equations, which is more efficient than implicit schemes. Moreover, the new scheme is four-layer in time and three-layer in space, so that boundary conditions can be imposed in the classical way. The computations of the initial conditions for the three intermediate time layers are explicitly constructed. Furthermore, the stability condition is derived and the restriction for the time step is given explicitly. Finally, numerical examples are completed to show the performance of our new method.

Numerical approximations to a new phase field model for two phase flows of complex fluids

We derive a new phase field theory for immiscible mixtures of nematic liquid crystals and viscous fluids using the variational principle coupled with the generalized Onsager principle. A novel phase transition mechanism is implemented to couple the nematic liquid crystal phase with the viscous fluid phase to arrive at the dissipative hydrodynamic model for incompressible fluid mixtures. Through a delicate explicit–implicit numerical discretization, we develop a decoupled, linear scheme for a simplified version of the phase field model, as well as a coupled, nonlinear scheme for the full model. Both schemes are shown as unconditionally energy stable with consistent, discrete dissipative energy laws. Several numerical examples are presented to show the effectiveness of the new model and the new numerical schemes developed for it.