双时间步法在二维浅水方程求解中的应用

针对传统算法效率低的问题,将隐式双时间步法应用于求解二维浅水方程,建立了非结构网格下高效的有限体积模型。在应用双时间步法时,虚拟时间层中的定常问题采用高效的隐式LU-SGS（Lower-Upper Symmetric Gauss-Seidel）方法进行迭代求解。通过模拟计算4个典型算例以及与传统显式算法进行比较,对模型精度、效率及处理实际问题能力进行检验,分析了时间步长、内迭代次数对模型性能的影响。结果表明,双时间步法放宽了稳定性对时间步长的限制,时间步长可取到显式格式10倍以上,计算耗时减少了50%以上,模型具有良好精度与适应性,具有较好的推广应用价值。     

An interpretation of initial stress and strain methods,and numerical stability

This paper explores similarities and differences between initial stress and initial strain methods of analyses. The author takes the position that both methods, although perhaps conceptually different, are actually numerically the same. A simple approach for providing closed-form stability criteria for explicit and implicit time-marching procedures is described. For the case of associated, non-hardening viscoplasticity the approach provides criteria which are identical to those of Cormeau for explicit schemes. An implicit time-marching scheme which avoids the compliance matrix inversions at the constitutive level is also presented. Like many other implicit methods, this implicit time-marching scheme is unconditionally stable for θ≥½. The connection between initial strain plasticity and an explicit viscoplastic time-marching scheme is also addressed.     

Development of an implicit material point method for geotechnical applications

An implicit material point method (MPM), a variant of the finite element method (FEM), is presented in this paper. The key feature of MPM is that the spatial discretisation uses a set of material points, which are allowed to move freely through the background mesh. All history-dependent variables are tracked on the material points and these material points are used as integration points similar to the Gaussian points. A mapping and re-mapping algorithm is employed, to allow the state variables and other information to be mapped back and forth between the material points and background mesh nodes during an analysis. In contrast to an explicit time integration scheme utilised by most researchers, an implicit time integration scheme has been utilised here. The advantages of such an approach are twofold: firstly, it addresses the limitation of the time step size inherent in explicit integration schemes, thereby potentially saving significant computational costs for certain types of problems; secondly, it enables an improved algorithm accuracy, which is important for some constitutive behaviours, such as elasto-plasticity. The main purpose of this paper is to provide a unified MPM framework, in which both quasi-static and dynamic analyses can be solved, and to demonstrate the model behaviour. The implementation closely follows standard FEM approaches, where possible, to allow easy conversion of other FEM codes. Newton’s method is used to solve the equation of motion for both cases, while the formation of the mass matrix and the required updating of the kinematic variables are unique to the dynamic analysis. Comparisons with an Updated Lagrangian FEM and an explicit MPM code are made with respect to the algorithmic accuracy and time step size in a couple of representative examples, which helps to illustrate the relative performance and advantages of the implicit MPM. A geotechnical application is then considered, illustrating the capabilities of the proposed method when applied in the geotechnical field.     

Comparison of transmitting boundaries in dynamic finite element analyses using explicit time integration

Explicit time integration schemes provide an efficient solution to non-linear dynamic finite element analyses of geotechnical problems especially when high frequency response is important. Such explicit time integration schemes require one of two distinct transmitting boundary formulations to overcome the problem of radiation damping. These are the superposition boundary approach, which involves the cancellation of the reflected waves by combining the solutions of two different boundary conditions and the viscous boundary approach, which involves the absorption of incident wave energy by frequency independent viscous dashpots. The theoretical justification of these two approaches and their means of implementation are reviewed. The solutions obtained using the two different boundary approaches to the problem of a rigid massless circular footing vibrating on an elastic half-space are compared with an independent theoretical solution. The performance of the boundaries for problems involving step loading is also examined and the implications for any loading pattern with a non-zero time average are discussed.     

On the application of the maximum entropy meshfree method for elastoplastic geotechnical analysis

In this study, the Maximum Entropy Meshfree (MEM) method is employed for analysing geotechnical problems involving material nonlinearity, assuming small strains. The efficiency of the MEM method is evaluated through several solution schemes for the global governing equations as well as the local constitutive equations. The conventional implicit approach involving the Newton-Raphson method and an explicit adaptive dynamic relaxation technique are employed for solving the governing equations, while local constitutive equations are solved numerically as well as analytically. Two- and three-dimensional numerical experiments are performed to study the efficiency of different configurations of the solution scheme, which leads to some important conclusions about application of the MEM method in geotechnical problems.     

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.     

Computational procedures for non-linear three-dimensional analysis with some advanced constitutive laws

Computational procedures for implementing some constitutive models are described and introduced in three-, and two-dimensional finite element procedures; here variable moduli, Drucker–Prager, critical state and cap models are considered. Consistent numerical schemes are presented with applications to a number of example problems. These procedures can provide successful results with advanced constitutive laws for three-dimensional analysis of a wide range of non-linear problems.     

Numerical integration and FE implementation of a hypoplastic constitutive model

Hypoplastic constitutive equation based on nonlinear tensor functions possesses a failure surface but no yield surface. In this paper, we consider the numerical integration and FE implementation of a simple hypoplastic constitutive equation. The accuracy of several integration methods, including implicit and explicit methods, is examined by performing a set of triaxial compression tests. Adaptive explicit schemes show the best performance. In addition, the stress drift away from the failure surface is corrected with a predictor-corrector scheme, which is verified by two boundary value problems, i.e. rigid footing tests and slope stability.     

Numerical study on finite element implementation of hypoplastic models

A comprehensive numerical study on finite element implementation of hypoplastic models is presented. Two crucial aspects, local integration of the constitutive equations (the local problem) and forming tangent operators for Newton–Raphson iteration (the global problem), are investigated. For solving the local problem, different integration algorithms, including explicit and implicit methods, are examined using tri-axial compression tests and incremental stress response envelopes, as well as typical boundary value problems. For solving global problems, three different ways of generating the tangent operator are compared. The numerical evidences indicate that, in terms of accuracy, efficiency and robustness, explicit methods with substepping and error control are the best choices for constitutive integration of hypoplastic models while the so-called continuum tangent operators have certain advantages over two other types of numerically-generated consistent tangent operators.     

滑坡稳定分析传递系数法若干问题探讨

滑坡稳定性分析传递系数法有隐式和显式两种解法。当前对这两种解法的认识在稳定系数定义、满足的静力平衡方程类型、稳定系数计算公式表达方式、稳定系数大小关系和稳定系数计算误差等方面存在着不清楚或不一致之处。本文对此逐一进行了分析,并在此基础上提出了传递系数法解法选择的建议。     

A quasi-explicit modification of the two time level family of schemes for non-linear transient analysis

A quasi-explicit modification of well-known two time level implicit schemes for first-order, non-linear initial value problems is presented. Stable and accurate quadratic families for the quasi-explicit generalized trapezoidal method (GTM) and the generalized midpoint method (GMM) are derived. Accuracy, stability and oscillatory behaviour of the proposed algorithms are studied. Much higher stability and oscillatory levels are obtained for quasi-explicit modification of the generalized trapezoidal method than those of the generalized midpoint method.     

Evaluation of numerical stress-point algorithms on elastic–plastic models for unsaturated soils with hardening dependent on the degree of saturation

Constitutive models of unsaturated soils, and in particular those based on constitutive variables which include both degree of saturation and suction, are characterised by strong non linearities due to hydromechanical coupling. In this paper, a refined Runge–Kutta–Dormand–Prince explicit algorithm and a fully implicit Euler scheme are compared for the integration of the latter class of models. The explicit and implicit procedures have been tested along different hydromechanical paths, involving various hydraulic and mechanical external control conditions. Accuracy and efficiency of the algorithms have been investigated. The results confirm that substepping is mandatory for the explicit algorithm to converge regardless the initial step size and to remain sufficiently accurate. The value of the incremental hydromechanical work per unit volume was calculated during the explicit integration procedure. The numerical results show that the maximum size of the substep which can be adopted to meet a given tolerance depends on the gradient of the incremental work per unit volume. Therefore, the latter appears a good candidate to identify problematic integration steps in terms of convergence. Accuracy of the implicit algorithm also depends on the chosen step size, although the algorithm proved to be convergent in all the paths analysed.     

A coupled time integration algorithm for discontinuous deformation analysis using the numerical manifold method

To improve the computational efficiency of the numerical manifold method for discontinuous deformation simulations, a spatial-domain coupled explicit-implicit time integration algorithm is proposed. A subdomain partition algorithm based on a super manifold element is developed for the numerical manifold method to simulate dynamic motions of blocky rock mass. In different subdomains, explicit or implicit time integration method is employed respectively based on its contact and motion status. These subdomains interact through assembling the corresponding explicit or implicit time integration-based matrices of different rock blocks. The computational efficiency of the discontinuity system under dynamic loading is improved by partially diagonalizing the global matrices. Two verification examples of a sliding block along an inclined plane under a horizontal acceleration excitation and a multiblock system acted on by dynamic forces are studied to examine the accuracy of the proposed numerical method, respectively. A highly fractured rock mass situated on an inclined slope subjected to seismic excitations is then studied to show the computational efficiency of the developed algorithm. The simulated results are in good agreement with those from the versions using purely implicit or explicit time integration algorithm for the numerical manifold method. The computational efficiency is shown to be higher using the proposed algorithm, which demonstrates its potential for application in dynamic analysis of highly fractured rock masses.     

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.     

Evaluating the performance of an explicit dynamic relaxation technique in analyzing non-linear geotechnical engineering problems

Explicit dynamic relaxation is an efficient tool that has been used to solve problems involving highly non-linear differential equations. The key feature of this method is the ability to use explicit dynamic algorithms in solving static problems. Few attempts have been made to date to apply this technique in conventional geotechnical engineering. In this study, an algorithm that incorporates the application of a stiffness dependent time step scheme is proposed. The algorithm has been successfully used to solve 2D and 3D non-linear geotechnical engineering problems. To calibrate the developed algorithm, numerical simulations have been conducted for a strip and square footings supported by Mohr–Coulomb material. Performance of four different types of brick elements used in collapse load calculation is examined in terms of convergence speed and accuracy. In addition, the role of employing adaptive time steps in reducing the number of iterations needed for convergence is also evaluated.     

CS2: A PIECEWISE-LINEAR MODEL FOR LARGE STRAIN CONSOLIDATION

This paper presents a piecewise-linear finite-difference model for one-dimensional large strain consolidation called CS2. CS2 is developed using a fixed Eulerian co-ordinate system and constitutive relationships which are defined by discrete data points. The model is dimensionless such that solutions are independent of the initial height of the compressible layer and the absolute magnitude of the hydraulic conductivity of the soil. The capability of CS2 is illustrated using four example problems involving small strain, large strain, self-weight, and non-linear constitutive relationships. In each case, the performance of the model is comparable to other available analytical and numerical solutions. Using CS2, correction factors are developed for the conventional Terzaghi theory which account for the effect of vertical strain on computed values by elapsed time and maximum excess pore pressure during consolidation. © 1997 John Wiley & Sons, Ltd.     

A consistent formulation of a dilatant interface element

The paper considers a plane joint or interface element suitable for implementation into a standard non-linear finite element code. The element is intended to model discontinuities with rough contact surfaces, such as rock joints, where dilatant behaviour is present. Of particular concern is the formulation of a constitutive model which fully caters for all possible histories of opening, closing and sliding (accompained by dilation or contraction) in any direction. The non-linear incremental constitutive equations are formulated in a manner appropriate for a back-ward difference discretization in time along the path of loading. The advantage of such an approach is that no essential distinction need be drawn between opening, closing and sliding. Further, a convenient formulation of the constitutive equations is facilitated by representing the different contact conditions in relative displacement space. The state diagram in relative displacement space, however, changes from one time step to the next, and evolution equations for the updating must be formulated. These concepts are illustrated for two rock-joint models: a sawtooth asperity model and a limited dilation model. The models are based on a penalty formulation to enforce the contact constraints, and explicit equations for the tangent stiffness matrix and for the corrector step of the standard Newton–Raphson iterative algorithm are derived. These equations have been implemented as an user element into the finite element code ABAQUS⁷. Three examples are presented to illustrate the predictions of the formulation.     

Numerical simulation of undrained insertion problems in geotechnical engineering with the Particle Finite Element Method (PFEM)

The paper presents total-stress numerical analyses of large-displacement soil-structure interaction problems in geomechanics using the Particle Finite Element Method (PFEM). This method is characterized by frequent remeshing and the use of low order finite elements to evaluate the solution. Several important features of the method are: (i) a mixed formulation (displacement-mean pressure) stabilized numerically to alleviate the volumetric locking effects that are characteristic of low order elements when the medium is incompressible, (ii) a penalty method to prescribe the contact constraints between a rigid body and a deformable media combined with an implicit scheme to solve the tangential contact constraint, (iii) an explicit algorithm with adaptive substepping and correction of the yield surface drift to integrate the finite-strain multiplicative elasto-plastic constitutive relationship, and (iv) the mapping schemes to transfer information between successive discretizations. The performance of the method is demonstrated by several numerical examples, of increasing complexity, ranging from the insertion of a rigid strip footing to a rough cone penetration test. It is shown that the proposed method requires fewer computational resources than other numerical approaches addressing the same type of problems.     

A parallelizable method for two‐phase flows in naturally‐fractured reservoirs

A parallelizable, semi‐implicit numerical method is proposed for the study of naturally‐fractured reservoir systems. It has proved to be computationally efficient in producing accurate numerical solutions for the dual‐porosity model for immiscible, two‐phase flow in such reservoirs. The method combines hybridized mixed finite elements, a new version of the modified method of characteristics, a sophisticated operator‐splitting procedure for separating the pressure calculation in the fractures from that of the saturation, another operator splitting to handle the interaction of the matrix blocks and the fractures, and domain decomposition iterative procedures for both the pressure and the saturation. It permits moderately long time steps for the pressure and the saturation in the fractures and matrix blocks by using short, inexpensive microsteps to treat the transport portion of the saturation equation in the fractures. This paper is devoted to the formulation of the method and a discussion of numerical results for five‐spot and vertical cross‐section examples.     

Robust iterative schemes for non-linear poromechanics

We consider a non-linear extension of Biot’s model for poromechanics, wherein both the fluid flow and mechanical deformation are allowed to be non-linear. Specifically, we study the case when the volumetric stress and the fluid density are non-linear functions satisfying certain assumptions. We perform an implicit discretization in time (backward Euler) and propose two iterative schemes for solving the non-linear problems appearing within each time step: a splitting algorithm extending the undrained split and fixed stress methods to non-linear problems, and a monolithic L-scheme. The convergence of both schemes are shown rigorously. Illustrative numerical examples are presented to confirm the applicability of the schemes and validate the theoretical results.