一种二阶锥线性化方法在上限有限元中的应用研究

对于平面应变条件下岩土稳定性分析,基于线性规划的上限有限元需对常用的摩尔-库仑屈服准则形成的二阶锥约束进行线性化,直接地处理方法是以外接多边形替代锥体投影形成的圆域。为了提高线性化精度往往需直接增加外接多边形边数,从而造成线性规划模型中决策变量包含大量的塑性乘子变量,使计算难度大为增加甚至变得不可行。为此,引入Ben-Tal和Nemirovsky提出的一种二阶锥线性化方法,并将其嵌入到自编的上限有限元程序。经算例分析发现,该法与外接多边形线性化方法所获计算结果相互印证,且其通过适量的增加决策变量和等式约束数目,能保证摩尔-库仑屈服准则线性化精度,同时形成的线性规划规模更小,可望应用于基于线性规划模型的上限有限元中。     

Development of an infinite element boundary to model gravity for subsurface civil engineering applications

The accurate modelling of gravity is of crucial importance for a variety of issues including, but not restricted to, the identification of buried objects. Gravity is an unbounded problem, which causes challenges when applying numerical models, i..e.., it results in computational difficulties when specifying the relevant boundary conditions. In order to address this, previous research has tended to generate artificial boundary conditions, e.g., truncating the simulated domain and adding many unrealistic zero-density layers, which introduces more unknown parameters and unnecessarily excessive computational time. In order to overcome such inaccuracies, this paper proposes an innovative development of the finite element modelling technique, which represents a step change in the field of gravity forward modelling. A comprehensive formulation of an infinite element to reproduce the far-field boundary effect using only one layer of infinite elements is presented. The developed model considerably reduces the computational time while obtaining high degrees of accuracy. The model is validated against the exact solution of the problem, and its results show an excellent performance. The proposed method can significantly improve the postprocessing and interpretation stages of data analysis relevant to micro-gravity sensors. The new method is applied to subsurface civil engineering although its applicability is manifold.     

On the problem of boundary conditions in isentropic coordinates

The main advantages of constant potential enthalpy as a vertical coordinate are weaker horizontal velocity gradients in frontal regions and a higher vertical resolution. A disadvantage is the intersection of isentropes with the ground and folding of these surfaces. A numerical model is proposed to overcome the difficulties imposed by the intersection of isentropes with the ground. The model contains a physical and computational domain. The top and bottom surfaces of the computational domain are isentropes whereas the physical domain of the flow confined above by a free surface of constant pressure, and the bottom of this domain is the surface of the earth. In the present study the top surfaces of these two domains coincide (θ _T, P_Tare constants in space and time). The model was tested for the study of frontogenesis and cyclogenesis and integrated for 7 days. The results correspond to enstrophy-conserving finite difference scheme.     

Robust optimisation of CO2 sequestration strategies under geological uncertainty using adaptive sparse grid surrogates

Geologic CO₂ sequestration in deep saline aquifers is a promising technique to mitigate the effect of greenhouse gas emissions. Designing optimal CO₂ injection strategy becomes a challenging problem in the presence of geological uncertainty. We propose a surrogate assisted optimisation technique for robust optimisation of CO₂ injection strategies. The surrogate is built using Adaptive Sparse Grid Interpolation (ASGI) to accelerate the optimisation of CO₂ injection rates. The surrogate model is adaptively built with different numbers of evaluation points (simulation runs) in different dimensions to allow automatic refinement in the dimension where added resolution is needed. This technique is referred to as dimensional adaptivity and provides a good balance between the accuracy of the surrogate model and the number of simulation runs to save computational costs. For a robust design, we propose a utility function which comprises the statistical moment of the objective function. Numerical testing of the proposed approach applied to benchmark functions and reservoir models shows the efficiency of the method for the robust optimisation of CO₂ injection strategies under geological uncertainty.     

基于距离和局部Delaunay三角化控制的颗粒离散元模型填充方法研究

根据颗粒离散元Kelvin 接触力计算模型,分析了圆形颗粒体模拟材料力学特性应具备的条件,在此基础上提出了一种新颗粒模型构建方法。该方法首先在复杂模型域内随机生成种子,然后利用相切条件逐步扩展填充整个区域。填充过程中借助局部Delaunay三角化网格控制新颗粒的生成,采用复杂几何体距离控制颗粒与模型边界的相对位置,对靠近模型边界的颗粒进行容忍性优化填充,从而增加模型颗粒与边界的耦合性。同时对模型孔隙进行再填充,保证每个填充颗粒至少与3个颗粒相切,提高了模型内颗粒间的耦合性和模型的密度。最后采用任意多边形控制材料边界,将模型材料的设置简化为判断点是否在多边形内,简化了复杂模型材料属性的设置过程。结果表明：与膨胀颗粒生成法相比,该方法生成模型重叠量小、颗粒间及颗粒-边界相互耦合、填充率高。因此,颗粒黏结力破坏后不会造成飞溢现象,可适用于任意连通域模型的生成,能更好地实现复杂岩土细观介质变形破坏机制的模拟与研究。     

An approach for providing quasi‐convexity to yield functions and a generalized implicit integration scheme for isotropic constitutive models based on 2 unknowns

Yield and plastic potential surfaces are often affected by problems related to convexity. One such problem is encountered when the yield surface that bounds the elastic domain is itself convex; however, convexity is lost when the surface expands to pass through stress points outside the current elastic domain. In this paper, a technique is proposed, which effectively corrects this problem by providing linear homothetic expansion with respect to the centre of the yield surface. A very compact implicit integration scheme is also presented, which is of general applicability for isotropic constitutive models, provided that their yield and plastic potential functions are based on a separate mathematical definition of the meridional and deviatoric sections and that stress invariants are adopted as mechanical quantities. The elastic predictor‐plastic corrector algorithm is based on the solution of a system of 2 equations in 2 unknowns only. This further reduces to a single equation and unknown in the case of yield and plastic potential surfaces with a linear meridional section. The effectiveness of the proposed convexification technique and the efficiency and stability of the integration scheme are investigated by running numerical analyses of a notoriously demanding boundary value problem.     

Topic 3: State of the art of mathematical modelling of agricultural impact on groundwater pollution

The state of the art of mathematical modelling of groundwater pollution by agricultural activities is well advanced at the present time and a large number of models are available for a wide variety of water pollution problems. Also a variety of computational methods are available from very simple ones to the most advanced. In general, there are two ways of obtaining a water quality model: “from below,” when water quality data based on actual observation or measurement are explicitly analyzed by the calculus; and “from above,” when the calculus is available in advance as a theory and the user tries to find its interpretation for a given case and from rather limited data. Classifications may be based on the type of modelled processes, the mathematical approach, or on the purpose of the models. Eleven papers were presented in the Symposium on this topic, with most of the models being of the dynamic type.     

Shrinking Sphere Kinetics for Batch Dissolution of Mixed Particles of a Single Substance at High Under-Saturation: Validation with Sodium Chloride, but with Biogenic Silica in Mind

The cubic equation recently derived for the increase in concentration of a solute with time, as the solid dissolves in batch according to the shrinking sphere model at high under-saturation, is extended to dissolutions of mixtures of differently sized particles. This problem needs to be solved if batch dissolutions are to play their part in the proposed amelioration of global warming and associated climate change by accelerated ‘re-burial’ of excess CO₂ in ocean sediment. The upgraded model was tested using sodium chloride dissolved in 50% aqueous propanone, whence the model fitted two separate runs with 500 and 212 μm, and 212 and 38 μm, diameter crystals, respectively. The key to simulating dissolution in this way lies in the dissolutions being independent of each other. It is further shown that although this condition was implicit in the recent derivation of the cubic equation, it was not recognised at the time. The work should be applicable to any batch dissolution of mixed particles at high under-saturation, and hence, may find use in many industrial and laboratory dissolutions. Simulations show how agglomerated mixtures can yield a straight line on the plot of ln(1 − C/C _T) versus time, as was reported to occur recently with sodium chloride taken ‘straight from the bottle’. It is shown that this probably explains why exponential dissolutions may have seemed appropriate to the dissolution of biogenic silica in earlier literature. This study suggests that a new round of biogenic silica dissolutions, but with sieved samples, would be worthwhile, with the likelihood that shrinking sphere behaviour might well be found to characterise the kinetics. The opportunity is taken to investigate a number of aspects of the shrinking sphere model not generally discussed before, e.g. the graph for the change in surface area with time. The limitations of using cubic salt crystals with the shrinking sphere model are discussed.     

A quadrature-free method for simulation and inversion of 1.5D direct current (DC) borehole measurements

Resistivity inverse problems are routinely solved in order to characterize hydrocarbon bearing formations. They often require a large number of forward problems simulations. When considering a one dimensional (1D) planarly layered media, semi-analytical methods can be employed in order to solve a single forward problem in a fraction of a second. However, in some situations, a large number of (over one million) simulations is required, preventing this method to be used as a real time (logging) alternative. In this paper, we propose a novel semi-analytical method that dramatically reduces the total computational time, so it can be employed for real time inversion. In our proposed method, we select an ad hoc basis representation for the spectral solution such that its inverse Hankel transform can be computed analytically. The proposed method requires a pre-process that is expensive when compared with a single evaluation in classical semi-analytical methods. However, subsequent evaluations can be rapidly obtained, decreasing thus the total computational time by orders of magnitude when the number of required forward simulations is large.     

Discrete element method simulations of the seismic response of shallow foundations including soil‐foundation‐structure interaction

A novel three‐dimensional particle‐based technique utilizing the discrete element method is proposed to analyze the seismic response of soil‐foundation‐structure systems. The proposed approach is employed to investigate the response of a single‐degree‐of‐freedom structure on a square spread footing founded on a dry granular deposit. The soil is idealized as a collection of spherical particles using discrete element method. The spread footing is modeled as a rigid block composed of clumped particles, and its motion is described by the resultant forces and moments acting upon it. The structure is modeled as a column made of particles that are either clumped to idealize a rigid structure or bonded to simulate a flexible structure of prescribed stiffness. Analysis is done in a fully coupled scheme in time domain while taking into account the effects of soil nonlinear behavior, the possible separation between foundation base and soil caused by rocking, the possible sliding of the footing, and the dynamic soil‐foundation interaction as well as the dynamic characteristics of the superstructure. High fidelity computational simulations comprising about half a million particles were conducted to examine the ability of the proposed technique to model the response of soil‐foundation‐structure systems. The computational approach is able to capture essential dynamic response patterns. The cyclic moment–rotation relationships at the base center point of the footing showed degradation of rotational stiffness by increasing the level of strain. Permanent deformations under the foundation continued to accumulate with the increase in number of loading cycles. Copyright © 2011 John Wiley & Sons, Ltd.     

On numerical computation of impedance functions for rigid soil-structure interfaces embedded in heterogeneous half-spaces

We present fully-discrete procedures for computing the impedance functions of rigid massless soil-structure interfaces that are embedded in arbitrarily heterogeneous half-spaces. The finite element method (FEM) is used for obtaining the wave responses of (visco-)elastic half-spaces truncated by Perfectly Matched Layers (PMLs), which provide the wave absorbing boundary conditions. The devised FEM-PML approach is verified in both time and frequency domains by using various benchmark solutions. Requirements on the prescribed input excitations for obtaining accurate impedances in the time domain as well as the relative computational cost of time- and frequency domain solutions are investigated. Accuracy of the implemented PMLs in extracting the impedance functions is also examined in comparison to Lysmer–Kuhlemeyer dashpots; and it was found that this simplified boundary treatment is generally inadequate. The utility of the proposed method is demonstrated by extracting the impedance matrix of rectangular and circular voids embedded in a linearly stiffening half-space. Impedance functions for such complex soil-structure systems are shown to be highly coupled and frequency-dependent due to wave reflections and interference caused by the soil heterogeneity and interface geometry. Fully discrete approaches, such as those proposed herein, are necessary to adequately capture these effects.     

Semi‐analytical far field model for three‐dimensional finite‐element analysis

A challenging computational problem arises when a discrete structure (e.g. foundation) interacts with an unbounded medium (e.g. deep soil deposit), particularly if general loading conditions and non‐linear material behaviour is assumed. In this paper, a novel method for dealing with such a problem is formulated by combining conventional three‐dimensional finite‐elements with the recently developed scaled boundary finite‐element method. The scaled boundary finite‐element method is a semi‐analytical technique based on finite‐elements that obtains a symmetric stiffness matrix with respect to degrees of freedom on a discretized boundary. The method is particularly well suited to modelling unbounded domains as analytical solutions are found in a radial co‐ordinate direction, but, unlike the boundary‐element method, no complex fundamental solution is required. A technique for coupling the stiffness matrix of bounded three‐dimensional finite‐element domain with the stiffness matrix of the unbounded scaled boundary finite‐element domain, which uses a Fourier series to model the variation of displacement in the circumferential direction of the cylindrical co‐ordinate system, is described. The accuracy and computational efficiency of the new formulation is demonstrated through the linear elastic analysis of rigid circular and square footings. Copyright © 2004 John Wiley & Sons, Ltd.     

电法测井的三维有限元模拟

研究了不同测量方式(井-地,地-井,井-井)下点源场井中电法的三维有限元数值模拟。考虑到深度方向上大范围的网格剖分和井眼的影响及井-井测量等因素,采用放射状三棱柱单元的网格剖分方式,以提高网格质量,减少剖分单元数;给出了三棱柱单元的坐标变换公式,进行精确的单元积分,减少了单元积分时间;结合非结构化网格技术,实现了复杂模型的模拟;开发出相应的程序实现了复杂条件下(如考虑井眼影响、井井测量、倾斜井情形、地形起伏等)电法测井的三维有限元模拟,数值算例验证了方法的可靠性及计算效率,并对不同情形下的异常响应进行了分析,为进一步的反演工作奠定了基础。     

Adaptive POD model reduction for solute transport in heterogeneous porous media

We study the applicability of a model order reduction technique to the solution of transport of passive scalars in homogeneous and heterogeneous porous media. Transport dynamics are modeled through the advection-dispersion equation (ADE) and we employ Proper Orthogonal Decomposition (POD) as a strategy to reduce the computational burden associated with the numerical solution of the ADE. Our application of POD relies on solving the governing ADE for selected times, termed snapshots. The latter are then employed to achieve the desired model order reduction. We introduce a new technique, termed Snapshot Splitting Technique (SST), which allows enriching the dimension of the POD subspace and damping the temporal increase of the modeling error. Coupling SST with a modeling strategy based on alternating over diverse time scales the solution of the full numerical transport model to its reduced counterpart allows extending the benefit of POD over a prolonged temporal window so that the salient features of the process can be captured at a reduced computational cost. The selection of the time scales across which the solution of the full and reduced model are alternated is linked to the Péclet number (P e), representing the interplay between advective and dispersive processes taking place in the system. Thus, the method is adaptive in space and time across the heterogenous structure of the domain through the combined use of POD and SST and by way of alternating the solution of the full and reduced models. We find that the width of the time scale within which the POD-based reduced model solution provides accurate results tends to increase with decreasing P e. This suggests that the effects of local-scale dispersive processes facilitate the POD method to capture the salient features of the system dynamics embedded in the selected snapshots. Since the dimension of the reduced model is much lower than that of the full numerical model, the methodology we propose enables one to accurately simulate transport at a markedly reduced computational cost.     

Limit analysis of 2-D and 3-D structures based on an ellipsoid yield criterion

In this paper, a nonlinear numerical technique is developed to calculate the limit load and failure mode of structures obeying an ellipsoid yield criterion by means of the kinematic limit theorem, nonlinear programming theory and displacement-based finite element method. Using an associated flow rule, a general yield criterion expressed by an ellipsoid equation can be directly introduced into the kinematic theorem of limit analysis. The yield surface is not linearized and instead a nonlinear purely kinematic formulation is obtained. The nonlinear formulation has a smaller number of constraints and requires less computational effort than a linear formulation. By applying the finite element method, the kinematic limit analysis with an ellipsoid yield criterion is formulated as a nonlinear mathematical programming problem subject to only a small number of equality constraints. The objective function corresponds to the dissipation power which is to be minimized and an upper bound to the plastic limit load of a structure can then be calculated by solving the minimum optimization problem. An effective, direct iterative algorithm has been developed to solve the resulting nonlinear programming formulation. The calculation is based purely on kinematically admissible velocities. The stress field does not need to be calculated and the failure mode of structures can be obtained. The proposed method can be used to calculate the bearing capacity of clay soils in a direct way. Some examples are given to illustrate the validity and effectiveness of the proposed method.     

Lower bound limit analysis using finite elements and linear programming

This paper describes a technique for computing lower bound limit loads in soil mechanics under conditions of plane strain. In order to invoke the lower bound theorem of classical plasticity theory, a perfectly plastic soil model is assumed, which may be either purely cohesive or cohesive-frictional, together with an associated flow rule. Using a suitable linear approximation of the yield surface, the procedure computes a statically admissible stress field via finite elements and linear programming. The stress field is modelled using linear 3-noded traingles and statically admissible stress discontinuities may occur at the edges of each triangle. Imposition of the stress-boundary, equilibrium and yield conditions leads to an expression for the collapse load which is maximized subject to a set of linear constraints on the nodal stresses. Since all of the requirements for a statically admissible solution are satisfied exactly (except for small round-off errors in the optimization computations), the solution obtained is a strict lower bound on the true collapse load and is therefore ‘safe’. A major drawback of the technique, as first described by Lysmer,¹ is the large amount of computer time required to solve the linear programming problem. This paper shows that this limitation may be avoided by using an active set algorithm, rather than the traditional simplex or revised simplex strategies, to solve the resulting optimization problem. This is due to the nature of the constraint matrix, which is always very sparse and typically has many more rows that columns. It also proved that the procedure can, without modification, be used to derive strict lower bounds for a purely cohesive soil which has increasing strength with depth. This important class of problem is difficult to tackle using conventional methods. A number of examples are given to illustrate the effectiveness of the procedure.     

An improved implicit numerical integration of a non‐associated,three‐invariant cap plasticity model with mixed isotropic–kinematic hardening for geomaterials

The Sandia GeoModel is a continuum elastoplastic constitutive model that captures many features of the mechanical response for geological materials over a wide range of porosities and strain rates. Among the specific features incorporated into the formulation are a smooth compression cap, isotropic/kinematic hardening, nonlinear pressure dependence, strength differential effect, and rate sensitivity. This study attempts to provide enhancements regarding computational tractability, domain of applicability, and robustness of the model. A new functional form is presented for the yield and plastic potential functions. This reformulation renders a more accurate, robust, and efficient model as it eliminates spurious solutions attributed to the original form. In addition, we achieve a high‐performance implementation, because the local iterative method is allowed to recast residual vectors with a uniform dimensionality. The model is also furnished with a smooth, elliptical tension cap to account for the tensile yielding. Moreover, an efficient algorithm is introduced, which decreases the computational cost by differentiating the updated shear yield surface from the cap surface based on the trial relative stress state. Finally, various numerical examples including a large‐scale boundary value problem are presented to demonstrate the fidelity of the modified model and to analyze its numerical performance. Copyright © 2015 John Wiley & Sons, Ltd.     

Three-dimensional temporally dependent dispersion through porous media: analytical solution

A three-dimensional model for non-reactive solute transport in physically homogeneous subsurface porous media is presented. The model involves solution of the advection-dispersion equation, which additionally considered temporally dependent dispersion. The model also account for a uniform flow field, first-order decay which is inversely proportional to the dispersion coefficient and retardation factor. Porous media with semi-infinite domain is considered. Initially, the space domain is not solute free. Analytical solutions are obtained for uniform and varying pulse-type input source conditions. The governing solute transport equation is solved analytically by employing Laplace transformation technique (LTT). The solutions are illustrated and the behavior of solute transport may be observed for different values of retardation factor, for which simpler models that account for solute adsorption through a retardation factor may yield a misleading assessment of solute transport in ‘‘hydrologically sensitive’’ subsurface environments.     

Using global node-based velocity in random walk particle tracking in variably saturated porous media: application to contaminant leaching from road constructions

Precise and efficient numerical simulation of transport processes in subsurface systems is a prerequisite for many site investigation or remediation studies. Random walk particle tracking (RWPT) methods have been introduced in the past to overcome numerical difficulties when simulating propagation processes in porous media such as advection-dominated mass transport. Crucial for the precision of RWPT methods is the accuracy of the numerically calculated ground water velocity field. In this paper, a global node-based method for velocity calculation is used, which was originally proposed by Yeh (Water Resour Res 7:1216–1225, 1981). This method is improved in three ways: (1) extension to unstructured grids, (2) significant enhancement of computational efficiency, and (3) extension to saturated (groundwater) as well as unsaturated systems (soil water). The novel RWPT method is tested with numerical benchmark examples from the literature and used in two field scale applications of contaminant transport in saturated and unsaturated ground water. To evaluate advective transport of the model, the accuracy of the velocity field is demonstrated by comparing several published results of particle pathlines or streamlines. Given the chosen test problem, the global node-based velocity estimation is found to be as accurate as the CK method (Cordes and Kinzelbach in Water Resour Res 28(11):2903–2911, 1992) but less accurate than the mixed or mixed-hybrid finite element methods for flow in highly heterogeneous media. To evaluate advective–diffusive transport, a transport problem studied by Hassan and Mohamed (J Hydrol 275(3–4):242–260, 2003) is investigated here and evaluated using different numbers of particles. The results indicate that the number of particles required for the given problem is decreased using the proposed method by about two orders of magnitude without losing accuracy of the concentration contours as compared to the published numbers.