Monotonicity conditions for control volume methods on uniform parallelogram grids in homogeneous media

Control volume methods are frequently used in porous media flow. This article gives an example on how one method, the Multipoint Flux Approximation method (MPFA), fails to satisfy the maximum principle for strong anisotropies or grid skewnesses, and develops conditions for when the monotonicity property holds for uniform parallelogram grids in homogeneous media. The conditions developed are applicable to any nine-point scheme in 2D or 27-point scheme in 3D, and is useful when the method produces a system matrix which is not an M-matrix.     

A segregated flow scheme to control numerical dispersion for multi-component flow simulations

One of the challenges for reservoir simulation is numerical dispersion. For waterflooding applications the effect is controlled due to the self-sharpening nature of a Buckley–Leverett shock. However, for multi-component flow simulations, incorrect wavespeeds can develop leading to the excessive smearing of fronts because of the coupling of compositional dispersion with the fractional flow. Rather than implementing a higher-order discretization method, we propose a simple scheme based on segregation-in-flow within a gridblock to control numerical dispersion. We extend the method originally proposed for polymer flooding to augmented waterflooding simulations in general as well as simulations of miscible or near miscible gas injection. For compositional simulations of gas injection, this is done through a coupled limited-flash/upstream-exclusion assumption. To test the scheme, an in-house streamline simulator has been modified and validated for modeling low-salinity floods as well as ternary two-phase displacements. Simulation results presented with and without segregation demonstrate the potential of the approach as a heuristic method to control numerical dispersion in multi-component flow simulations.     

Modeling piecewise planar fault discontinuities without element-partitioning in 3D reservoir-geomechanical models

Using the extended finite element method, faults can be introduced into a three-dimensional reservoir-geomechanical model without meshing to assess the potential for fault reactivation associated with industrial activities such as disposal of CO₂ or fluid extraction. The residuals to the governing equations include basis functions formed from the enrichment functions for strong and weak discontinuities. The traditional spatial integration scheme is based on an algorithm that partitions each enriched element into a collection of subtetrahedra, which is time and storage intensive. To avoid element partitioning, we adopt the homogeneous numerical integration scheme (referred as HNI hereafter) where the integration of homogeneous monomials over each polyhedron is converted into the integration of the same monomials over the one-dimensional edges of the polyhedron by using Stokes's theorem and Euler's homogeneous function theorem. The integrands for the strong and weak discontinuities are derived to implement the integration scheme in the three-dimensional reservoir-geomechanical model. An example of a jagged fault is presented to reveal the advantages of HNI scheme over the standard element partitioning approach. Several other examples that involve fluid flow, fault sliding, and fault sliding triggered by injection pressure are also presented to demonstrate that accurate and efficient computations are realized by the new integration scheme.     

非均匀层状介质一维波动方程精确解的有限差分算法

平面波的传播问题通常可以归结为一维波动方程的定解问题。在非均匀介质中,即使简单的一维波动方程也需要借助于数值方法获得近似解。3层5点古典差分格式是计算偏微分方程一种常用算法,作为一种显式迭代格式,需要满足稳定性条件 ,其中 为波速, 为空间采样间隔, 为时间采样间隔。当 时, ,古典差分格式达到临界稳定状态。在这种情况下,平面波在 时间内的传播距离恰好等于空间采样间隔,差分格式真实地反映了平面波的传播原理,因而可以得到一维波动方程的精确解。但是,由于在非均匀介质中存在不连续的波阻抗界面,此方法不适于计算非均匀介质的波场。为了将临界稳定情况下的古典差分格式推广应用至非均匀层状介质,提出了一种能够处理波阻抗界面的有限差分格式,并应用傅里叶分析法得到其稳定性条件。模型算例验证了此算法的正确性。     

Dynamic soil–structure interaction analysis using Lanczos vectors with adaptive time integration technique

An analytical procedure to obtain the response of soil–structure interaction problems, time domain is described. The procedure makes use of large domain for descritization along with co-ordinate transformation using Lanczos vectors. The responses are obtained in time domain using an adaptive direct integration method. The scheme has the ability to estimate errors due to temporal discretization as well as co-ordinate transformation. The procedure has been applied to half-space problems and non-convex domains for validation of the scheme, and the scheme obeys causality condition in both the situations. The present method has all the advantages of time domain scheme which is local both in space and time with small computational effort. Copyright © 1999 John Wiley & Sons, Ltd.     

孕震构造块体与相应地震区划分方法

可靠地划分地震区可奠定地震预测与地震危险性评价的地质基础,具有十分重要的意义。笔者等通过研究分析指出板内孕震构造块体侧向边界可由区域性大断层或由区域性大断层与板块边界界定,底边界为康拉德面或低速高导层;板间孕震构造块体为俯冲板块,可由区域性大断层和(或)板块边界约束;在同一个孕震构造块体和同一轮地震周期的地震具有内在联系。因此,地震区可定义为代表相应孕震构造块体地震活动的区域,其可表征该块体内源自锁固段破裂的地震活动。基于笔者等提出的孕震构造块体和相应地震区边界确定原则,把全球两大地震带(环太平洋地震带和欧亚地震带)划分为62个地震区;每个地震区的分区方案均通过了多锁固段脆性破裂理论的检验,这说明方案可靠。进而,笔者等归纳总结了地震区划分方法。     

Isogeometric analysis for unsaturated flow problems

Unsaturated flow problems in porous media often described by Richards’ equation are of great importance in many engineering applications. In this contribution, we propose a new numerical flow approach based on isogeometric analysis (IGA) for modeling the unsaturated flow problems. The non-uniform rational B-spline (NURBS) basis is utilized for spatial discretization whereas the stable implicit backward Euler method for time discretization. The nonlinear Richards’ equation is iteratively solved with the aid of the Newton–Raphson scheme. Owing to some desirable features of an efficient numerical flow approach, major advantages of the present formulation involve: (a) numerical oscillation at the wetting front can be avoided or facilitated, simply by using either an h-refinement or a lumped mass matrix technique; (b) higher-order exactness can be obtained due to the nature of the IGA features; (c) the approach is straightforward to implement and it does not need any transformation, e.g., Kirchhoff transformation or filter algorithm; and (d) in contrast to the Picard iteration scheme, which forms linear convergences, the proposed approach can however yield quadratic convergences by using the Newton–Raphson method for solving resultant nonlinear equations. Numerical model validation is analyzed by solving a three-dimensional unsaturated flow problem in soil, and its derived results are verified against analytical solutions. Numerical applications are then studied by considering three extensive examples with simple and complex configurations to further show the accuracy and applicability of the present IGA.     

Numerical modeling of the flow and transport of radionuclides in heterogeneous porous media

This paper is concerned with numerical methods for the modeling of flow and transport of contaminant in porous media. The numerical methods feature the mixed finite element method over triangles as a solver to the Darcy flow equation and a conservative finite volume scheme for the concentration equation. The convective term is approximated with a Godunov scheme over the dual finite volume mesh, whereas the diffusion–dispersion term is discretized by piecewise linear conforming triangular finite elements. It is shown that the scheme satisfies a discrete maximum principle. Numerical examples demonstrate the effectiveness of the methodology for a coupled system that includes an elliptic equation and a diffusion–convection–reaction equation arising when modeling flow and transport in heterogeneous porous media. The proposed scheme is robust, conservative, efficient, and stable, as confirmed by numerical simulations.      

Comparison between the Gauss’ law method and the zero current method to calculate multi-species ionic diffusion in saturated uncharged porous materials

A novel numerical method based on the finite element approach is established for the zero current method approach for calculating multi-species ionic diffusion. The proposed numerical method uses the direct calculation of the coupled set of equations in favor of the staggering approach. A one-step truly implicit time stepping scheme is adopted together with an implementation of a modified Newton–Raphson iteration scheme for search of equilibrium at each considered time step calculation. Results from the zero current case are compared with existing results from the solutions of the more general Gauss’ law method.     

Variational method for objective analysis of scalar variable and its derivative

In this study real time data have been used to compare the standard and triangle method by performing the objective analysis of mean sea level pressure. In the standard method, derivative fields are obtained from the grid point data using finite difference scheme whereas in the triangle method, a set of non-overlapping triangles are formed from the observations and the scalar and the spatial derivatives are computed directly at the centroid of each of the non-overlapping triangles. These scalars and their derivatives are then mapped to uniform grids by using the standard method. It has been found that objectively analysed scalar field obtained using standard method is superior to the scalar field derived by the triangle method, whereas the derivative fields produced by triangle method are superior to the derivative fields produced using standard method. A variational objective analysis scheme has been developed and an experiment has been carried out with depression case of June (11–15) 2004. It is found that the new scheme (variational) is able to extract the better parts of both triangle and standard methods. The results of this study will be useful in carrying out diagnostic calculations that involve derivative estimates.     

边坡治理群决策的二维足码定位法研究

在边坡治理决策过程中,不仅涉及边坡本身的条件和状态等因素,而且涉及经济、技术、环境及专家和决策者的个人素质等因素,因此边坡治理方案选择是一个复杂的系统工程。在实际工程中,往往同时有几个治理方案供选择,必须对它们进行优化评价,从中选出较满意的方案。针对边坡治理方案评价是一个群决策问题,提出了边坡治理群决策的二维足码定位法,利用简单的二维足码坐标图或评价函数,可以确定边坡治理的最优方案。以溪沟口滑坡治理为例,分析了该方法的实际应用情况。     

Using a shifted Grünwald-Letnikov scheme for the Caputo derivative to study anomalous solute transport in porous medium

We propose an extension of the shifted Grünwald-Letnikov method to solve fractional partial differential equations in the Caputo sense with arbitrary fractional order derivative α and with an advective term. The method uses the relation between Caputo and Riemann-Liouville definitions, the shifted Grünwald-Letnikov, and the traditional backward and forward finite difference method. The stability of the method is investigated for the implicit and explicit scheme with homogeneous boundary conditions, and a stability criterion is found for the advective-dispersive equation. An application of the method is used to solve contaminant diffusion and advective-dispersive problems. The numerical solution for the fractional diffusion and fractional advection-dispersion is compared with their respective analytical solutions for different time and space grid refinements. The diffusion simulation exhibited a good fit between the analytical and numerical solutions, with the explicit scheme going from stable to unstable as the time and space refinement changes. The fractional advection-dispersion application produced small deviations from the analytical solution. These deviations, however, are analogous to the numerical dispersions encountered in conventional finite difference solutions of the advection-dispersion equation. The new method is also compared with the traditional L2 method. Notably, an example that involves asymmetrical fractional conditions, a fractional diffusivity that depends on time, and a source term show how the methods compare. Overall, this study assesses the quality and easiness of use of the numerical method.     

基于三维有限元的地下连续墙深基坑逆作法施工方案设计

基于三维有限元,对采用"逆作法"施工的地下连续墙深基坑进行了设计方案研究,按变形控制提出了设计方法,并以天津某基坑为例进行了计算.首先,对结构进行了竖向荷载作用下的结构设计计算,然后,在结构计算的基础上提出了3种基坑施工方案,分别建立了三维有限元模型.对各方案的不同工况进行计算模拟,获得了各方案在不同工况下的连续墙和楼板的应力应变结果,找出了施工的控制工况和受力的薄弱部位.最后,结合实际情况分析了各方案的优点和不足之处,对方案进行了优选,可为最终设计方案的确定提供参考.     

Short Communication: A numerical method for diffusive transport with moving boundaries and discontinuous material properties

A method for the numerical simulation of diffusive transport with moving boundaries is developed and tested. The variable domain is mapped onto a fixed region, which introduces a term of convective form to the transformed governing equation. The resulting convection/diffusion equation is solved by a finite-difference method. An ‘Immersed Interface’ Method (IIM) is introduced in order to retain second-order accuracy near discontinuities in material properties, where the solution is not smooth. The method performs well in benchmark calculations against an analytical solution. The IIM scheme is capable of treating a strong discontinuity in the gradient, and it is readily extended to two or three dimensions. The methods are illustrated through a calculation for the temperature profile in a growing continental ice sheet, in which the thermal properties are discontinuous at the rock/ice interface. © 1997 by John Wiley & Sons, Ltd.     

Analysis of Slope Stability and Software Development Based on Single-Grid and Two-Grid Finite Element Methods

Slope stability has been a long-standing problem in geotechnical practice. The difficulty in analyzing slope stability may be attributed to geological complexity and versatile environments, but the limitations of analysis methods should not be ignored. In this work, the conventional limit equilibrium methods and the FE stress-based limit equilibrium methods with special emphasis on their limitations are reviewed, and the advantages of the finite element method with the capability of strength reduction are pointed out and illustrated by using some slope examples. Particularly, a flowchart of decision tree on selecting single-grid or two-grid search scheme is proposed for the shear strength reduction finite element method of slope stability analysis. Based on some slope examples, it is demonstrated that the automatic selection or switch between the single-grid scheme and the two-grid scheme can be remarkably superior to the single mesh scheme from the perspectives of computational cost and software development.     

High‐order ADE scheme for solving the fluid diffusion equation in non‐uniform grids and its application in coupled hydro‐mechanical simulation

To improve the stability and efficiency of explicit technique, one proposed method is to use an unconditionally stable alternating direction explicit (ADE) scheme. However, the standard ADE scheme is only moderately accurate and restricted to uniform grids. This paper derives a novel high‐order ADE scheme capable of solving the fluid diffusion equation in non‐uniform grids. The new scheme is derived by performing a fourth‐order finite difference approximation to the spatial derivatives of the diffusion equation in non‐uniform grid. The implicit Crank‐Nicolson technique is then applied to the resulting approximation, and the subsequent equation is split into two alternating direction sweeps, giving rise to a new high‐order ADE scheme. Because the new scheme can be potentially applied in coupled hydro‐mechanical (H‐M) simulation, the pore pressure solutions from the new scheme are then sequentially coupled with an existing geomechanical simulator in the computer program Fast Lagrangian Analysis of Continua. This coupling procedure is called the sequentially explicit coupling technique based on the fourth‐order ADE scheme (SEA‐4). Verifications of well‐known consolidation problems showed that the new ADE scheme and SEA‐4 can reduce computer runtime by 46% to 75% to that of Fast Lagrangian Analysis of Continua's basic scheme. At the same time, the techniques still maintained average percentage error of 1.6% to 3.5% for pore pressure and 0.2% to 1.5% for displacement solutions and were still accurate under typical grid non‐uniformities. This result suggests that the new high‐order ADE scheme can provide an efficient explicit technique for solving the flow equation of a coupled H‐M problem, which will be beneficial for large‐scale and long‐term H‐M problems in geoengineering.     

Modeling of fluid–solid interaction in granular media with coupled lattice Boltzmann/discrete element methods: application to piping erosion

In this article, we present a numerical method to deal with fluid–solid interactions and simulate particle–fluid systems as encountered in soils. This method is based on a coupling between two methods, now widely used in mechanics of granular media and fluid dynamics respectively: the discrete element (DE) method and the lattice Boltzmann (LB) method. The DE method is employed to model interactions between particles, whereas the LB method is used to describe an interstitial Newtonian fluid flow. The coupling presented here is a full one in the sense that particle motions act on fluid flow and reciprocally. This article presents in details each of the two methods and the principle of the coupling scheme. Determination of hydrodynamic forces and torques is also detailed, and the treatment of boundaries is explained. The coupled method is finally illustrated on a simple example of piping erosion, which puts in evidence that the combined LB–DE scheme constitutes a promising tool to study coupled problems in geomechanics. Copyright © 2011 John Wiley & Sons, Ltd.     

Impact of difference accuracy on computational properties of vertical grids for a nonhydrostatic model

Starting from nonhydrostatic Boussinesq approximation equations, a general method is introduced to deduce the dispersion relationships. A comparative investigation is performed on inertia-gravity wave with horizontal lengths of 100, 10 and 1 km. These are examined using the second-order central difference scheme and the fourth-order compact difference scheme on vertical grids that are currently available from the perspectives of frequency, horizontal and vertical component of group velocity. These findings are compared to analytical solutions. The obtained results suggest that whether for the second-order central difference scheme or for the fourth-order compact difference scheme, Charny–Phillips and Lorenz (L) grids are suitable for studying waves at the above-mentioned horizontal scales; the Lorenz time-staggered and Charny–Phillips time staggered (CPTS) grids are applicable only to the horizontal scales of less than 10 km, and N grid (unstaggered grid) is unsuitable for simulating waves at any horizontal scale. Furthermore, by using fourth-order compact difference scheme with higher difference precision, the errors of frequency and group velocity in horizontal and vertical directions produced on all vertical grids in describing the waves with horizontal lengths of 1, 10 and 100 km cannot inevitably be decreased. So in developing a numerical model, the higher-order finite difference scheme, like fourth-order compact difference scheme, should be avoided as much as possible, typically on L and CPTS grids, since it will not only take many efforts to design program but also make the calculated group velocity in horizontal and vertical directions even worse in accuracy.     

基于“数字水调”规划的黄河水量调度方案编制

在“数字水调”规划下,将自适应控制思想引入黄河水调方案编制工作中,研究开发了黄河水量调度方案的自适应编制法。该法能够跟踪适时数据,计算分析系统状态,及时调整调度期内水库泄流计划以及用户引水计划,以达到“防凌、防断流、公平分水”的调度目标。运用该法编制了2001至2002调度年的水调方案并与当年实际调度方案进行了对比,结果显示自适应方案编制法能够快速高效的编制年、月、旬方案,编制的方案完全满足调度要求。     

On the use of enriched finite element method to model subsurface features in porous media flow problems

In this paper, a new enrichment scheme is proposed to model fractures and other conduits in porous media flow problems. Inserting this scheme into a partition of unity based method results in a new numerical method that does not require the mesh to honor the specific geometry of these subsurface features. The new scheme involves a specially designed integration procedure and enrichment functions, which can capture effects of local heterogeneity introduced by subsurface features on the pressure solution. The new method is also capable of modeling fractures with low as well as high conductivity. Another feature of the proposed scheme is that, even though two enrichment functions are used to model the permeability change at the two rock/fracture interfaces of a fracture, only one element partition is made for numerical integration. To demonstrate the accuracy and effectiveness of the proposed approach, production problems for wells that were stimulated or completed by longitudinal fracture, transverse fractures, and perforations are studied.