Asymptotic formulae for correcting finite element predicted natural frequencies of gravity and embankment dams

Extraction of natural frequencies of a gravity dam or an embankment dam plays an important role in the seismic design of the dam because the seismic response of a dam is dependent largely on the dynamic characteristics of the dam. Owing to the lack of exact solutions and the geometry of a dam, numerical methods such as finite element methods have been often used to extract the natural frequencies of the dam. Since the finite element method is an approximate one, the resulting finite element solution to the natural frequency of a dam cannot be safely used unless its accuracy is evaluated within the acceptable range for the seismic design of the dam. To solve this problem, some asymptotic formulae for correcting the finite element predicted natural frequencies of a gravity dam and an embankment dam have been developed in this paper. Since the present asymptotic formulae are derived from the fact that the finite element solution tends to the exact one if the finite element size used approaches zero, they provide a corrected solution of higher accuracy for the natural frequency of a dam so that the accuracy of a finite element solution can be evaluated against this corrected solution. After the correctness and usefulness of the present formulae are assessed, two practical examples have been given to show how the asymptotic formulae can be used to correct and evaluate the discretization error for the finite element predicted natural frequencies of gravity dams and embankment dams.     

Adaptive-grid simulation of groundwater flow in heterogeneous aquifers

The prediction of contaminant transport in porous media requires the computation of the flow velocity. This work presents a methodology for high-accuracy computation of flow in a heterogeneous isotropic formation, employing a dual-flow formulation and adaptive gridding. The dual equations, describing the hydraulic head and the streamfunction, are numerically solved through finite element approximations. The application of classic finite-element methods requires a rather large number of nodes to represent suitably the flow in high-contrast formations. We present a mesh-adaptive approach that enhances the accuracy of the numerical flow solution for a given computational effort. We rely on an a posteriori error estimator to identify areas where refinements of the finite element mesh are needed or unrefinements are acceptable. We also demonstrate through numerical experiments that the developed methodology efficiently enhances accuracy through successive mesh adaptation.     

Geoid Determination through Ellipsoidal Stokes Boundary-Value Problem

Martinec and Grafarend (1997) have shown how the construction of Green's function in the Stokes boundary-value problem with gravity data distributed on an ellipsoid of revolution is approached in the O(e ₀ ² )-approximation. They have also expressed the ellipsoidal Stokes function describing the effect of ellipticity of the boundary as a finite sum of elementary functions. We present an effective method of avoiding the singularity of spherical and the ellipsoidal Stokes functions, and also an analytical expression for the ellipsoidal Stokes integral around the computational point suitable for numerical solution. We give the numerical results of solving the ellipsoidal Stokes boundary-value problem and their difference with respect to the spherical Stoke boundary-value problem.     

A new finite element technique for the solution of two-phase flow through porous media

A new upstream weighting finite element technique is developed for improved solution of the two-phase immiscible flow equations. Unlike the upstream weighting technique used by previous investigators, the new technique does not employ finite difference concepts to achieve the required upstream weighting of relative permeabilities or mobilities. Instead, upstream weighting is achieved by (1) representing the relative permeabilities or mobilities as continuous functions expressed in terms of the shape functions and nodal values (2) using asymmetric weighting functions to weight the spatial terms in the flow equations. These weighting functions are constructed such that they are dependent on the flow direction along each side of an element.In conjunction with the proposed technique, two solution schemes for treating the resulting set of non-linear algebraic equations are presented. These are the fully-implicit chord slope incremental solution scheme and the Newton-Raphson solution scheme. Both schemes allow the use of large time steps without being unstable.The proposed numerical technique is applied to two problems (1) the one-dimensional Buckley-Leverett problem (2) the two-dimensional five-spot well flow problem. Results indicate that this technique is superior to not only earlier finite element schemes but also five-point upstream finite difference formulae.     

Stability and accuracy of two-dimensional kinematic wave overland flow modeling

A two-dimensional finite element based overland flow model was developed and used to study the accuracy and stability of three numerical schemes and watershed parameter aggregation error. The conventional consistent finite element scheme results in oscillations for certain time step ranges. The lumped and the upwind finite element schemes are tested as alternatives to the consistent scheme. The upwind scheme did not improve on the stability or the accuracy of the solution, while the lumped scheme provided stable and accurate solutions for time steps twice the size of time steps needed for the consistent scheme. A new accuracy based dynamic time step estimate for the two-dimensional overland flow kinematic wave solution is developed for the lumped scheme. The newly developed dynamic time step estimates are functions of the mesh size, and time of concentration of the watershed hydrograph. Due to lack of analytical solutions, the time step was developed by comparing numerical solutions of various levels of discretization to a reference solution using a very fine mesh and a very small time step. The time step criteria were tested on a different set of problems and proved to be adequate for accurate and stable solutions. A sensitivity analysis for the watershed slope, Manning’s roughness coefficient and excess rainfall rate was conducted in order to test the effect of parameter aggregation on the stability and accuracy of the solution. The results of this analysis show that aggregation of the slope data resulted in the highest error. The roughness coefficient had a smaller effect on the solution while the rainfall intensity did not show any significant effect on the flow rate solution for the range of rainfall intensity used. This work pioneers the challenge of providing guidelines for accurate and stable numerical solutions of the two-dimensional kinematic wave equations for overland flow.     

Error analysis of finite difference methods for two-dimensional advection–dispersion–reaction equation

In this paper, the numerical errors associated with the finite difference solutions of two-dimensional advection–dispersion equation with linear sorption are obtained from a Taylor analysis and are removed from numerical solution. The error expressions are based on a general form of the corresponding difference equation. The variation of these numerical truncation errors is presented as a function of Peclet and Courant numbers in X and Y direction, a Sink/Source dimensionless number and new form of Peclet and Courant numbers in X–Y plane. It is shown that the Crank–Nicolson method is the most accurate scheme based on the truncation error analysis. The effects of these truncation errors on the numerical solution of a two-dimensional advection–dispersion equation with a first-order reaction or degradation are demonstrated by comparison with an analytical solution for predicting contaminant plume distribution in uniform flow field. Considering computational efficiency, an alternating direction implicit method is used for the numerical solution of governing equation. The results show that removing these errors improves numerical result and reduces differences between numerical and analytical solution.     

Grid cell distortion and MODFLOW's integrated finite-difference numerical solution

The ground water flow model MODFLOW inherently implements a nongeneralized integrated finite-difference (IFD) numerical scheme. The IFD numerical scheme allows for construction of finite-difference model grids with curvilinear (piecewise linear) rows. The resulting grid comprises model cells in the shape of trapezoids and is distorted in comparison to a traditional MODFLOW finite-difference grid. A version of MODFLOW-88 (herein referred to as MODFLOW IFD) with the code adapted to make the one-dimensional DELR and DELC arrays two dimensional, so that equivalent conductance between distorted grid cells can be calculated, is described. MODFLOW IFD is used to inspect the sensitivity of the numerical head and velocity solutions to the level of distortion in trapezoidal grid cells within a converging radial flow domain. A test problem designed for the analysis implements a grid oriented such that flow is parallel to columns with converging widths. The sensitivity analysis demonstrates MODFLOW IFD's capacity to numerically derive a head solution and resulting intercell volumetric flow when the internal calculation of equivalent conductance accounts for the distortion of the grid cells. The sensitivity of the velocity solution to grid cell distortion indicates criteria for distorted grid design. In the radial flow test problem described, the numerical head solution is not sensitive to grid cell distortion. The accuracy of the velocity solution is sensitive to cell distortion with error <1% if the angle between the nonparallel sides of trapezoidal cells is <12.5 degrees. The error of the velocity solution is related to the degree to which the spatial discretization of a curve is approximated with piecewise linear segments. Curvilinear finite-difference grid construction adds versatility to spatial discretization of the flow domain. MODFLOW-88's inherent IFD numerical scheme and the test problem results imply that more recent versions of MODFLOW 2000, with minor modifications, have the potential to make use of a curvilinear grid.     

Effects of the geometry of two‐dimensional fractures on their hydraulic aperture and on the validity of the local cubic law

Flow through rough fractures is investigated numerically in order to assess the validity of the local cubic law for different fracture geometries. Two‐dimensional channels with sinusoidal walls having different geometrical properties defined by the aperture, the amplitude, and the wavelength of the walls' corrugations, the corrugations asymmetry, and the phase shift between the two walls are considered to represent different fracture geometries. First, it is analytically shown that the hydraulic aperture clearly deviates from the mean aperture when the walls' roughness, the phase shift, and/or the asymmetry between the fracture walls are relatively high. The continuity and the Navier–Stokes equations are then solved by means of the finite element method and the numerical solutions compared to the theoretical predictions of the local cubic law. Reynolds numbers ranging from 0.066 to 66.66 are investigated so as to focus more particularly on the effect of flow inertial effects on the validity of the local cubic law. For low Reynolds number, typically less than 15, the local cubic law properly describes the fracture flow, especially when the fracture walls have small corrugation amplitudes. For Reynolds numbers higher than 15, the local cubic law is valid under the conditions that the fracture presents a low aspect ratio, small corrugation amplitudes, and a moderate phase lag between its walls.     

Free surface and saturated-unsaturated analyses of borehole infiltration tests above the water table

Constant head borehole infiltration tests are widely used for the in situ evaluation of saturated hydraulic conductivities of unsaturated soils above the water table. The formulae employed in analysing the results of such tests disregard the fact that some of the infiltrating water may flow under unsaturated conditions. Instead, these formulae are based on various approximations of the classical free surface theory which treats the flow region as if it were fully saturated and enclosed within a distinct envelope, the so-called ‘free surface’. A finite element model capable of solving free surface problems is used to examine the mathematical accuracy of the borehole infiltration formulae. The results show that in the hypothetical case where unsaturated flow does not exist, the approximate formulae are reasonably accurate within·a practical range of borehole conditions. To see what happens under conditions closer to those actually encountered in the field, the effect of unsaturated flow on borehole infiltration is investigated by means of two different numerical models: a mixed explicit-implicit finite element model, and a mixed explicit-implicit integrated finite difference model. Both of these models give nearly identical results; however, the integrated finite difference model is considerably faster than the finite element model. The relatively low computational efficiency of the finite element scheme is attributed to the large number of operations required in order to re-evaluate the conductivity (stiffness) matrix at each iteration in this highly non-linear saturated-unsaturated flow problem. The saturated-unsaturated analysis demonstrates that the classical free surface approach provides a distorted picture of the flow pattern in the soil. Contrary to what one would expect on the basis of this theory, only a finite region of the soil in the immediate vicinity of the borehole is saturated, whereas a significant percentage of the flow takes place under unsaturated conditions. As a consequence of disregarding unsaturated flow, the available formulae may underestimate the saturated hydraulic conductivity of fine grained soils by a factor of two, three, or more. Our saturated-unsaturated analysis leads to an improved design of borehole infiltration tests and a more accurate method for interpreting the results of such tests. The analysis also shows how one can predict the steady state rate of infiltration from data collected during the early transient period of the test.     

Nonlinear seismic response analysis of earth dams

The objective of this paper is to propose a general and efficient numerical procedure for analysing the dynamic response of geotechnical structures, which are considered as both nonlinear and two phase systems. In Section 2, the appropriate coupled dynamic field equations for the response of a two-phase soil system are briefly reviewed. The finite element spatial discretization of the field equations is described and time integration for the resulting nonlinear semi-discrete finite element equations is discussed. In Section 3, iterative techniques are examined for the solution of the global nonlinear system of finite element equations. A large amount of computational effort is expended in the iterative phase of the solution and so the iterative procedure used must be both reliable and efficient. The performance of three iterative procedure is examined: Newton Raphson, Modified Newton Raphson and Quasi-Newton methods, including BGFS and Broyden updates. Finally, in Section 4, the elasto-plastic earthquake response analysis of a two phase nonhomogeneous earth dam is presented. Extensive documentation exists¹ for the particular problem selected including recorded earthquake motions at the base and crest of the dam. The results of the numerical calculations are compared to the recorded response of the dam.     

A two-grid method for coupled free flow with porous media flow

This paper presents a two-grid method for solving systems of partial differential equations modelling incompressible free flow coupled with porous media flow. This work considers both the coupled Stokes and Darcy as well as the coupled Navier-Stokes and Darcy problems. The numerical schemes proposed are based on combinations of the continuous finite element method and the discontinuous Galerkin method. Numerical errors and convergence rates for solutions obtained from the two-grid method are presented. CPU times for the two-grid algorithm are shown to be significantly less than those obtained by solving the fully coupled problem.     

Coupling of boundary and finite elements for soil-structure interaction problems

The investigation of complex soil-structure interaction problems is usually carried out with numerical solution procedures such as the finite element or the boundary element method. It must be noted, however, that the choice of one or the other of these approaches is not just a matter of preferences; depending on the type of the problem under consideration, either boundary or finite elements may be more advantageous. A considerable expansion in the computational power can be obtained, on the other hand, if one resorts to hybrid schemes which retain the main advantages of the two methods and eliminate their respective disadvantages. This paper presents results obtained with a boundary element-finite element coupling procedure, and discusses its applicability to some representative soil-structure interaction problems. The structures considered are elastic systems, such as foundations, tunnels and filled trenches (modelled by finite elements), which are coupled with homogeneous elastic halfspaces (modelled by boundary elements). The examples demonstrate the importance of using a model that includes wave radiation effects. The coupling approach is formulated entirely in the time domain so that an extension of the algorithm to non-linear analyses seems to present no further difficulties.     

Finite element analysis of three-dimensional groundwater flow problems

A numerical method is presented for analysing either steady state or transient three-dimensional groundwater flow problems. The governing equation is formulated in terms of the finite element process using the Galerkin approach, and cubic isoparametric elements are used to simulate the flow domain as these permit accurate modelling of curved boundaries. Particular attention is paid to the time dependent movement of the phreatic surface where an iterative technique based on the replacement of the original transient problem by a discrete number of steady state problems is used to effect a solution. Furthermore, in tracing the movement of the surface use is made of the element formulation theory in order to compute the normal to the boundary.The validity of the technique is first established by analysing a radially symmetrical problem for which an alternative analytical solution is available. Finally, a general three-dimensional flow system is studied for which there is no known analytical solution. It is shown that relatively few elements are required to yield practical solutions.     

面向目标自适应海洋可控源电磁三维矢量有限元正演

本文实现了一种面向目标自适应海洋可控源电磁三维矢量有限元方法.为满足三维复杂电性结构模拟的需求,网格剖分采用非结构化六面体.在组装刚度矩阵之后,形成的大型复数线性方程组分解为等价的实数形式,利用带预条件的广义最小残差法进行求解.在获得微分方程的解之后,为提高解的准确性,通过面向目标的自适应误差估计来指示网格细化,重点加密能使观测点数值模拟精度提高的网格.对于大规模三维数据,为了使模型空间的并行计算达到均衡负载的效果,我们使用METIS函数库来进行网格计算任务量的划分.最后,通过对比一维解析解与三维自适应矢量有限元计算结果,验证了程序的正确性;通过自适应过程中误差指示子的分布,验证了面向目标自适应的有效性;通过对三维复杂模型进行均衡负载下的并行计算,测试了程序的可扩展性.     

电导率各向异性的海洋电磁三维有限单元法正演

本文提出了一种基于非结构化网格的海洋电磁有限单元正演算法.为了回避场源奇异性,文中选用二次场算法,将背景电阻率设置为水平层状且各向异性,场源在水平层状各向异性介质中所激发的一次场通过汉克尔积分得到.基于Coulomb规范得到二次矢量位和标量位所满足的Maxwell方程组,通过Galerkin加权余量法形成大型稀疏有限元方程,采用不完全LU分解(ILU)预条件因子的quasi-minimum residual(QMR)迭代解法对有限元方程进行求解得到二次矢量位和标量位;进而,利用滑动平均方法得到二次矢量位和标量位在空间的导数,由此得到二次电磁场;通过一维模型对算法的可靠性进行验证,与此同时,针对实际复杂海洋电磁模型,比较有限元模拟结果与积分方程模拟结果,进一步验证算法精度.若干计算结果均表明,文中算法具有良好的通用性,适用于井中电磁、航空电磁,环境地球物理等非均匀且各向异性介质中的电磁感应基础研究.     

基于Lamb问题的谱元法和有限元法模拟精度比较

在地震动数值模拟方法中,谱元法和有限元法是应用较广泛的两种方法。基于经典的Lamb问题模型,首先推导给出地表竖向位移的解析解答。然后分别利用常用的四阶谱元法和线性有限元法,模拟了地表脉冲力源作用下模型的位移响应。考虑有意义的最短波长内的采样点个数及单元高宽比的变化,对比了两种方法的模拟精度;结果表明:对于谱元法,观测点与波源之间需至少包含两个网格,在此条件下,最短波长内包含一个网格（最短波长内5个采样点）时,数值解与解析解的误差小于1%,已达很高的精度;对于有限元法,最短波长内需包含10个网格时才能达到这一精度。此外,在满足网格尺寸要求的前提下,单元水平向与垂直向尺寸的比值在1∶1到5∶1的范围内时,谱元法和有限元法的模拟精度均变化不大。因此,单位波长内采样点个数相同时,谱元法的模拟精度比有限元法高的多,同时,在一定范围内两种方法的模拟结果对于宽高比的变化不敏感。     

An efficient probabilistic finite element method for stochastic groundwater flow

We present an efficient numerical method for solving stochastic porous media flow problems. Single-phase flow with a random conductivity field is considered in a standard first-order perturbation expansion framework. The numerical scheme, based on finite element techniques, is computationally more efficient than traditional approaches because one can work with a much coarser finite element mesh. This is achieved by avoiding the common finite element representation of the conductivity field. Computations with the random conductivity field only arise in integrals of the log conductivity covariance function. The method is demonstrated in several two- and three-dimensional flow situations and compared to analytical solutions and Monte Carlo simulations. Provided that the integrals involving the covariance of the log conductivity are computed by higher-order Gaussian quadrature rules, excellent results can be obtained with characteristic element sizes equal to about five correlation lengths of the log conductivity field. Investigations of the validity of the proposed first-order method are performed by comparing nonlinear Monte Carlo results with linear solutions. In box-shaped domains the log conductivity standard deviation σ_Y may be as large as 1.5, while the head variance is considerably influenced by nonlinear effects as σ_Y approaches unity in more general domains.     

Numerically quantifying energy loss caused by squirt flow

In interconnected microcracks, or in microcracks connected to spherical pores, the deformation associated with the passage of mechanical waves can induce fluid flow parallel to the crack walls, which is known as squirt flow. This phenomenon can also occur at larger scales in hydraulically interconnected mesoscopic cracks or fractures. The associated viscous friction causes the waves to experience attenuation and velocity dispersion. We present a simple hydromechanical numerical scheme, based on the interface-coupled Lamé–Navier and Navier–Stokes equations, to simulate squirt flow in the frequency domain. The linearized, quasi-static Navier–Stokes equations describe the laminar flow of a compressible viscous fluid in conduits embedded in a linear elastic solid background described by the quasi-static Lamé–Navier equations. Assuming that the heterogeneous model behaves effectively like a homogeneous viscoelastic medium at a larger spatial scale, the resulting attenuation and stiffness modulus dispersion are computed from spatial averages of the complex-valued, frequency-dependent stress and strain fields. An energy-based approach is implemented to calculate the local contributions to attenuation that, when integrated over the entire model, yield results that are identical to those based on the viscoelastic assumption. In addition to thus validating this assumption, the energy-based approach allows for analyses of the spatial dissipation patterns in squirt flow models. We perform simulations for a series of numerical models to illustrate the viability and versatility of the proposed method. For a 3D model consisting of a spherical crack embedded in a solid background, the characteristic frequency of the resulting P-wave attenuation agrees with that of a corresponding analytical solution, indicating that the dissipative viscous flow problem is appropriately handled in our numerical solution of the linearized, quasi-static Navier–Stokes equations. For 2D models containing either interconnected cracks or cracks connected to a circular pore, the results are compared with those based on Biot's poroelastic equations of consolidation, which are solved through an equivalent approach. Overall, our numerical simulations and the associated analyses demonstrate the suitability of the coupled Lamé–Navier and Navier–Stokes equations and of Biot's equations for quantifying attenuation and dispersion for a range of squirt flow scenarios. These analyses also allow for delineating numerical and physical limitations associated with each set of equations.     

Mixing cell method for solving the solute transport equation with spatially variable coefficients

The advection–dispersion equation with spatially variable coefficients does not have an exact analytical solution and is therefore solved numerically. However, solutions obtained with several of the traditional finite difference or finite element techniques typically exhibit spurious oscillation or numerical dispersion when advection is dominant. The mixing cell and semi-analytical solution methods proposed in this study avoid such oscillation or numerical dispersion when advection dominates. Both the mixing cell and semi-analytical solution methods calculate the spatial step size by equating numerical dispersion to physical dispersion. Because of the spatial variability of the coefficients the spatial step size varies in space. When the time step size Δt→0, the mixing cell method reduces to the semi-analytical solution method. The results of application to two cases show that the mixing cell and semi-analytical solution methods are better than a finite difference method used in the study. © 1998 John Wiley & Sons, Ltd.     

三角网格有限元法声波与弹性波模拟频散分析

本文对声波与弹性波方程进行有限元法离散,构造有限元法频散关系的一般特征值问题,分析了时间离散格式为中心差分的三角网格有限元法声波与弹性波模拟的频散特性. 比较了三种质量矩阵即分布式质量矩阵、集中质量矩阵和混合质量矩阵对有限元法频散的影响;选取四种典型三角网格,分析了混合质量矩阵有限元（MFEM）频散的方向各向异性;数值频散、方向各向异性随插值阶数的增加逐渐减弱,当空间为三阶插值时,频散主要表现为随采样率的变化而几乎无明显方向各向异性, 其频散幅值也较小. 控制其他影响因素不变的情况下,研究了不同波速比介质中弹性波的数值频散. 最后给出了三角网格MFEM的数值耗散性.