Mixed finite element methods and higher order temporal approximations for variably saturated groundwater flow

Richards’ equation (RE) is commonly used to model flow in variably saturated porous media. However, its solution continues to be difficult for many conditions of practical interest. Among the various time discretizations applied to RE, the method of lines (MOL) has been used successfully to introduce robust, accurate, and efficient temporal approximations. At the same time, a mixed-hybrid finite element method combined with an adaptive, higher order time discretization has shown benefits over traditional, lower order temporal approximations for modeling single-phase groundwater flow in heterogeneous porous media. Here, we extend earlier work for single-phase flow and consider two mixed finite element methods that have been used previously to solve RE using lower order time discretizations with either fixed time steps or empirically based adaption. We formulate the two spatial discretizations within a MOL context for the pressure head form of RE as well as a fully mass-conservative version. We conduct several numerical experiments for both spatial discretizations with each formulation, and we compare the higher order, adaptive time discretization to a first-order approximation with formal error control and adaptive time step selection. Based on the numerical results, we evaluate the performance of the methods for robustness and efficiency.     

An evaluation of temporally adaptive transformation approaches for solving Richards' equation

Developing robust and efficient numerical solution methods for Richards' equation (RE) continues to be a challenge for certain problems. We consider such a problem here: infiltration into unsaturated porous media initially at static conditions for uniform and non-uniform pore size media. For ponded boundary conditions, a sharp infiltration front results, which propagates through the media. We evaluate the resultant solution method for robustness and efficiency using combinations of variable transformation and adaptive time-stepping methods. Transformation methods introduce a change of variable that results in a smoother solution, which is more amenable to efficient numerical solution. We use adaptive time-stepping methods to adjust the time-step size, and in some cases the order of the solution method, to meet a constraint on nonlinear solution convergence properties or a solution error criterion. Results for three test problems showed that adaptive time-stepping methods provided robust solutions; in most cases transforming the dependent variable led to more efficient solutions than untransformed approaches, especially as the pore-size uniformity increased; and the higher-order adaptive time integration method was robust and the most efficient method evaluated.     

A spatially and temporally adaptive solution of Richards’ equation

Efficient, robust simulation of groundwater flow in the unsaturated zone remains computationally expensive, especially for problems characterized by sharp fronts in both space and time. Standard approaches that employ uniform spatial and temporal discretizations for the numerical solution of these problems lead to inefficient and expensive simulations. In this work, we solve Richards’ equation using adaptive methods in both space and time. Spatial adaption is based upon a coarse grid solve and a gradient error indicator using a fixed-order approximation. Temporal adaption is accomplished using variable order, variable step size approximations based upon the backward difference formulas up to fifth order. Since the advantages of similar adaptive methods in time are now established, we evaluate our method by comparison with a uniform spatial discretization that is adaptive in time for four different one-dimensional test problems. The numerical results demonstrate that the proposed method provides a robust and efficient alternative to standard approaches for simulating variably saturated flow in one spatial dimension.     

Calculating actual crop evapotranspiration under soil water stress conditions with appropriate numerical methods and time step

Calculation of actual crop evapotranspiration under soil water stress conditions is crucial for hydrological modeling and irrigation water management. Results of actual evapotranspiration depend on the estimation of water stress coefficient from soil water storage in the root zone, which varies with numerical methods and time step used. During soil water depletion periods without irrigation or precipitation, the actual crop evapotranspiration can be calculated by an analytical method and various numerical methods. We compared the results from several commonly used numerical methods, including the explicit, implicit and modified Euler methods, the midpoint method, and the Heun's third‐order method, with results of the analytical method as the bench mark. Results indicate that relative errors of actual crop evapotranspiration calculated with numerical methods in one time step are independent of the initial soil water storage in the range of soil water stress. Absolute values of relative error decrease with the order of numerical methods. They also decrease with the number of time step, which can ensure the numerical stability of successive simulation of soil water balance. Considering the calculation complexity and calculation errors caused by numerical approximation for different time step and maximum crop evapotranspiration, the explicit Euler method is recommended for the time step of 1 day (d) or 2 d for maximum crop evapotranspiration less than 5 mm/d, the midpoint method or the modified Euler method for the time step of up to one week or 10 d for maximum crop evapotranspiration less than 5 mm/d, and the Heun's third‐order method for the time step of up to 15 d. Copyright © 2011 John Wiley & Sons, Ltd.     

A posteriori local error estimation and adaptive time-stepping for newmark integration in dynamic analysis

A simple a posteriori local error estimate for Newmark time integration schemes in dynamic analysis is presented, based on the concept of a so called ‘post-processing’ technique. In conjunction with the error estimate, an adaptive time-stepping algorithm is described, which adjusts the time step size so that the local error of each time step is within a prescribed error tolerance. Numerical examples given in the paper indicate that the error estimate is asymptotically convergent, computationally efficient and convenient, and the adaptive time-stepping scheme can predict a nearly optimal step size from time to time, thus making the numerical solution reliable in an efficient manner.     

Implementation and application of the unconditionally stable explicit parametrically dissipative KR‐α method for real‐time hybrid simulation

In real‐time hybrid simulations (RTHS) that utilize explicit integration algorithms, the inherent damping in the analytical substructure is generally defined using mass and initial stiffness proportional damping. This type of damping model is known to produce inaccurate results when the structure undergoes significant inelastic deformations. To alleviate the problem, a form of a nonproportional damping model often used in numerical simulations involving implicit integration algorithms can be considered. This type of damping model, however, when used with explicit integration algorithms can require a small time step to achieve the desired accuracy in an RTHS involving a structure with a large number of degrees of freedom. Restrictions on the minimum time step exist in an RTHS that are associated with the computational demand. Integrating the equations of motion for an RTHS with too large of a time step can result in spurious high‐frequency oscillations in the member forces for elements of the structural model that undergo inelastic deformations. The problem is circumvented by introducing the parametrically controllable numerical energy dissipation available in the recently developed unconditionally stable explicit KR‐α method. This paper reviews the formulation of the KR‐α method and presents an efficient implementation for RTHS. Using the method, RTHS of a three‐story 0.6‐scale prototype steel building with nonlinear elastomeric dampers are conducted with a ground motion scaled to the design basis and maximum considered earthquake hazard levels. The results show that controllable numerical energy dissipation can significantly eliminate spurious participation of higher modes and produce exceptional RTHS results. Copyright © 2014 John Wiley & Sons, Ltd.     

Efficient simulation of geothermal processes in heterogeneous porous media based on the exponential Rosenbrock–Euler and Rosenbrock-type methods

Simulation of geothermal systems is challenging due to coupled physical processes in highly heterogeneous media. Combining the exponential Rosenbrock–Euler method and Rosenbrock-type methods with control-volume (two-point flux approximation) space discretizations leads to efficient numerical techniques for simulating geothermal systems. In terms of efficiency and accuracy, the exponential Rosenbrock–Euler time integrator has advantages over standard time-discretization schemes, which suffer from time-step restrictions or excessive numerical diffusion when advection processes are dominating. Based on linearization of the equation at each time step, we make use of matrix exponentials of the Jacobian from the spatial discretization, which provide the exact solution in time for the linearized equations. This is at the expense of computing the matrix exponentials of the stiff Jacobian matrix, together with propagating a linearized system. However, using a Krylov subspace or Léja points techniques make these computations efficient.The Rosenbrock-type methods use the appropriate rational functions of the Jacobian of the ODEs resulting from the spatial discretization. The parameters in these schemes are found in consistency with the required order of convergence in time. As a result, these schemes are A-stable and only a few linear systems are solved at each time step. The efficiency of the methods compared to standard time-discretization techniques are demonstrated in numerical examples.     

Real‐time hybrid testing using the unconditionally stable explicit CR integration algorithm

Real‐time hybrid testing combines experimental testing and numerical simulation, and provides a viable alternative for the dynamic testing of structural systems. An integration algorithm is used in real‐time hybrid testing to compute the structural response based on feedback restoring forces from experimental and analytical substructures. Explicit integration algorithms are usually preferred over implicit algorithms as they do not require iteration and are therefore computationally efficient. The time step size for explicit integration algorithms, which are typically conditionally stable, can be extremely small in order to avoid numerical stability when the number of degree‐of‐freedom of the structure becomes large. This paper presents the implementation and application of a newly developed unconditionally stable explicit integration algorithm for real‐time hybrid testing. The development of the integration algorithm is briefly reviewed. An extrapolation procedure is introduced in the implementation of the algorithm for real‐time testing to ensure the continuous movement of the servo‐hydraulic actuator. The stability of the implemented integration algorithm is investigated using control theory. Real‐time hybrid test results of single‐degree‐of‐freedom and multi‐degree‐of‐freedom structures with a passive elastomeric damper subjected to earthquake ground motion are presented. The explicit integration algorithm is shown to enable the exceptional real‐time hybrid test results to be achieved. Copyright © 2008 John Wiley & Sons, Ltd.     

Second-order characteristic methods for advection–diffusion equations and comparison to other schemes

We develop two characteristic methods for the solution of the linear advection diffusion equations which use a second order Runge–Kutta approximation of the characteristics within the framework of the Eulerian–Lagrangian localized adjoint method. These methods naturally incorporate all three types of boundary conditions in their formulations, are fully mass conservative, and generate regularly structured systems which are symmetric and positive definite for most combinations of the boundary conditions. Extensive numerical experiments are presented which compare the performance of these two Runge–Kutta methods to many other well perceived and widely used methods which include many Galerkin methods and high resolution methods from fluid dynamics.     

A high order Godunov scheme with constrained transport and adaptive mesh refinement for astrophysical and geophysical MHD

We propose to extend the well-known MUSCL-Hancock scheme for Euler equations to the induction equation modeling the magnetic field evolution in kinematic dynamo problems. The scheme is based on an integral form of the underlying conservation law which, in our formulation, results in a “finite-surface” scheme for the induction equation. This naturally leads to the well-known “constrained transport” method, with additional continuity requirement on the magnetic field representation. The second ingredient in the MUSCL scheme is the predictor step that ensures second order accuracy both in space and time. We explore specific constraints that the mathematical properties of the induction equations place on this predictor step, showing that three possible variants can be considered. We show that the most aggressive formulations reach the same level of accuracy than the other ones, at a lower computational cost. More interestingly, these schemes are compatible with the Adaptive Mesh Refinement (AMR) framework. It has been implemented in the AMR code RAMSES. It offers a novel and efficient implementation of a second order scheme for the induction equation. The scheme is then adaptated to solve for the full MHD equations using the same methodology. Through a series of test problems, we illustrate the performances of this new code using two different MHD Riemann solvers (Lax–Friedrich and Roe) and the need of the Adaptive Mesh Refinement capabilities in some cases. Finally, we show its versatility by applying it to the ABC dynamo problem and to the collapse of a magnetized cloud core.     

Modified precise time step integration method of structural dynamic analysis

The precise time step integration method proposed for linear time-invariant homogeneous dynamic systems can provide precise numerical results that approach an exact solution at the integration points. However, difficulty arises when the algorithm is used for non-homogeneous dynamic systems, due to the inverse matrix calculation and the simulation accuracy of the applied loading. By combining the Gaussian quadrature method and state space theory with the calculation technique of matrix exponential function in the precise time step integration method, a new modified precise time step integration method （e.g., an algorithm with an arbitrary order of accuracy） is proposed. In the new method, no inverse matrix calculation or simulation of the applied loading is needed, and the computing efficiency is improved. In particular, the proposed method is independent of the quality of the matrix H. If the matrix H is singular or nearly singular, the advantage of the method is remarkable. The numerical stability of the proposed algorithm is discussed and a numerical example is given to demonstrate the validity and efficiency of the algorithm.     

ASSESSMENT OF MACCORMACK SCHEME IN UNSTEADY SEDIMENT LADEN FLOW COMPUTATION

１ＩＮＴＲＯＤＵＣＴＩＯＮＲｉｖｅｒｓｉｎＴａｉｗａｎａｒｅｒｅｌａｔｉｖｅｌｙｓｔｅｅｐｅｒｃｏｍｐａｒｅｄｔｏｔｈｏｓｅｉｎｏｔｈｅｒｃｏｎｔｉｎｅｎｔ．Ｌｏｃａｌｙｏｃｕｒｅｄｓｕｐｅｒｃｒｉｔｉｃａｌｆｌｏｗａｒｅｆａｉｒｌｙｃｏｍｍｏｎｉｎ...     

求解运动方程的一种不等时间步长的显式数值积分方法

在波动有限元模拟中, 若采用传统的显式数值积分方法求解运动方程, 计算时间步长需采用计算区内满足稳定条件要求的最小时间步长. 然而, 对于大部分计算区域, 这一时间步长过小, 是不必要的. 本文提出了一种不等时间步长的显式数值积分方法, 其基本思想是不同的计算区域采用满足各自稳定条件的计算时间步长. 最后, 本文通过数值试验检验了这一方法的可行性及其对数值计算精度的影响.      

THREE-DIMENSIONAL NUMERICAL MODELING OF COHESIVE SEDIMENT TRANSPORT BY TIDAL CURRENT IN PEARL RIVER ESTUARV

IINTRODUCTIONEstUariesareprominentcoastalfeatUres.Estuariesareofgreateconomicssignificancetomankind.Attheseareas,manyharborsandwaterchannelshavetobebuiltforeconomicpurposes.ThedesignandconstrUctionofcoastalstrUctUresinestUariesrequireknowledgeofhydrodynamicsaswellassedimenttransportinsuchregions.ThenatUreofestuariesiscontrolledbyvariouscoastalhydrodynamicprocesses.Undertheactionofhydrodynamics,sedimentdepositionsorerosionswilloccurinestuariesornearcoastalstrUCtures.Tomaintainnavigati…     

Numerical mass balance for soil-moisture transport problems

The Galerkin finite element method coupled with the Crank-Nicolson time advance procedure is often used as a numerical analog for unsaturated soil-moisture transport problems. The Crank-Nicolson procedure leads to numerical mass balance problems which results in instability. A new temporal and spatial integration procedure is proposed that exactly satisfies mass balance for the approximating function used. This is accomplished by fitting polynomials continuously throughout the time and space domain and integrating the governing differential equations. To reduce computational effort, the resulting higher order polynomials are reduced to quadratic and linear piece-wise continuous polynomial approximation functions analogous to the finite element approach. Results indicate a substantial improvement in accuracy over the combined Galerkin and Crank-Nicolson methods when comparing to simplified problems where analytical solutions are available.     

Computational methods for reactive transport modeling: A Gibbs energy minimization approach for multiphase equilibrium calculations

We present a numerical method for multiphase chemical equilibrium calculations based on a Gibbs energy minimization approach. The method can accurately and efficiently determine the stable phase assemblage at equilibrium independently of the type of phases and species that constitute the chemical system. We have successfully applied our chemical equilibrium algorithm in reactive transport simulations to demonstrate its effective use in computationally intensive applications. We used FEniCS to solve the governing partial differential equations of mass transport in porous media using finite element methods in unstructured meshes. Our equilibrium calculations were benchmarked with GEMS3K, the numerical kernel of the geochemical package GEMS. This allowed us to compare our results with a well-established Gibbs energy minimization algorithm, as well as their performance on every mesh node, at every time step of the transport simulation. The benchmark shows that our novel chemical equilibrium algorithm is accurate, robust, and efficient for reactive transport applications, and it is an improvement over the Gibbs energy minimization algorithm used in GEMS3K. The proposed chemical equilibrium method has been implemented in Reaktoro, a unified framework for modeling chemically reactive systems, which is now used as an alternative numerical kernel of GEMS.     

Evaluation of numerical time‐integration schemes for real‐time hybrid testing

We present a comparison of methods for the analysis of the numerical substructure in a real‐time hybrid test. A multi‐tasking strategy is described, which satisfies the various control and numerical requirements. Within this strategy a variety of explicit and implicit time‐integration algorithms have been evaluated. Fully implicit schemes can be used in fast hybrid testing via a digital sub‐step feedback technique, but it is shown that this approach requires a large amount of computation at each sub‐step, making real‐time execution difficult for all but the simplest models. In cases where the numerical substructure poses no harsh stability condition, it is shown that the Newmark explicit method offers advantages of speed and accuracy. Where the stability limit of an explicit method cannot be met, one of the several alternatives may be used, such as Chang's modified Newmark scheme or the α‐operator splitting method. Appropriate methods of actuator delay compensation are also discussed. Copyright © 2008 John Wiley & Sons, Ltd.     

Perfectly matched layer boundary conditions for the second-order acoustic wave equation solved by the rapid expansion method

We derive a governing second-order acoustic wave equation in the time domain with a perfectly matched layer absorbing boundary condition for general inhomogeneous media. Besides, a new scheme to solve the perfectly matched layer equation for absorbing reflections from the model boundaries based on the rapid expansion method is proposed. The suggested scheme can be easily applied to a wide class of wave equations and numerical methods for seismic modelling. The absorbing boundary condition method is formulated based on the split perfectly matched layer method and we employ the rapid expansion method to solve the derived new perfectly matched layer equation. The use of the rapid expansion method allows us to extrapolate wavefields with a time step larger than the ones commonly used by traditional finite-difference schemes in a stable way and free of dispersion noise. Furthermore, in order to demonstrate the efficiency and applicability of the proposed perfectly matched layer scheme, numerical modelling examples are also presented. The numerical results obtained with the put forward perfectly matched layer scheme are compared with results from traditional attenuation absorbing boundary conditions and enlarged models as well. The analysis of the numerical results indicates that the proposed perfectly matched layer scheme is significantly effective and more efficient in absorbing spurious reflections from the model boundaries.     

A forward particle tracking Eulerian–Lagrangian Localized Adjoint Method for solution of the contaminant transport equation in three dimensions

The contaminant transport equation is solved in three dimensions using the Eulerian–Lagrangian Localized Adjoint Method (ELLAM). Trilinear and finite volume test functions defined by the characteristics of the governing equation are employed and compared. Integrations are simplified by forward tracking of integration points along the characteristics. The resulting equations are solved using a preconditioned conjugate gradient method. The algorithm is coupled to a block-centered finite difference approximation of the groundwater flow equation similar to that used in the popular MODFLOW code. The ELLAM is tested by comparison with 1D and 3D analytic solutions. The method is then applied with random, spatially correlated hydraulic conductivities in a simulation of a tracer experiment performed on Cape Cod, Massachusetts. The linear test function ELLAM was found to perform better than the finite volume ELLAM. Both ELLAM formulations were found to be robust, computationally efficient and relatively straightforward to implement. When compared to traditional particle tracking and characteristics codes commonly used with MODFLOW, the ELLAM retains the computational advantages of traditional characteristic methods with the added advantage of good mass conservation.     

关于Newmark-β法机理的一种解释

一般认为,Newmark-β法属于积分类型的动力数值分析方法,和基于荷载分段插值类型的数值方法不是相同类型的方法.在本文中,研究了这两类方法之间的关系,以最常使用的两种Newmark方法-平均常加速度法和线性加速度法为例,从Newmark基本假定出发推导出这两种方法所具有的荷载分布模式.结果发现:平均常加速度法和线性加...