Solving density driven flow problems with efficient spatial discretizations and higher-order time integration methods

Modelling density driven flow problems requires an excessive computational time and/or heavy equipments due to the non-linear coupling between flow and transport equations. In this work, we develop a robust numerical model with efficient advanced approximations for both spatial and temporal discretizations in order to reduce the excessive computational requirement while maintaining accuracy.     

Non-linear dynamic earth dam–foundation interaction using a BE–FE method

A general, rigorous, coupled Boundary Element–Finite Element (BE–FE) formulation is presented for non-linear seismic soil–structure interaction in two dimensions. The BE–FE method is applied to investigate the inelastic response of earth dams to transient SV waves. The dam body, consisting of heterogeneous materials modelled with a simple non-linear hysteretic model, is discretized with finite elements, whereas the elastic half-space is discretized with boundary elements. The study focuses on the combined effects of the material non-linearity and foundation flexibility. The results show the significant effect of the foundation flexibility in reducing the response through radiation of energy. For excitations with peak ground accelerations from 0·2gto 0·6g, the crest acceleration amplification ranges from 2·5 to 1·4 and seems to be comparable with field observations and results from other studies. Deamplification increasing with strain is reported at the lower part of the dam. The method is computationally powerful and can be used for efficient non-linear analysis of complex soil–structure systems. The efficiency of the BE–FE method allows further improvements with incorporation of a more advanced constitutive model and consideration of the generation and dissipation of pore-water pressures during the earthquake. © 1998 John Wiley & Sons, Ltd.     

Behaviour of a hydrodynamic finite element model

A finite element model which solves the vertically integrated momentum and continuity equations is described. Linear triangular elements are used to describe the geometry and parameter variations. The Galerkin method of weighted residuals is employed to cast the equations in a form amenable to numerical solution. The model is based on a fully implicit formulation using finite differences for the temporal derivatives.Means of evaluating the non-linear terms of the governing equations are described, and model results are presented for a frictionless tidal channel. The example is chosen such that the non-linearities have a large influence on the solution, and as a result the linearization scheme significantly affects the model's behaviour.Suppression of the non-linear instabilities generated by the convective terms in the momentum equations is examined for the case of flow around a 180° bend. Both the imposition of artificially high roughness coefficients and the use of an effective eddy viscosity are examined in terms of their ability to damp the oscillations which arise for this example.Finally, model results are presented for a case study involving determination of remedial measures to improve flow conditions at a river outfall in Southern Ontario.     

Investigating the Behaviour of Two-Dimensional Finite Element Models of Compound Channel Flow

This paper describes results from a recent study concerning the numerical modelling of compound channel flow using two generalized two-dimensional finite element codes specifically adapted to floodplain studies: RMA-2 and TELEMAC-2D. By application to an 11 km reach of the River Culum, Devon, UK, simulations are developed to investigate the impact of numerical technique, mesh resolution and topographic parameterization on model results. The research is shown to raise a number of issues concerning the construction, calibration and validation of two-dimensional finite element models for this flow problem. © 1997 by John Wiley & Sons, Ltd.     

Estimates of flow resistance and eddy viscosity coefficients for 2D modelling on vegetated floodplains

The problem of quantifying the effects of flexible plants on flow resistance and eddy viscosity by vegetated floodplains is first addressed with a one‐dimensional (1D) approximation based upon the so‐called lateral distribution method. The estimates so obtained are then tested with two‐dimensional (2D) numerical simulations based on the full shallow water equations through the use of the computational code Telemac‐2D. Data obtained on a physical model of the Besòs River (Spain), whose floodplains were covered with plastic ornamental plants to mimic the effect of flexible vegetation, is used for the validation of the numerical results. Additionally, the values of flow resistance estimated numerically with the 1D and 2D simulations are compared with values obtained in a rectangular flume under flow conditions (slope, water depth and artificial lining) similar to those used on the reduced model. It is then established that as more physical mechanisms are included in the mathematical model used to study the problem, the ratio between the floodplain and the main channel flow resistance coefficient increases. The approach demonstrates that whenever enough flow data is available, the lateral distribution method delivers values of flow resistance and eddy viscosity which are highly consistent with 2D numerical modelling. This finding could mean considerable savings in the burdensome task of specifying flow resistance and turbulence dissipation values for 2D modelling of large compound channel systems. Copyright © 2004 John Wiley & Sons, Ltd.     

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.     

Hillslope hydrographs by the finite element method

Simulation models may be used to explore the implications of making specific assumptions about the nature of a real world system, and then to make predictions of the behaviour of that system under a set of naturally occurring conditions. It is important that understanding generated by the former should be gained before predictive use of the system model. This paper describes and uses a finite-element model of transient, partially saturated water flow within a hillslope soil mantle overlying an impermeable bedrock, to make an investigation into the effects of parameter variations and initial conditions on the hillslope hydrograph. The results clearly demonstrate that the response of the hillslope system to rainfall is highly non-linear and that the initial conditions, particularly in the unsaturated zone, are of paramount importance in governing the timing and magnitude of the hydrograph peak. Hillslope convergence appears as the dominant topographic parameter but the non-linearity of the response and the complex interdependence between the soil and topographic parameters restrict the possibility of further definite conclusions about the relative sensitivity of the simulated hillslope hydrograph to changes in these parameters.     

Weakly non-linear analysis of small-amplitude internal gravity waves forced by isolated topography

Asymptotic methods and numerical simulations are used to examine the evolution of an internal gravity wave packet comprising a continuous spectrum of horizontal wavenumbers and propagating upwards in a continuously stratified shear flow. In the multiple-scale framework for a horizontally localized wave packet generated by stratified flow over a localized mountain range with multiple peaks, there are in general two horizontal scales: the “fast” scale which is defined by the oscillations within the packet, i.e. the number of peaks, and the “slow scale” which is defined by the horizontal extent of the packet, i.e. the width of the mountain range. The focus here is on the specific case of an isolated mountain where the spectrum of horizontal wavenumbers is centred at zero and the multiple-scaling procedure is thus simplified by the absence of the fast spatial scale. The background flow is vertically sheared and critical-level interactions occur. The time frame within which non-linear critical-level effects become significant is determined by the magnitude of the non-linear terms in the governing equations. With the isolated mountain forcing this time frame is significantly longer than in the case of a multiple-peak mountain range forcing and it depends on the horizontal scale of the forcing, as well as on the amplitude. At leading-order, the non-linear asymptotic solution approaches a steady state in the outer region at late time, but the zero-wavenumber component of the solution continues to evolve with time in the vicinity of the critical level.     

Simulation of damage in RC frames with variable axial forces

The purpose of this paper is to describe how the methods and concepts of continuum damage and fracture mechanics can be applied to the modelling of the hysteretic behaviour of Reinforced Concrete (RC) frame members under variable axial loads. A frame member is considered as the assemblage of an elastic beam-column and two inelastic hinges, as in conventional nonlinear analysis of RC frames. As a result of the combination of damage mechanics and standard RC theory, a simplified model of damage is proposed and implemented as a finite element (FE). This new element can be used with any non-linear commercial FE program. The numerical simulation of several experiments, for which data were available in the literature, verifies the accuracy of the model. Copyright © 1999 John Wiley & Sons, Ltd.     

Effects of chemical reactions on density-dependent fluid flow: on the numerical formulation and the development of instabilities

A three-dimensional, reactive numerical flow model is developed that couples chemical reactions with density-dependent mass transport and fluid flow. The model includes equilibrium reactions for the aqueous species, kinetic reactions between the solid and aqueous phases, and full coupling of porosity and permeability changes that result from precipitation and dissolution reactions in porous media. A one-step, global implicit approach is used to solve the coupled flow, transport and reaction equations with a fully implicit upstream-weighted control volume discretization. The Newton–Raphson method is applied to the discretized non-linear equations and a block ILU-preconditioned CGSTAB method is used to solve the resulting Jacobian matrix equations. This approach permits the solution of the complete set of governing equations for both concentration and pressure simultaneously affected by chemical and physical processes. A series of chemical transport simulations are conducted to investigate coupled processes of reactive chemical transport and density-dependent flow and their subsequent impact on the development of preferential flow paths in porous media. The coupled effects of the processes driving flow and the chemical reactions occurring during solute transport is studied using a carbonate system in fully saturated porous media. Results demonstrate that instability development is sensitive to the initial perturbation caused by density differences between the solute plume and the ambient groundwater. If the initial perturbation is large, then it acts as a “trigger” in the flow system that causes instabilities to develop in a planar reaction front. When permeability changes occur due to dissolution reactions occurring in the porous media, a reactive feedback loop is created by calcite dissolution and the mixed convective transport of the system. Although the feedback loop does not have a significant impact on plume shape, complex concentration distributions develop as a result of the instabilities generated in the flow system.     

The unsteady plane flow of ice-sheets: A parabolic problem with two moving boundaries

Abstract

Finite difference algorithms have been developed to solve a one-dimensional non-linear parabolic equation with one or two moving boundaries and to analyse the unsteady plane flow of ice-sheets. They are designed to investigate the response of an ice-sheet to changes in climate, and to reconstruct climatic changes implied by past ice-sheet variations inferred from glacial geological data. Two algorithms are presented and compared. The first, a fixed domain method, replaces time as an independent variable with span. The grid interval in real space is kept constant, and thus the number of grid points changes with span. The second, a moving mesh method, retains time as one of the independent variables, but normalises the spatial variable relative to the span, which now enters the diffusion and advection coeficients in the parabolic equation for the surface profile.

Crank-Nicholson schemes for the solution of the equations are constructed, and iterative schemes for the solution of the resulting non-linear equations are considered.

Boundary (margin) motion is governed by the surface slope at the margin. Differentiation of the evolution equations results in an evolution equation for the margin slopes. It is shown that incorporation of this evolution equation, while not formally increasing the accuracy of the finite difference schemes, in practice increases accuracy of the solution.     

OpenGeoSys-Gem: A numerical tool for calculating geochemical and porosity changes in saturated and partially saturated media

Reactive transport codes that use a Gibbs Energy Minimization (GEM) to solve chemical equilibria are uncommon. We present a new coupling of the Richards flow module of the Finite Element (FE) based OpenGeoSys code with the GEM based chemical solver GEMS3K. The coupled code is highly parallelized using an overlapping domain decomposition approach in combination with execution of multiple threads that solve chemical equilibria in parallel. FE reactive transport schemes are often affected by spurious concentration oscillations. We effectively suppress these oscillations with a linearized algebraic flux corrected transport (FCT) algorithm. An application example is presented which investigates the evolution of material interfaces in a deep geological repository for nuclear waste. The example uses all features of the new coupled code: flow and multi-component transport in variably saturated media, and a very complex chemical setup which makes extensive use of (non-linear) solid solution formulations for mineral phases.     

A simple numerical experiment concerning the general circulation in the lower stratosphere

Summary By means of highly truncated spherical harmonic expansions, an extended four-level quasi-geostrophic model with variable Coriolis parameter is transformed into a set of ordinary non-linear differential equations. Non-adiabatic effects, frictional dissipation, and boundary effects are approximately included in the equations. A numerical experiment made with the equations succeeds in producing many realistic statistical gross features, especially in the lower stratosphere, e. g., a poleward temperature incrase, the up-gradient horizontal transports of heat and momentum due to large-scale eddies, the upward energy flux of extra-long waves, and the trapping of the upward energy flux of tropospheric unstable waves near the tropopause. The mean energy flow in the lower stratosphere and in the troposphere are analyzed and compared with each other, indicating very clearly the baroclinical activness of the troposphere and the passiveness of the lower stratosphere. The dynamics in the lower stratosphere are discussed. the mean meridional circulation is also studied.     

Predicting floodplain inundation: raster‐based modelling versus the finite‐element approach

We compare two approaches to modelling floodplain inundation: a raster‐based approach, which uses a relatively simple process representation, with channel flows being resolved separately from the floodplain using either a kinematic or diffusive wave approximation, and a finite‐element hydraulic model aiming to solve the full two‐dimensional shallow‐water equations. A flood event on a short (c. 4 km) reach of the upper River Thames in the UK is simulated, the models being validated against inundation extent as determined from satellite synthetic aperture radar (SAR) imagery. The unconstrained friction parameters are found through a calibration procedure, where a measure of fit between predicted and observed shorelines is maximized. The raster and finite‐element models offer similar levels of performance, both classifying approximately 84% of the model domain correctly, compared with 65% for a simple planar prediction of water surface elevation. Further discrimination between models is not possible given the errors in the validation data. The simple raster‐based model is shown to have considerable advantages in terms of producing a straightforward calibration process, and being robust with respect to channel specification. Copyright © 2001 John Wiley & Sons, Ltd.     

Free and steady state vibration of non-linear structures using a finite element–non-linear eigenvalue technique

The problem of free vibration of non-linear structures is considered initially. It is shown that this problem can be represented as a non-linear eigenvalue problem. Variational principles for non-linear eigenvalue problems are defined. These variational principles are implemented with finite element models to define numerical approximations for the free vibration problem. The solution of these approximate equations provides a set of non-linear modal vectors and natural frequencies which vary with the amplitude of the solution. The non-linear eigenvalue parameters can be used in modal expansion approximations for the non-linear transient or steady state response of structural systems. To demonstrate the proposed techniques the free vibration and steady state vibration characteristics of a geometrically non-linear circular plate are determined.     

Taylor columns in horizontally sheared flow

Abstract

Solutions of the steady, inviscid, non-linear equations for the conservation of potential vorticity are presented for linearly sheared geostrophic flow over a right circular cylinder. The indeterminancy introduced by the presence of closed streamline regions is removed by requiring that the steady flow retains above topography a given fraction of that fluid initially present there, assuming the flow to have been started from rest. Those solutions which retain the largest fraction in uniform and negatively sheared streams satisfy the Ingersoll (1969) criterion (that, in the limit of vanishingly small viscosity, closed streamline regions are stagnant) and so are unaffected by Ekman pumping. These flows are set up on the advection time scale. In positively sheared flows the maximum retention solutions do not satisfy the Ingersoll criterion and thus would be slowly spun down on the far longer viscous spin-up time.

For arbitrary isolated topography, both the partial retention and Ingersoll problems are reduced to a one-dimensional non-linear integral equation and the solution of the Ingersoll problem obtained in the limit of strong positive shear. The stagnant region is symmetric about the zero velocity line and extends to infinity in the streamwise direction. Its cross-stream width is proportional to the rotation rate and fractional height occupied by the obstacle and inversely proportional to the strength of the shear, decreasing inversely as the square of distance upstream and downstream.     

一维波动方程波阻抗反演的同伦方法

文中从地震勘探一维波动方程反问题出发,研究了一种反演地层参数的同伦方法,该方法把非线性方程组的求解转化成常微分方程初值问题的数值求解,从而给出一种稳定的计算速度快、抗噪能力强的全局收敛的反演方法．理论模型和实例试算的结果表明了同伦方法是一种有效的反演算法,特别适用于非线性、多极值的地球物理反演问题,在地球物理非线性反演中具有广泛的应用前景．     

Non-linear earthquake-response analysis of long-span cable-stayed bridges: Theory

The non-linear dynamic analysis of three-dimensional long-span cable-stayed bridges when subjected to seismic loading is formulated. All possible sources of non-linearity, such as cable sag, axial force-bending moment interaction in bridge towers and girders and change of geometry of the whole bridge due to large displacements are considered in the analysis. Both cases of uniform and multiple-support seismic excitations are considered in the non-linear formulation of the problem. A tangent stiffness, iterative procedure is utilized to capture the the non-linear seismic response. The non-linear equations of motion are solved using a step-by-step integration technique in the real displacement coordinate space as well as in the modal coordinate space to save computational time.     

A modal procedure for seismic analysis of non-linear base-isolated multistorey structures

A modal procedure for non-linear analysis of multistorey structures with high-damping base-isolation systems was proposed. Two different isolation devices were considered in the analysis: an high-damping laminated rubber bearing and a lead-rubber bearing. Starting from deformational properties verified by tests, the isolation systems were characterized using three different analytical models (an Elastic Viscous, a Bilinear Hysteretic and a Wen's Model) with parameters depending from maximum lateral strain. After non-linear modelling of isolation and lateral-force-resisting systems, the effects of material non-linearities were considered as pseudo-forces applied to the equivalent linear system (Pseudo-Force Method) and the formally linearized equations of motion were uncoupled by the transformation defined by the complex mode shapes. The modal responses were finally obtained with an extension of Nigam–Jennings technique to non-linear and non-classically damped systems, in conjunction with an iterative technique searching for non-linear contributions satisfying equations of motion and constitutive laws. Since the properties of the isolated structure usually change with maximun lateral strain of isolation bearings, the integration of a new set of governing equations was required for each design-displacement value. The procedure proposed was described in detail and then applied for the determination of modal and total seismic responses in some real cases. At first, a very good agreement between non-linear responses obtained with the proposed mode superposition and with a direct integration method was observed. Then a comparison of results obtained with the three different analytical models of the isolation bearings was carried out. At last, the exact modal response obtained with analytical models depending from the design displacement of the isolation bearings was compared with two different approximated solutions, evaluated using mode shapes and isolation properties, respectively, calculated under simplified hypothesis.© 1998 John Wiley & Sons, Ltd.     

Non-linear steady state vibration and dynamic snap through of shallow arch beams

The non-linear steady state vibration of shallow arch beams is studied by a finite element method based on the principle of virtual work. Both the free and forced periodic vibrations are considered. The axial and flexural deformations are coupled by the induced axial force along the beam element. The spatial discretization is achieved by the usual finite element method and the steady state nodal displacements are expanded into a Fourier series. The harmonic balance method gives a set of non-linear algebraic equations in terms of the vibrating frequency and the Fourier coefficients of nodal displacements. The non-linear algebraic equations are solved by the Newtonian algorithm iteratively. The combined algorithm is called the incremental harmonic balance method. The importance of the conditions of completeness and balanceability is presented. Since the non-linearity is essentially softening, different orders of internal resonances between two modes can occur repeatedly. Isolated response curves are possible and are connected to the bifurcation of a particular excited mode.