Convergence Difficulties Research Articles

Contribution to Numerical Modeling of Water Flow in Variably Saturated, Heterogeneous Porous Media

Solutions to the Richards equation for water flow in variably saturated porous media are the focus of this paper. Working with field conditions, the extreme variability and complexity of soil, initial and boundary conditions can make the flow problem difficult to solve. This paper proposes to improve the computational efficiency of the mixed hybrid finite element (MHFE) method coupled with the variables transformation. The transform variables were introduced in order to simulate problems with convergence difficulty attributed to the presence of sharp wetting fronts. Furthermore, for better convergence behaviour, a technique that switches between the mixed-form and the pressure-head based form of the Richard’s equation was applied. Special attention was given to the top boundary conditions dealing with ponding or evaporation problems. In order to avoid non-physical oscillation problems, a mass condensation scheme was implemented in the model. Performance indicators in time and error of different options of the numerical model are defined, analyzed and classified. Thus, for each test case, a suitable numerical method that identifies which form of the Richards equation is best suited, the relevance of the switching technique as well as the utility of the transformation of the primary variable is possible. The results for the 1D Numerical test cases that have been performed matched those from the literature results. For evaporation and infiltration problem’s, the number of iterations needed to get the solution decrease when using the method of transformed pressure. Finally, knowing the soil heterogeneity, initial and boundary conditions, an agglomerative hierarchical clustering allows to analyze the need or not to transform variables and to use other options.

Simulations of two-dimensional steady isothermal and non-isothermal steady flows with slip for a viscoelastic memory-integral fluid

Purpose – The purpose of this paper is to compute flow effects of the transition from adherence-to-slip in two-dimensional flows, for a polymer melt obeying a memory-integral viscoelastic equation, in isothermal and non-isothermal cases. Design/methodology/approach – Temperature dependence is expressed by Arrhenius and William-Landel-Ferry models. A coupling approach is defined. For the dynamic equations, the Stream-Tube Method (STM) is used with finite differences in a mapped rectangular domain of the real domain, where streamlines are parallel and straight. STM avoids particle-tracking problems and allows simple formulae to evaluate stresses resulting from the constitutive equation. For the temperature field, a finite-element method is carried out to solve the energy equation in the real domain. Findings – The approach avoids numerical problems arising with classical formulations and proves to be robust and efficient. Large elasticity levels are attained without convergence and refinement difficulties that may arise close to the “stick-slip” transition section. The method highlights the role of temperature conditions and reveals interesting differences for the ducts considered. Practical implications – The results of the study are of interest for polymer processing where slip at the wall can be encountered, in relation with the physical properties of the materials. Originality/value – The paper presents a simple approach that limits considerably numerical problems coming from stick-slip boundary conditions and avoids particle-tracking. Results are obtained at flow rates encountered in industrial conditions.

A three dimensional integral equation approach for fluids under confinement: Argon in zeolites.

In this work, we explore the ability of an inhomogeneous integral equation approach to provide a full three dimensional description of simple fluids under conditions of confinement in porous media. Explicitly, we will consider the case of argon adsorbed into silicalite-1, silicalite-2, and an all-silica analogue of faujasite, with a porous structure composed of linear (and zig-zag in the case of silicalite-1) channels of 5-8 Å diameter. The equation is based on the three dimensional Ornstein-Zernike approximation proposed by Beglov and Roux [J. Chem. Phys. 103, 360 (1995)] in combination with the use of an approximate fluid-fluid direct correlation function furnished by the replica Ornstein-Zernike equation with a hypernetted chain closure. Comparison with the results of grand canonical Monte Carlo/molecular dynamics simulations evidences that the theory provides an accurate description for the three dimensional density distribution of the adsorbed fluid, both at the level of density profiles and bidimensional density maps across representative sections of the porous material. In the case of very tight confinement (silicalite-1 and silicalite-2), solutions at low temperatures could not be found due to convergence difficulties, but for faujasite, which presents substantially larger channels, temperatures as low as 77 K are accessible to the integral equation. The overall results indicate that the theoretical approximation can be an excellent tool to characterize the microscopic adsorption behavior of porous materials.

Modeling and simulation of a frictional translational joint with a flexible slider and clearance

A finite element method for the dynamic description of multibody systems with a frictional translational joint is presented and discussed in this work. By considering the slider deformation and the clearance size between the slider and the guide, Coulomb and Stribeck models are used to describe friction forces. The geometric constraints of the translational joints are treated as multi-unilateral constraints, and the impacts between the slider and guide are neglected when the clearance sizes are very small. A model of contact forces is established based on the penalty method and a criterion is proposed to evaluate the influence of clearance size on nodal contact forces. Time-varying contact states of nodes are determined by applying two contact detection algorithms for stick and slip, respectively. Based on kineto-elasto-dynamic method and varying penalty factors, the method avoids convergence difficulties and enforces the accuracy. Numerical results reveal the capability of the new method and the influences of small deformations of the slider, clearance size and friction forces on the dynamic response of the multibody systems with a translational joint.

Numerical analysis of single pad of thrust bearing with the energy equation solved by the characteristic-based split method

The solution of the energy equation of thermo-elasto-hydrodynamic analysis for bearings by the finite element method usually leads to convergence difficulties due to the presence of convection terms inherited from the Navier–Stokes equations. In this work, the numerical analysis is performed with finite element method universally by adopting the characteristic-based split method to solve the energy equation. Five case studies of fixed pad thrust bearings have been set up with different geometries, loads, and lubricants. The two-dimensional film pressure is obtained by solving the Reynolds equation with pre-defined axial load on the pad. The energy equation of the lubricant film and the heat transfer equation of the bearing pad are handled by characteristic-based split method and conventional finite element method in three-dimensional space, respectively. Hot oil carry-over effect and variable lubricant viscosity are considered in the simulations. The results of the temperature distributions in the lubricant film and the bearing pad are presented. The possible usability of characteristic-based split method for future thermo-elasto-hydrodynamic analysis is discussed.

A Robust, Rapidly Convergent Method That Solves the Water Distribution Equations for Pressure-Dependent Models

AbstractIn the past, pressure-dependent models (PDMs) have suffered from convergence difficulties. In this paper conditions are established for the existence and uniqueness of solutions to the PDM problem posed as two optimization problems, one based on weighted least squares (WLS) and the other based on the co-content function. A damping scheme based on Goldstein’s algorithm is used and has been found to be both reliable and robust. A critical contribution of this paper is that the Goldstein theorem conditions guarantee convergence of the new method. The new methods have been applied to a set of eight challenging case study networks, the largest of which has nearly 20,000 pipes and 18,000 nodes, and are shown to have convergence behavior that mirrors that of the global gradient algorithm on demand-dependent model problems. A line search scheme based on the WLS optimization problem is proposed as the preferred option because of its smaller computational cost. Additionally, various consumption functions, i...

High-performance simulation of fracture in idealized ‘brick and mortar’ composites using adaptive Monte Carlo minimization on the GPU

Simulation of the nonlinear mechanical response of materials with explicit representation of microstructural features is extremely challenging. These models typically involve a very large number of degrees of freedom, and are prone to convergence difficulties when searching for roots to nonlinear equilibrium equations. We focus on an idealized material model that is motivated by the microstructure of synthetic nacre: individual ‘bricks’ (representing ceramic platelets) interact through nonlinear cohesive springs (representing a small volume fraction of polymer that bonds the platelets). The model simulates composite fracture through rupture of the cohesive springs. The problem is cast in terms of energy minimization and is essentially described by ‘nearest neighbor’ interactions. The principal focus of this paper is to illustrate the computational gains achievable by the strategic marriage of robust, global Monte Carlo minimization algorithms to the graphics processing unit architecture, and to describe how they were realized on the Nvidia GPU. Results comparing the computation times for graphics processing unit and central processing unit implementations demonstrate that a new adaptive version of the simulated annealing algorithm yields a speedup of approximately 5 times, whereas the graphics processing unit implementation yields a speed-up of about 16 times over conventional four-core central processing unit implementations. The resulting speed enhancement for adaptive graphics processing unit minimization of a factor of 80 enables a far broader range of simulations than has previously been possible. Simulations involving as many as 300,000 bricks can be performed in hours, as compared to weeks required by central processing unit implementation. Many aspects of this approach are translatable to other physical problems involving energy minimization in systems with large numbers of degrees of freedom.

Assessing Approaches to Learning in School Readiness

This study examines the psychometric properties of two assessments of children’s approaches to learning: the Devereux Early Childhood Assessment (DECA) and a 13-item approaches to learning rating scale (AtL) derived from the Arizona Early Learning Standards (AELS). First, we administered questionnaires to 1,145 randomly selected parents/guardians of first-time kindergarteners. Second, we employed confirmatory factor analysis (CFA) with parceling for DECA to reduce errors due to item specificity and prevent convergence difficulties when simultaneously estimating DECA and AtL models. Results indicated an overlap of 55% to 72% variance between the domains of the two instruments and suggested that the new AtL instrument is an easily administered alternative to the DECA for measuring children’s approaches to learning. This is one of the first studies that investigated DECA’s approaches to learning dimension and explored the measurement properties of an instrument purposely derived from a state’s early learning guidelines.

Multi-variate Hardy-type lattice point summation and Shannon-type sampling

The famous Shannon sampling theorem gives an answer to the question of how a one-dimensional time-dependent bandlimited signal can be reconstructed from discrete values in lattice points. In this work, we are concerned with multi-variate Hardy-type lattice point identities from which space-dependent Shannon-type sampling theorems can be obtained by straightforward integration over certain regular regions. An answer is given to the problem of how a signal bandlimited to a regular region in q-dimensional Euclidean space allows a reconstruction from discrete values in the lattice points of a (general) q-dimensional lattice. Weighted Hardy-type lattice point formulas are derived to allow explicit characterizations of over- and undersampling, thereby specifying not only the occurrence, but also the type of aliasing in a thorough mathematical description. An essential tool for the proof of Hardy-type identities in lattice point theory is the extension of the Euler summation formula to second order Helmholtz-type operators involving associated Green functions with respect to the “boundary condition” of periodicity. In order to circumvent convergence difficulties and/or slow convergence in multi-variate Hardy-type lattice point summation, some summability methods are necessary, namely lattice ball and Gaus–Weierstras averaging. As a consequence, multi-variate Shannon-type lattice sampling becomes available in a proposed summability context to accelerate the summation of the associated cardinal-type series. Finally, some aspects of constructive approximation in a resulting Paley–Wiener framework are indicated, such as the recovery of a finite set of lost samples, the reproducing Hilbert space context of spline interpolation.

A Priori Identifiability of Target-Mediated Drug Disposition Models and Approximations.

A priori identifiability of mathematical models assures that for a given input/output experiment, the parameter set has one unique solution within a defined space, independent of the experimental design. Many biologic therapeutics exhibit target-mediated drug disposition (TMDD), and use of the full compartmental model describing this system is well documented. In practice, estimation of the full parameter set for TMDD models, given real-world clinical data, is characterized by convergence difficulties and unstable solutions. Still, the formal assessment of the a priori identifiability of these systems has yet to be reported. The exact arithmetic rank (EAR) approach was used to test the a priori identifiability of a TMDD model as well as model approximations. The full TMDD and quasi-equilibrium/rapid binding (QE/RB), quasi-steady state (QSS), and Michaelis-Menten (MM) approximations were fully identifiable, a priori, regardless of whether observations were taken from a single or multiple compartments. The results of these identifiability analyses indicated that the difficulty with TMDD model convergence, a posteriori, lies in the experimental design, not in the mathematical identifiability in the lack of samples from several compartments. Experiments can be tailored to resolve these structurally non-identifiable parameters, notwithstanding practical implementation challenges. This work highlights the importance of identifiability analyses, specifically how they can influence experimental design and selection of the appropriate model structure to describe a dynamic biological system.

A finite strain quadrilateral based on least-squares assumed strains

When compared with advanced triangle formulations (e.g. Allman triangle and Arnold MINI), specially formulated low order quadrilateral elements still present performance advantages for bending-dominated and quasi-incompressible problems. However, simultaneous mesh distortion insensitivity and satisfaction of the Patch test is difficult. In addition, many enhanced-assumed (EAS) formulations show hourglass patterns in finite strains for large values of compression or tension; EAS elements often present convergence difficulties in Newton iteration, particularly in the presence of high bulk modulus or nearly-incompressible plasticity. Alternatively, we discuss the adequacy of a new assumed-strain 4-node quadrilateral for problems where high strain gradients are present. Specifically, we use relative strain projections to obtain three versions of a selectively-reduced integrated formulation complying a priori with the patch test. Assumed bending behavior is directly introduced in the higher-order strain term. Elements make use of least-square fitting and are generalization of classical B‾ and F‾ techniques. We avoid ANS (assumed natural strains) by defining the higher-order strain in contravariant/contravariant coordinates with a fixed frame. The kinematical part of the constitutive updating is based on quadratic incremental Green–Lagrange strains. Linear tests and both hyperelastic and elasto-plastic constitutive laws are used to test the element in realistic cases.

Nonlinear analysis of multiphase transport in porous media in the presence of viscous, buoyancy, and capillary forces

Nonlinear convergence problems in numerical reservoir simulation can lead to unacceptably large computational time and are often the main impediment to performing simulation studies of large-scale problems. We analyze the nonlinearity of the discrete transport (mass conservation) equation for immiscible, incompressible, two-phase flow in porous media in the presence of viscous, buoyancy, and capillary forces. Although simulation problems are multi-dimensional with large numbers of cells and variables, we find that the essence of the nonlinear behavior can be understood by studying the discretized (numerical) flux function for the interface between two cells. The numerical flux is expressed in terms of the saturations of the two cells. Discontinuities in the first-order derivative of the flux function (referred to as kinks) and inflection lines are identified as the cause of convergence difficulty. These critical features (kinks and inflections) change the curvature of the numerical flux function abruptly, and can lead to overshoots, oscillations, or divergence in Newton iterations.Based on our understanding of the nonlinearity, a nonlinear solver is developed, referred to as the Numerical Trust Region (NTR) solver. The solver is able to guide the Newton iterations safely and efficiently through the different saturation ‘trust-regions’ delineated by the kinks and inflections. Specifically, overshoots and oscillations that often lead to convergence failure are avoided. Numerical examples demonstrate that our NTR solver has superior convergence performance compared with existing methods. In particular, convergence is achieved for a wide range of timestep sizes and Courant–Friedrichs–Lewy (CFL) numbers spanning several orders of magnitude. In addition, a discretization scheme is proposed for handling heterogeneities in capillary-pressure–saturation relationship. The scheme has less degree of nonlinearity compared with the standard Single-point Phase-based Upstream weighting scheme, leading to an improved nonlinear convergence performance especially when used together with our NTR solver. Our proposed numerical solution strategy that is based on the numerical flux and handles capillarity extends the previous work by Jenny et al. (2009) [6] and Wang and Tchelepi (2013) [7] significantly.

Adequacy assessment of mathematical models in the dynamics of litter decomposition in a tropical forest Mosaic Atlantic, in southeastern Brazil.

The study of litter decomposition and nutrient cycling is essential to know native forests structure and functioning. Mathematical models can help to understand the local and temporal litter fall variations and their environmental variables relationships. The objective of this study was test the adequacy of mathematical models for leaf litter decomposition in the Atlantic Forest in southeastern Brazil. We study four native forest sites in Parque Estadual do Rio Doce, a Biosphere Reserve of the Atlantic, which were installed 200 bags of litter decomposing with 20 × 20 cm nylon screen of 2 mm, with 10 grams of litter. Monthly from 09/2007 to 04/2009, 10 litterbags were removed for determination of the mass loss. We compared 3 nonlinear models: 1 - Olson Exponential Model (1963), which considers the constant K, 2 - Model proposed by Fountain and Schowalter (2004), 3 - Model proposed by Coelho and Borges (2005), which considers the variable K through QMR, SQR, SQTC, DMA and Test F. The Fountain and Schowalter (2004) model was inappropriate for this study by overestimating decomposition rate. The decay curve analysis showed that the model with the variable K was more appropriate, although the values of QMR and DMA revealed no significant difference (p > 0.05) between the models. The analysis showed a better adjustment of DMA using K variable, reinforced by the values of the adjustment coefficient (R2). However, convergence problems were observed in this model for estimate study areas outliers, which did not occur with K constant model. This problem can be related to the non-linear fit of mass/time values to K variable generated. The model with K constant shown to be adequate to describe curve decomposition for separately areas and best adjustability without convergence problems. The results demonstrated the adequacy of Olson model to estimate tropical forest litter decomposition. Although use of reduced number of parameters equaling the steps of the decomposition process, no difficulties of convergence were observed in Olson model. So, this model can be used to describe decomposition curves in different types of environments, estimating K appropriately.

Hybrid upwind discretization of nonlinear two-phase flow with gravity

Multiphase flow in porous media is described by coupled nonlinear mass conservation laws. For immiscible Darcy flow of multiple fluid phases, whereby capillary effects are negligible, the transport equations in the presence of viscous and buoyancy forces are highly nonlinear and hyperbolic. Numerical simulation of multiphase flow processes in heterogeneous formations requires the development of discretization and solution schemes that are able to handle the complex nonlinear dynamics, especially of the saturation evolution, in a reliable and computationally efficient manner. In reservoir simulation practice, single-point upwinding of the flux across an interface between two control volumes (cells) is performed for each fluid phase, whereby the upstream direction is based on the gradient of the phase-potential (pressure plus gravity head). This upwinding scheme, which we refer to as Phase-Potential Upwinding (PPU), is combined with implicit (backward-Euler) time discretization to obtain a Fully Implicit Method (FIM). Even though FIM suffers from numerical dispersion effects, it is widely used in practice. This is because of its unconditional stability and because it yields conservative, monotone numerical solutions. However, FIM is not unconditionally convergent. The convergence difficulties are particularly pronounced when the different immiscible fluid phases switch between co-current and counter-current states as a function of time, or (Newton) iteration. Whether the multiphase flow across an interface (between two control-volumes) is co-current, or counter-current, depends on the local balance between the viscous and buoyancy forces, and how the balance evolves in time. The sensitivity of PPU to small changes in the (local) pressure distribution exacerbates the problem. The common strategy to deal with these difficulties is to cut the timestep and try again. Here, we propose a Hybrid-Upwinding (HU) scheme for the phase fluxes, then HU is combined with implicit time discretization to yield a fully implicit method. In the HU scheme, the phase flux is divided into two parts based on the driving force. The viscous-driven and buoyancy-driven phase fluxes are upwinded differently. Specifically, the viscous flux, which is always co-current, is upwinded based on the direction of the total-velocity. The buoyancy-driven flux across an interface is always counter-current and is upwinded such that the heavier fluid goes downward and the lighter fluid goes upward. We analyze the properties of the Implicit Hybrid Upwinding (IHU) scheme. It is shown that IHU is locally conservative and produces monotone, physically-consistent numerical solutions. The IHU solutions show numerical diffusion levels that are slightly higher than those for standard FIM (i.e., implicit PPU). The primary advantage of the IHU scheme is that the numerical overall-flux of a fluid phase remains continuous and differentiable as the flow regime changes between co-current and counter-current conditions. This is in contrast to the standard phase-potential upwinding scheme, in which the overall fractional-flow (flux) function is non-differentiable across the boundary between co-current and counter-current flows.

Numerical simulation of non-adiabatic capillary tubes. Special emphasis on the near-saturation zone

The aim of this article is to present a distributed numerical model that simulates the thermal and fluid-dynamic phenomena inside non-adiabatic capillary tubes. The resolution approach is based on a two-phase flow model where the fluid domain is discretized in a one-dimensional way, and the governing equations (continuity, momentum, and energy) are solved by means of a step-by-step algorithm. The model explained herein consists of an improved and extended version of previous works (Escanes et al., 1995; García-Valladares et al., 2002a,b; Ablanque et al., 2010) including two additional features. On the one hand, it allows the simulation of the two typical geometric arrangements found in capillary-tube/suction-line heat exchangers (i.e. concentric and lateral). On the other hand, it has an enhanced capability to address the convergence difficulties found in distributed models at the near-saturation zone. This document presents the major numerical adaptations done to the model, a comprehensive validation of the two geometric configurations, the model performance when tackling the aforementioned numerical difficulties and finally, some numerical studies.

Selection of Intermodal Conductivity Averaging Scheme for Unsaturated Flow in Homogeneous Media

The nonlinear solvers in numerical solution of water flow in variably saturated soils are prone to convergence difficulties. Many aspects can give rise to such difficulties, like very dry initial conditions, which causes a steep pressure gradient and great variation of hydraulic conductivity occurs across the wetting front during the infiltration of water. So, the averaging method applied to compute hydraulic conductivity between two adjacent nodes of the computational grid is one of the most important issues influencing the accuracy of the numerical solution of one-dimension al unsaturated flow equation i.e., Richards’ equation. A number of averaging schemes such as arithmetic, geometric, harmonic and arithmetic mean saturation have been proposed in the literature for homogeneous soil. The resulting numerical schemes are evaluated in terms of accuracy and computational time. It can be seen that the averaging scheme in the framework of arithmetic approaches favorably to other methods for a range of test cases.

Substructural parameters and dynamic loading identification with limited observations

Convergence difficulty and available complete measurement information have been considered as two primary challenges for the identification of large-scale engineering structures. In this paper, a time domain substructural identification approach by combining a weighted adaptive iteration (WAI) algorithm and an extended Kalman filter method with a weighted global iteration (EFK-WGI) algorithm was proposed for simultaneous identification of physical parameters of concerned substructures and unknown external excitations applied on it with limited response measurements. In the proposed approach, according to the location of the unknown dynamic loadings and the partially available structural response measurements, part of structural parameters of the concerned substructure and the unknown loadings were first identified with the WAI approach. The remaining physical parameters of the concerned substructure were then determined by EFK-WGI basing on the previously identified loadings and substructural parameters. The efficiency and accuracy of the proposed approach was demonstrated via a 20-story shear building structure and 23 degrees of freedom (DOFs) planar truss model with unknown external excitation and limited observations. Results show that the proposed approach is capable of satisfactorily identifying both the substructural parameters and unknown loading within limited iterations when both the excitation and dynamic response are partially unknown.

Analysis of spudcan–footprint interaction in a single soil with nonlinear FEM

The footprints that remain on the seabed after offshore jack-up platforms completed operations and moved out provide a significant risk for any future jack-up installation at that site. Detrimental horizontal and/or rotational loads will be induced on the base cone of the jack-up platform leg (spudcan) in the preloading process where only vertical loads are normally expected. However, there are no specific guidelines on design of spudcan re-installation very close to or partially overlapping existing footprints. This paper presents a rational design approach for assessing spudcan–footprint interaction and the failure process of foundation in a single layer based on nonlinear finite element method. The relationship between the distance between the spudcan and the footprint and the horizontal sliding force has been obtained. Comparisons of simulation and experimental results show that the model in this paper can deal well with the combined problems of sliding friction contact, fluid–solid coupling, and convergence difficulty. The analytical results may be useful to jack-up installation workovers close to existing footprints.

A surface hopping algorithm for nonadiabatic minimum energy path calculations.

The article introduces a robust algorithm for the computation of minimum energy paths transiting along regions of near-to or degeneracy of adiabatic states. The method facilitates studies of excited state reactivity involving weakly avoided crossings and conical intersections. Based on the analysis of the change in the multiconfigurational wave function the algorithm takes the decision whether the optimization should continue following the same electronic state or switch to a different state. This algorithm helps to overcome convergence difficulties near degeneracies. The implementation in the MOLCAS quantum chemistry package is discussed. To demonstrate the utility of the proposed procedure four examples of application are provided: thymine, asulam, 1,2-dioxetane, and a three-double-bond model of the 11-cis-retinal protonated Schiff base.

A Straightforward Convergence Method for ICCG Simulation of Multiloop and Time-Stepping FE Model of Synchronous Generators with Simultaneous AC and Rectified DC Connections

Now electric machines integrate with power electronics to form inseparable systems in lots of applications for high performance. For such systems, two kinds of nonlinearities, the magnetic nonlinearity of iron core and the circuit nonlinearity caused by power electronics devices, coexist at the same time, which makes simulation time-consuming. In this paper, the multiloop model combined with FE model of AC-DC synchronous generators, as one example of electric machine with power electronics system, is set up. FE method is applied for magnetic nonlinearity and variable-step variable-topology simulation method is applied for circuit nonlinearity. In order to improve the simulation speed, the incomplete Cholesky conjugate gradient (ICCG) method is used to solve the state equation. However, when power electronics device switches off, the convergence difficulty occurs. So a straightforward approach to achieve convergence of simulation is proposed. At last, the simulation results are compared with the experiments.