二阶声波方程频域PML边界条件及频域变网格步长并行计算

研究了二阶声波方程频域PML边界条件和频域变网格并行计算技术。PML边界是一种较为理想的吸收边界方法,多用在求解时域应力速度方程中,但对于频域声波正演,二阶位移方程更常用。从一阶声波方程PML吸收边界条件导出频域二阶位移方程PML边界条件,模拟算例得到的频率切片、时间切片和地震记录对比都说明该边界条件吸收效果很好。频域单炮正演不同频率间是独立的,据此低频部分采用大网格计算,高频采用小网格,实现变网格步长计算技术,这是较时间域正演的一个优势,在保证模拟质量的同时,减少计算量和内存消耗。     

基于联系累积前景理论的基坑支护方案决策

基坑支护方案决策受决策者的风险态度和指标的随机模糊性影响与控制,是一个复杂的不确定性决策难题。应用联系数与累积前景理论的耦合方法,探讨了基坑支护方案的联系累积前景理论优选新模型,即首先基于原始规范化的决策信息确定正、负理想方案,再基于联系数理论定量统一描述待选方案与正、负理想方案间的确定和不确定关系,以构建基于联系前景效用价值函数的综合决策模型,并求解最优的权向量,确定合理的基坑支护方案。实例应用和与其他方法结果对比表明,文中方法应用于基坑支护方案优选是有效可行的,且能定量表达待选方案与理想方案间的联系和转化态势,避免了基坑支护方案决策过程的非理性。     

二维浅水明流的一种二阶高性能算法

为了精确求解二维浅水方程组,在由任意三角形构成的无结构网格上建立了有限体积MUSCL算法,这是文献[1]中有限体积一阶Osher格式的一类二阶推广.通过采用预测校正二步时间积分和单元内引入坡度限制,该算法在时空上均可达到二阶精度.跨单元边界的法向数值通量采用通量向量分裂公式(FVS)计算.最后,通过长江口南支潮流计算和一维瞬时溃坝模型算例来说明该格式的优良性能及在浅水流动计算中的应用.     

Dispersion and cost analysis of some finite difference schemes in one-parameter acoustic wave modeling

Given a precision threshold to be imposed on the group velocity error and a finite difference scheme for the acoustic wave equation, it is possible to determine time-step and grid-spacing in an optimal manner, i.e., so as to minimize the computational cost. Using this optimal cost as a criterion, it becomes easy to compare schemes for efficiency in homogeneous media. Heterogeneous media with constant density can be accommodated to a certain extent by minimizing the cost over a range of Courant numbers. Such analysis shows that, amongst the second-order Taylor series schemes in time, higher-order schemes are generally more efficient than lower-order schemes. However, this result does not extend to very high order schemes.     

辅助微分方程完全匹配层在声波方程数值模拟中的应用

在地震波数值模拟中,为提高算法精度,需要使用高阶时间更新格式,而普通的非分裂完全匹配层（PML）吸收边界局限于低阶时间格式。辅助微分方程完全匹配层（ADE PML）是一种可以适应任意阶时间格式的非分裂完全匹配层技术,且可以直接应用复频移拉伸算子以提高PML在高角度入射时的效果。作者将ADE PML应用于声波方程四阶Runge Kutta时间格式的数值模拟中,对其吸收效能进行了检验。数值模拟表明,复频移ADE PML在高角度入射时表现优于非复频移ADE PML。另外,不同辅助变量更新格式的吸收效果存在微小差异,显格式下计算结果与解析解吻合较好。长时间能量衰减计算表明ADE PML可以稳定至2 × 10⁵时间步。     

One step back for a leap forward: toward operational measurements of elements at risk

This study presents the evaluation of 1 year of operational lightning forecasts provided for Europe, using the Weather Research and Forecasting model coupled with a cloud-top height-based lightning parameterization scheme. Three different convective parameterization schemes were employed for parameterizing sub-grid cloud-top heights and consequently driving the lightning scheme. Triggering of the lightning scheme was controlled by means of a model-resolved microphysics-based masking filter, while the formulation for deriving lightning flash rates was also modified, assuming a single “marine” equation instead of the original equations discriminating between continental and marine lightning. Gridded lightning observations were used for evaluating model performance on a dichotomous decision basis. Analysis showed that the lightning scheme is sensitive to the parameterization of convection. In particular, the Kain–Fritsch convective scheme was found to outperform the Grell–Devenyi and Grell–Freitas schemes, showing a statistically significant better performance with respect to lightning prediction. This was most evident during the warm season, while smaller differences among the schemes were recorded during the cold season. Further, for all examined convective schemes, it was found that the application of the masking filter is desirable for improving model performance in terms of lightning forecasting. Last, the reported results revealed that the refinement of the formulation of the lightning parameterization scheme, adhering to a “global” marine equation instead of distinguishing between land and sea lightning, may be necessary in order to obtain reliable lightning forecasts.     

Numerical simulation of seismic waves in models with anisotropic formations: coupling Virieux and Lebedev finite-difference schemes

Anisotropy is widespread in the Earth’s interior. However, there is a number of models where anisotropic formations comprise as few as 10–20?% of the volume, and this includes fractured reservoirs, thin-layered packs, etc. while the major part of the medium is isotropic. In this situation, the use of computationally intense anisotropy-oriented approaches throughout the computational domain is prodigal. So this paper presents an original advanced finite-difference algorithm based on the domain decomposition technique with individual scheme used inside subdomains. It means that the standard staggered grid scheme or the Virieux scheme is used in the main part of the model which is isotropic, while the anisotropy-oriented Lebedev scheme is utilized inside domains with anisotropic formations. Finite-difference consistency conditions at the artificial interface where the schemes are coupled are designed to make the artificial reflections as low as possible, namely, for the second-order scheme, the third order of convergence of the reflection coefficients is proved.     

Explicit versus implicit front advancing schemes for the simulation of hydraulic fracture growth

The proper computation of the time evolution of the fracture front is the main challenge of three-dimensional (3D) hydraulic fracture growth simulation. We discuss explicit and implicit variants of a hydraulic fracture propagation scheme based on a level set representation of the fracture. Such a scheme couples a finite discretization of the governing equations and the near-tip hydraulic fracture asymptotes. We benchmark the accuracy, robustness, and stability of these different front advancing schemes on a number of test cases. Our results indicate a large computational gain of the explicit scheme at the expense of a slightly less accurate solution (few percent less accuracy over few time steps) when crossing heterogeneities. The predictor corrector scheme combines at least an approximately 25% computational gain while retaining the stability and accuracy of the fully implicit version of the scheme in all cases.     

二维分数阶对流-弥散方程的数值解

对二维时间分数阶对流-弥散方程和二维空间分数阶对流-弥散方程分别建立了差分格式,实现了对其的数值求解。针对理想算例进行计算求解,分析了时间和空间分数阶阶数取不同值时的扩散变化规律,验证了各自所描述的时间相关性与空间相关性。同时与传统的二维整数阶对流-弥散方程的求解结果作了对比。当时间和空间分数阶阶数α与γ分别取整数时,二维时间分数阶对流-弥散方程和二维空间分数阶对流-弥散方程都与传统二维整数阶对流-弥散方程的计算结果相同,说明提出的对二维分数阶对流-弥散方程的数值求解方法是可行的。其结果对地下水溶质运移的进一步研究提供了有效的手段。     

A Godunov method for localization in elastoplastic granular flow

The equations governing the elastic-plastic deformation of granular materials are typically hyperbolic, or contain small-magnitude damping or rate effects. A finite element algorithm is the standard method for the numerical integration of these systems. In particular, finite elements allow great flexibility in the design of grid geometry. However, modern finite difference methods for hyperbolic systems have been successful in aerodynamics computations, resolving wave structures more sharply than finite element schemes. In this paper we develop a finite difference scheme for granular flow problems. We report on a second-order Godunov-type scheme for the integration of hyperbolic equations for the elastoplastic deformation of a simple model of granular flow. The Godunov method includes a characteristic tracing step in the integration, providing minimal wave dispersion, and a slope limiting step, preventing unphysical oscillations. The granular flow model we consider is hyperbolic, but hyperbolicity is lost at a large value of accumulated plastic strain. This loss of hyperbolicity is a tell-tale signal for the formation of a shear band within the sample. Typically, when systems lose hyperbolicity a regularization mechanism is added to the model equations in order to maintain the well posedness of the system. These regularizations include viscosity, viscoplasticity, higher-order gradient effects or stress coupling. Here we appeal to a very different kind of regularization. When the system loses hyperbolicity and a shear band forms, we treat the band as an internal boundary, and impose jump conditions at this boundary. Away from the band, the system remains hyperbolic and the integration step proceeds as usual.     

A fully-coupled discontinuous Galerkin method for two-phase flow in porous media with discontinuous capillary pressure

In this paper, we formulate and test numerically a fully-coupled discontinuous Galerkin (DG) method for incompressible two-phase flow with discontinuous capillary pressure. The spatial discretization uses the symmetric interior penalty DG formulation with weighted averages and is based on a wetting-phase potential/capillary potential formulation of the two-phase flow system. After discretizing in time with diagonally implicit Runge-Kutta schemes, the resulting systems of nonlinear algebraic equations are solved with Newton’s method and the arising systems of linear equations are solved efficiently and in parallel with an algebraic multigrid method. The new scheme is investigated for various test problems from the literature and is also compared to a cell-centered finite volume scheme in terms of accuracy and time to solution. We find that the method is accurate, robust, and efficient. In particular, no postprocessing of the DG velocity field is necessary in contrast to results reported by several authors for decoupled schemes. Moreover, the solver scales well in parallel and three-dimensional problems with up to nearly 100 million degrees of freedom per time step have been computed on 1,000 processors.     

Hybrid Gaussian-cubic radial basis functions for scattered data interpolation

Scattered data interpolation schemes using kriging and radial basis functions (RBFs) have the advantage of being meshless and dimensional independent; however, for the datasets having insufficient observations, RBFs have the advantage over geostatistical methods as the latter requires variogram study and statistical expertise. Moreover, RBFs can be used for scattered data interpolation with very good convergence, which makes them desirable for shape function interpolation in meshless methods for numerical solution of partial differential equations. For interpolation of large datasets, however, RBFs in their usual form, lead to solving an ill-conditioned system of equations, for which, a small error in the data can cause a significantly large error in the interpolated solution. In order to reduce this limitation, we propose a hybrid kernel by using the conventional Gaussian and a shape parameter independent cubic kernel. Global particle swarm optimization method has been used to analyze the optimal values of the shape parameter as well as the weight coefficients controlling the Gaussian and the cubic part in the hybridization. Through a series of numerical tests, we demonstrate that such hybridization stabilizes the interpolation scheme by yielding a far superior implementation compared to those obtained by using only the Gaussian or cubic kernels. The proposed kernel maintains the accuracy and stability at small shape parameter as well as relatively large degrees of freedom, which exhibit its potential for scattered data interpolation and intrigues its application in global as well as local meshless methods for numerical solution of PDEs.     

Automatic substepping integration of the subloading tij model with stress path dependent hardening

This paper discusses the numerical integration of the subloading tij model. This is an elastoplastic model with stress path dependent hardening, which can predict the behaviour of normally consolidated clays or loose sands, as well as of over-consolidated clays or dense sands, with a small number of material parameters. Three features distinguish the subloading tij model from the conventional ones: (a) the use of a modified stress space given by tensor t_ij; (b) the split of the plastic strain increments in two components leading to a stress path dependent hardening; and (c) the use of two yield surfaces (subloading yield surface and normal yield surface). This last feature is based on the concept of sub-yielding stress states and adds an extra internal strain-like hardening variable, related to the relative density state, which demands its own evolution law. The three characteristics above greatly improve the prediction capabilities of the model, with respect to those of the well-known Cam clay model, at the cost of only two additional parameters. Nonetheless, the numerical integration of the constitutive equations of subloading tij model is a bit challenging, mainly due to the stress path dependent hardening. In order to integrate the equations of subloading tij model in the same way as for any conventional model, the authors reformulated its equations in a simpler and direct manner. Here, these equations are integrated using multi-step explicit schemes, such as modified-Euler and Runge–Kutta–Dormand–Price, with automatic error control. Simple forward-Euler scheme is also used for the sake of comparison. The results show that the modified-Euler scheme is more accurate as well as faster than the other schemes analysed over a wide range of error tolerance. Besides, the automatic feature of these schemes is a great convenience for the users of numerical codes.     

声波散射数值模拟的两种新方案

声波散射的数值模拟问题一般用网格法或积分方程法解决。当模型的尺度很大时,两种方法都会遇到计算机资源不足所造成的困难。另外,在网格法中,场源的位置和场源附近的波场奇异性逼近精度都受网格点的控制,因此难以满足实际问题所提出的要求。针对这些问题,提出了两种处理声波散射问题的新方案。一种主要针对网格法,另外一种针对积分方程法。在针对网格法的方案中,通过模型分解和波场分裂,将原始的总场计算问题转化为散射场计算问题。由于背景场是由解析公式给出的,所以可以将场源放置在数值网格的任意位置,不一定非得在网格点上。基于同样的原因,场源附近的波场奇异性可以精确地算出。在针对积分方程法的方案中,通过引入拟线性近似,使得散射场的数值求解不必再借助于代数方程组,只要进行数值积分即可。所建立的数值计算方案具有普遍的适用性,其基本思想可以直接用于解决弹性波散射的数值模拟问题并用于反演密度和速度。     

Simulation of severe thunder storm event: a case study over Pune,India

Numerical simulation of a typical tropical thunder storm event at Pune (18.53°N, 73.85°E), India, has been performed using the three nested domain configuration of Weather Research and Forecasting-Advanced Research Weather Model (version 3.2). The model simulations have been compared with observations. Sensitivity to cumulus parameterization schemes, namely Betts–Miller (BM), Grell–Devenyi (GD), and Kain–Fritsch (KF), for simulation of vertical structure and time evolution of weather parameters has been evaluated using observations from automatic weather station and global positioning system radiosonde ascents. Comparison of spatial distribution of 24-h accumulated rain with Tropical Rainfall Measuring Mission data shows that BM scheme could simulate better rain than GD and KF schemes. The BM scheme could well simulate the development of storm and heavy rain as it could generate sufficiently humid and deep layer in the lower and middle atmosphere, along with co-existence of updrafts and downdrafts and frozen hydrometeors at the middle level and rain water near the surface.     

Fast computation of Hankel Transform using orthonormal exponential approximation of complex kernel function

The computation of electromagnetic (EM) fields, for 1-D layered earth model, requires evaluation of Hankel Transform (HT) of the EM kernel function. The digital filtering is the most widely used technique to evaluate HT integrals. However, it has some obvious shortcomings. We present an alternative scheme, based on an orthonormal exponential approximation of the kernel function, for evaluating HT integrals. This approximation of the kernel function was chosen because the analytical solution of HT of an exponential function is readily available in literature. This expansion reduces the integral to a simple algebraic sum. The implementation of such a scheme requires that the weights and the exponents of the exponential function be estimated. The exponents were estimated through a guided search algorithm while the weights were obtained using Marquardt matrix inversion method. The algorithm was tested on analytical HT pairs available in literature. The results are compared with those obtained using the digital filtering technique with Anderson filters. The field curves for four types (A-, K-, H-and Q-type) of 3-layer earth models are generated using the present scheme and compared with the corresponding curves obtained using the Anderson sc heme. It is concluded that the present scheme is more accurate than the Anderson scheme     

Rapid solution of kriging equations,using a banded Gauss elimination algorithm

Summary A crucial concern when implementing computer algorithms for geostatistical analyses on microcomputers is speed of execution. Kriging, in particular, typically relies on a Gauss elimination algorithm to solve for weights. Because such an alogrithm is required for each estimate, the solution for weights can result in very slow program execution speed on a microcomputer. One approach to enhancing the efficiency of Gauss elimination is demonstrated herein. The upper triangle plus diagonal of the intersample covariance matrix is used in a modified banded Gauss elimination algorithm. Results show that such an algorithm yields approximately a two-fold reduction in execution time for kriging when the number of nearest neighbours used for estimation is large.     

A technique to test finite difference schemes to model some geophysical processes in a geological structure

A preliminary problem to solve in the set-up of a mathematical model simulating a geophysical process is the choice of a suitable discrete scheme to approximate the governing differential equations. This paper presents a simple technique to test finite difference schemes used in the modeling of geophysical processes occurring in a geological structure. This technique consists in generating analytical solutions similar to the ones characterizing a geophysical process, given general information on some relevant parameters. Useful information for the choice of the discrete scheme to employ in the mathematical model simulating the original geophysical process can be obtained from the comparison between the analytical solution and the approximated numerical solutions generated by means of different discrete schemes. Two classes of numerical examples approximating the differential equation that governs the steady state earth's heat flow have been treated using three different finite differences schemes. The first class of examples deals with media whose phenomenological parameters vary as continuous space functions; the second class, instead, deals with media whose phenomenological parameters vary as discontinuous space functions. The finite difference schemes that have been utilized are: Centered Finite Difference Scheme (CDS), Arithmetic Mean Scheme (AMS), and Harmonic Mean Scheme (HMS).The numerical simulations showed that: the CDS may yield physically inconsistent solutions if the lattice internodal distance is too large, but in case of phenomenological parameters varying as a continuous function, this pitfall can be avoided increasing the lattice node refinement. In case of phenomenological parameters varying as a discontinuous function, instead, the CDS may yield physically inconsistent solutions for any lattice-node refinement. The HMS produced good results for both classes of examples showing to be a scheme suitable to model situations like these.     

Exponential time integrators for stochastic partial differential equations in 3D reservoir simulation

The transport of chemically reactive solutes (e.g. surfactants, CO₂ or dissolved minerals) is of fundamental importance to a wide range of applications in oil and gas reservoirs such as enhanced oil recovery and mineral scale formation. In this work, we investigate exponential time integrators, in conjunction with an upwind weighted finite volume discretisation in space, for the efficient and accurate simulation of advection–dispersion processes including non-linear chemical reactions in highly heterogeneous 3D oil reservoirs. We model sub-grid fluctuations in transport velocities and uncertainty in the reaction term by writing the advection–dispersion–reaction equation as a stochastic partial differential equation with multiplicative noise. The exponential integrators are based on the variation of constants solution and solve the linear system exactly. While this is at the expense of computing the exponential of the stiff matrix representing the finite volume discretisation, the use of real Léja point or the Krylov subspace technique to approximate the exponential makes these methods competitive compared to standard finite difference-based time integrators. For the deterministic system, we investigate two exponential time integrators, the second-order accurate exponential Euler midpoint (EEM) scheme and exponential time differencing of order one (ETD1). All our numerical examples demonstrate that our methods can compete in terms of efficiency and accuracy compared with standard first-order semi-implicit time integrators when solving (stochastic) partial differential equations that model mixing and chemical reactions in 3D heterogeneous porous media. Our results suggest that exponential time integrators such as the ETD1 and EEM schemes could be applied to typical 3D reservoir models comprising tens to hundreds of thousands unknowns.     

China summer precipitation simulations using an optimal ensemble of cumulus schemes

RegCM3 (REGional Climate Model) simulations of precipitation in China in 1991 and 1998 are very sensitive to the cumulus parameterization. Among the four schemes available, none has superior skills over the whole of China, but each captures certain observed signals in distinct regions. The Grell scheme with the Fritsch-Chappell closure produces the smallest biases over the North; the Grell scheme with the Arakawa-Schubert closure performs the best over the southeast of 100°E; the Anthes-Kuo scheme is superior over the northeast; and the Emanuel scheme is more realistic over the southwest of 100°E and along the Yangtze River Basin. These differences indicate a strong degree of independence and complementarity between the parameterizations. As such, an ensemble is developed from the four schemes, whose relative contributions or weights are optimized locally to yield overall minimum root-mean-square errors from observed daily precipitation. The skill gain is evaluated by applying the identical distribution of the weights in a different period. It is shown that the ensemble always produces gross biases that are smaller than the individual schemes in both 1991 and 1998. The ensemble, however, cannot eliminate the large rainfall deficits over the southwest of 100°E and along the Yangtze River Basin that are systematic across all schemes. Further improve-ments can be made by a super-ensemble based on more cumulus schemes and/or multiple models.