Numerical Difficulties Research Articles

Modeling of failure at the interface of ductile materials by applying the cohesive discontinuous Galerkin method

In this study, the failure behavior at the interface of ductile materials is investigated. In order to capture the degradation of the tractions at the interface, a cohesive zone (CZ) model is applied. The choice of the type of the CZ approach, i.e. either intrinsic or extrinsic, brings about different drawbacks. The former includes an elastic regime at the interface prior to the failure, which can result in numerical difficulties whereas the latter necessitates the re-meshing of the structure during crack propagation. In order to overcome these problems, the incomplete interior penalty Galerkin variant of the discontinuous Galerkin (DG) method is applied both at the interface and in the bulk instead of the standard conforming finite element method. In addition, the application of the DG method enables to use nonmatching meshes in the discretized model. To treat the bulk, an elastoplastic material model with isotropic hardening as well as different hardening rules for small strains is incorporated into the DG framework. Two numerical examples are computed to study the convergence behavior of the new cohesive discontinuous Galerkin (CDG) method in comparison to that of the conventional models. The new CDG method outperforms the conventional CZ continuous Galerkin elements in the presence of locking effects as well as hanging nodes.

Application of computational modeling for large wood dynamics with collisions on moveable channel beds

This study aims to analyze model reproducibility on large wood collisions on moveable channel beds. Flume experiments were conducted to observe the responses of moveable bed change with the behavior of large wood, including wood collisions. In the experiments, we employed time intervals of 5 and 9 s in wood supply, which modified the collision patterns. Subsequently, we conducted simulations based on the experimental conditions. Herein, to analyze the wood collision effect in the simulations, we employed the dashpot-spring model in a large wood dynamics model. Subsequently, we compared the simulations with and without wood collision and discussed the reproducibility of the model. The results of the simulations and experiments revealed that the wood pieces exhibited diverse behaviors, including floating, sliding, deposition, plane rotation, rolling, collision, remobilization, jam formation, and causing local bed changes. In particular, the results indicate that model reproducibility is dependent on wood collision because this can actively alter wood deposition patterns. We also described the piping and tunnel erosion processes linked to wood collision and moveable beds and suggested the corresponding, suitable computational treatment. This study presents an experimental and modeling methodology on wood collision and the interaction between bed change and large wood. Thus, useful information for solving numerical difficulties linked to large wood dynamics, including wood collisions, is provided.

Exponential stability and numerical analysis of a thermoelastic diffusion beam with rotational inertia and second sound

We study the dynamic behavior of a thermoelastic diffusion beam with rotational inertia and second sound, clamped at one end and free to move between two stops at the other. The contact with the stops is modeled with the normal compliance condition. The system, recently derived by Aouadi (2015), describes the behavior of thermoelastic diffusion thin plates under Cattaneo’s law for heat and mass diffusion transmission to remove the physical paradox of infinite propagation speeds of the classical Fourier’s and Fick’s laws. The system of equations is a coupling of a hyperbolic equation with four parabolic equations. It poses some new mathematical and numerical difficulties due to the lack of regularity and the nonlinear boundary conditions. The exponential stability of the solutions to the contact problem is obtained in the presence of rotational inertia thanks to a structural damping term. We propose a finite element approximation and we prove that the associated discrete energy decays to zero. Finally, we give an error estimate assuming extra regularity on the solution and we present some results of our numerical experiments.

Stabilised Variational Multi-scale Finite Element Formulations for Viscoelastic Fluids

The objective of this article is to summarise the work that we have been doing as a group in the context of stabilised finite element formulations for viscoelastic fluid flows. Viscoelastic fluids are complex non-Newtonian fluids, characterised by having an irreducible constitutive equation that needs to be solved coupled with the momentum and continuity equations. The finite element approximation of this kind of fluids presents several numerical difficulties. It inherits obviously the problems associated with the approximation of the incompressible Navier–Stokes equations. But, on top of that, now the constitutive equation is highly non-linear, with an advective term that may lead to both global and local oscillations in the numerical approximation. Moreover, even in the case of smooth solutions, it is necessary to meet some additional compatibility conditions between the velocity and the stress interpolation in order to ensure control over velocity gradients. The stabilised methods detailed in this work allow one to use equal order or even arbitrary interpolation for the problem unknowns ( $${\varvec{\sigma }}$$ - $${\varvec{u}}$$ -p) (elastic deviatoric stress-velocity-pressure) and to stabilise dominant convective terms, and all of them can be framed in the context of variational multi-scale methods. Some additional numerical ingredients that are introduced in this article are the treatment of the non-linearities associated with the problem and the possibility to introduce a discontinuity-capturing technique to prevent local oscillations. Concerning the constitutive equation, both the standard as the logarithmic conformation reformulation are discussed for stationary and time-dependent problems, and different versions of stabilised finite element formulations are presented in both cases.

Novel Nonlinear Complementarity Function Approach for Mechanical Analysis of Tensegrity Structures

Cable slacking may be frequently encountered in mechanical analysis of tensegrity structures. Traditional nonlinear finite element analysis in handling this issue may lead to numerical difficulty due to the dramatic changes of the stiffness matrix induced by cables switching between taut and slack states. To deal with the nonsmooth nature induced by cable slacking, we propose in this paper a novel nonlinear complementarity function-based framework for mechanical analysis of tensegrity structures. The core idea is that the constitutive relation of cables is reformulated as a set of equality equations. This is done by first introducing a parametric variable and expressing the constitutive relation with a complementarity condition, and then using the modified Fischer–Burmeister complementarity function. By considering the introduced parametric variables as additional degrees of freedom in the finite element model, it becomes possible to get rid of the nonsmooth nature, as well as to solve a continuous, smooth, differentiable problem by using the classical Newton–Raphson method. The proposed framework is first developed for static analysis and then extended to dynamic analysis, both under large displacement. The closed-form expressions of the Jacobian matrix are analytically derived. Multiple examples are presented to verify the robustness and effectiveness of the proposed framework.

Lagrangian particle model for 3D simulation of pellets and SPI fragments in tokamaks

A 3D numerical model for the ablation of pellets and shattered pellet injection fragments in tokamaks in the plasma disruption mitigation and fueling parameter space has been developed based on the Lagrangian particle (LP) method Samulyak et al (2018 J. Comput. Phys. 362 1–19). The pellet code implements the low magnetic Reynolds number MHD equations, kinetic models for the electronic heating, a pellet surface ablation model, an equation of state that supports multiple ionization states, radiation, and a model for grad-B drift of the ablated material across the magnetic field. The LP algorithm is highly adaptive, capable of simulating a large number of fragments in 3D while eliminating numerical difficulties of dealing with the tokamak background plasma. The code has achieved good agreement with theory for spherically symmetric ablation flows. Axisymmetric simulations of neon and deuterium pellets in magnetic fields ranging from 1 to 6 Tesla have been compared with previous simulations using the FronTier code, and very good agreement has also been obtained. The main physics contribution of the paper is a detailed study of the influence of 3D effects, in particular grad-B drift, on pellet ablation rates and properties of ablation clouds. Smaller reductions of ablation rates in magnetic fields compared to axially symmetric simulations have been demonstrated because the ablated material is not confined to narrowing channels in the presence of grad-B drift. Contribution of various factors in the grad-B drift model has also been quantified.

A First Violation Contact Algorithm that Correctly Captures History Dependence

Problems of elastic frictional contact are of great engineering relevance playing major roles in many fields such as tribology, attachment integrity, and even structural dynamics. These problems are historically extremely difficult to analyze either theoretically or numerically. A quick review of the literature provides insight as to a major cause of the numerical difficulty specifically that traditional algorithms search for contact solutions that satisfy all equality and inequality constraints using iterative processes that pass through states that do not satisfy those constraints and are physically inadmissible. A novel contact algorithm - the Method of First Violation - is introduced to address the deficiency identified in traditional methods by choosing load steps precisely so as to avoid iterations through inadmissible states. The new algorithm is demonstrated on multiple contact problems and the results compared to those of a conventional algorithm and to analytic solutions where available. Finally, the Method of First Violation is demonstrated on a contact problem with a nontraditional friction model (Dahls model) to demonstrate generality and robustness.

Friction Force During Lubricated Steady Sliding of a Rigid Cylinder on a Viscoelastic Substrate

We study the friction force during lubricated sliding of a rigid cylindrical indenter against a viscoelastic substrate in the iso-viscous visco-elasto-hydrodynamic lubrication (VEHL) regime. The substrate is represented by a foundation model. The solution is controlled by three dimensionless parameters. The first of these, λ, measures the time for the indenter to move one contact zone relative to the viscoelastic relaxation time; the second is the ratio of the long time to short time compliance of the substrate, $$c_{\infty } /c_{0}$$ ; the third parameter, β, is the ratio of average fluid flow rate to the sliding velocity. Although our solution works well for the full range of parameters, we focus on the “Hertz” regime (β >>1) where practically all the fluid in the contact region is squeezed out. This regime is quite common in soft contact lubrication problems and presents significant numerical difficulties. Our analysis gives insight into why these numerical difficulties arise. The friction force can be decomposed into two parts, one due to viscoelastic dissipation and the other from hydrodynamics. Although these two are generally coupled, in the Hertz limit, an important result is that the viscoelastic portion of the friction force can be well approximated by the solution of the corresponding “dry” sliding problem, in which there is no lubricating fluid layer. This provides a simple way to decouple the hydrodynamic portion of the friction force from the viscoelasticity of the substrate. We study how hydrodynamic pressure and film thickness vary with the controlling dimensionless parameters. Scaling laws for these relationships are given in closed form.

A normalization strategy for BESO-based structural optimization and its application to frequency response suppression

The bi-directional evolutionary structural optimization (BESO) method has been widely studied and applied due to its efficient iteration and clear boundaries. However, due to the use of the discrete design variable, numerical difficulties are more likely to occur with this method, especially in cases with strong nonlinearity. This limits the application of the BESO method in certain cases, such as the suppression of structural dynamic frequency response under high-frequency excitation. In this work, a normalization strategy is proposed for the BESO-based topology optimization, by which the magnitude of the sensitivities can be efficiently unified to the same order to avoid the possible numerical instabilities caused by the nonlinearity. To validate its merit in applications, the normalization-based BESO (NBESO) method is proposed for minimizing the structural frequency response. By means of the weighted sum method, a normalized weighted sum method is also proposed for multi-frequency involved problems. A series of 2D and 3D numerical examples is presented to illustrate the advantages of the NBESO. The effectiveness of the NBESO for multi-frequency response suppression is also demonstrated, in which the frequency ranges below and above the eigenfrequency are involved, respectively.

A new coupled multiphase flow–finite strain deformation–fault slip framework for induced seismicity

Production of hydrocarbons and water from subsurface reservoirs are known to cause permanent deformation of the reservoir and seismicity along faults both of which are detrimental to sustainable development of natural resources. Most of the prior studies on understanding fluid flow-induced plasticity and seismicity have focused on one or the other phenomenon due to the numerical difficulty associated with simultaneous modeling of the two failure phenomena because stress and deformation evolve non-linearly in both plasticity and seismicity. However, in reservoirs undergoing long-term production, plastic failure can alter the stress paths of points on a fault such that the onset, location, and magnitude of actual seismic events can no longer be predicted by a poroelastic simulation due to inaccurate stress and deformation history. We present a computational framework for coupled multiphase flow, finite strain poroplastic deformation, and dynamic fault slip and use it to understand the impact of plastic deformation on the onset, location, and magnitude of induced fault slip events. We evaluate the impact of plasticity on reservoir pressure, deformation, induced stress, and fault slip by comparing infinitesimal strain elastic and finite strain poro-elastoplastic models. For real-world applications, we consider different scenarios where the reservoir is either mechanically weaker or stronger than the caprock. We analyze the stress evolution as a function of the change in reservoir pressure to understand the role of contrast in reservoir and caprock elastic moduli on geomechanical stability. The results show that the poroplastic reservoir exhibits larger vertical deformation and delayed slip than the poroelastic reservoir after the same amount of oil production. For the same amount of pressure drop, a reservoir with a smaller modulus than the caprock displays a larger vertical displacement and an earlier onset of both plastic failure and fault slip. For a reservoir with a larger modulus than the caprock, vertical displacement is larger on the reservoir top boundary and smaller on the ground surface, and a higher pressure drop is needed to induce plastic failure and fault slip.

Modelling hollow pultruded FRP profiles under axial compression: Local buckling and progressive failure

Pultruded Fibre-Reinforced Polymer (PFRP) profiles are widely used as structural elements in many civil infrastructure applications. However, the anisotropic elasticity and the application-driven slenderness make these profiles prone to local buckling failure, well below their ultimate load capacity. In this paper, a numerical study was undertaken to characterise the local buckling and compressive failure of hollow PFRP profiles under axial compression. The Newton method were used along with the adaptive automatic stabilisation scheme and a controlled increment size in Abaqus 2019, to overcome the numerical difficulties in simulating local buckling. The numerical predictions were validated by experiments. The energy parameters and the constituent failure modes of the FEM models were used to explain the effect of dimension, layup, and slenderness ratio on the post-peak behaviour and failure modes of the PFRP profiles. Moreover, the reliance of the strain energy restoration after buckling and the axial and transverse deformations on these parameters was explained.

Structural topology optimization considering both performance and manufacturability: strength, stiffness, and connectivity

Structural topology optimization considering both performance and manufacturability is very attractive in engineering applications. This work proposes a formulation for structural topology optimization to achieve such a design, in which material strength, structural stiffness, and connectivity are simultaneously considered by integrating stress and simply-connected constraints into the compliance minimization problem. An effective solution algorithm consisting of different optimization techniques is introduced to handle various numerical difficulties resulted from this relatively complex multi-constraint and multi-field problem. Except for the stress penalization and aggregation techniques, the regional measure strategy is used together with the stability transformation method-based correction scheme to address stress constraints, which is also applied to the Poisson equation-based scalar field constraint in the simply-connected constraint. Numerical examples are presented to assess the features of the achieved design along with the performance of the employed algorithm. Comparisons with pure compliance, pure stress, and pure connectivity designs are provided to illustrate differences arising in the proposed design with respect to traditional approaches, also the necessity. Innovative manufacturing-oriented designs with consideration of the strength, stiffness, and connectivity are now available.

Implicit highly-coupled single-ion Hall-MHD formulation for hybrid particle-in-cell codes

The rudiments of a particle-based single-fluid two-temperature magnetohydrodynamic (MHD) algorithm have been outlined in Thoma et al. (2013). The extension of this algorithm to include the effect of Hall physics is described in this paper. An implicit leapfrog version of the algorithm, which allows timesteps large compared to the resistive decay time and other relevant timescales, has recently been added to a hybrid particle-in-cell code. In standard MHD the Hall term in the generalized Ohm’s law can often be neglected when the Hall parameter is small. This term must, however, be retained in regimes where it is non-negligible. The retention of displacement current in Maxwell’s equations avoids the numerical difficulties associated with the whistler mode, which are encountered in standard explicit Hall-MHD codes, and allows the algorithm to be incorporated into hybrid particle-in-cell codes, for which particles may migrate from a kinetic to fluid to MHD description based upon local ambient plasma conditions. A highly-coupled implicit Hall-MHD formalism is presented, in which displacement current can either be retained or neglected. Even when displacement current is neglected, the highly-coupled implicit formalism avoids the restrictive timesteps for the whistler mode in explicit Hall-MHD codes. A comparison of numerical and analytic dispersion analysis demonstrates the feasibility of this approach and establishes relevant constraints to assure numerical stability. The implementation of the algorithm is described, and test simulation results in 1D and 2D in both linear and nonlinear regimes are presented.

Cone-complementary manifold method for stability and failure analysis of jointed/fractured rock masses

Taking advantage of a unique dual-cover system, the numerical manifold method (NMM) is a powerful tool for analyzing stability and progressive failure of jointed/fractured rock masses. In this study, a new implicit cone complementary (CC) formulation is integrated into the NMM to capture the dynamic frictional-cohesive contact problems in jointed/fractured rock masses. The innovative numerical method, termed as the CC-NMM, circumvents numerical difficulties of using artificial contact penalty and open-close iteration in the original NMM. Benchmark examples demonstrate that both mechanical energy and translational momentum of the system are conserved in the proposed CC-NMM. Several examples are further designed to illustrate the excellent performance of the CC-NMM in simulating the progressive failure pattern and stability of the jointed/fractured rock masses.

Prestack waveform inversion based on analytical solution of the viscoelastic wave equation

Amplitude-variation-with-offset (AVO) inversion is based on single interface reflectivity equations. It involves some restrictions, such as the small-angle approximation, including only primary reflections, and ignoring attenuation. To address these shortcomings, the analytical solution of the 1D viscoelastic wave equation is used as the forward modeling engine for prestack inversion. This method can conveniently handle the attenuation and generate the full wavefield response of a layered medium. To avoid numerical difficulties in the analytical solution, the compound matrix method is applied to rapidly obtain the analytical solution by loop vectorization. Unlike full-waveform inversion, the proposed prestack waveform inversion (PWI) can be performed in a target-oriented way and can be applied in reservoir study. Assuming that a Q value is known, PWI is applied to synthetic data to estimate elastic parameters including compressional wave (P-wave) and shear wave (S-wave) velocities and density. After validating our method on synthetic data, this method is applied to a reservoir characterization case study. The results indicate that the reflectivity calculated by our approach is more realistic than that computed by using single interface reflectivity equations. Attenuation is an integral effect on seismic reflection; therefore, the sensitivity of seismic reflection to P-and S-wave velocities and density is significantly greater than that to Q, and the seismic records are sensitive to the low-frequency trend of Q. Thus, we can invert for the three elastic parameters by applying the fixed low-frequency trend of Q. In terms of resolution and accuracy of synthetic and real inversion results, our approach performs superiorly compared to AVO inversion.

Higher-order derivatives of the Green function in hyper-singular integral equations

Hyper-singular integral equations are often applied in the frequency-domain wave diffraction/radiation analyses of marine structures with thin plates or shell sub-structures. Their numerical solutions require the higher-order derivatives of the free-surface Green function featuring hyper-singularity, and hence the corresponding evaluation is very challenging. To circumvent the associated numerical difficulties, this paper will propose alternative formulations for the higher-order derivatives of both free-surface and Rankine-source parts of the Green function. For the free-surface term GF, the higher-order derivatives are analytically expressed by a combination of GF itself and its first-order horizontal radial derivative. Further, we derive an asymptotic representation, enabling us to deal exactly with a removable singularity in this representation. The superiority of the proposed formulation is demonstrated by comparing with a conventional direct differentiation. For the Rankine-source term, analytical expressions for the velocities induced by a uniform dipole distribution over a flat panel (involving second derivatives of the Rankine source term) are presented, which is directly relevant to numerical implementation based on constant panel methods. As illustrative examples, linear hydrodynamic coefficients of submerged circular impermeable and perforated plates are calculated for verification purposes. The proposed formulas are simple and easy to implement in the hyper-singular integral equations.

High precision implementation of Steck's recursion method for use in goodness-of-fit tests

Classical continuous goodness-of-fit (GOF) testing is employed for examining whether the data come from an assumed parametric model. In many cases, GOF tests assume a uniform null distribution and examine extreme values of the order statistics of the samples. Many of these statistics can be expressed by a function of the order statistics and the p-values amount to a joint probability statement based on the uniform order statistics. In this paper, we utilize Steck's recursion method and propose two high precision computing algorithms to compute the p-values for these GOF statistics. The numerical difficulties in implementing Steck's method are discussed and compared with solutions provided in high precision libraries.

A three-dimensional model of gas foil bearings and the effect of misalignment on the static performance of the first and second generation foil bearings

This paper presents a three-dimensional model of gas foil bearings with consideration of axial variation of the structural deflection. The top and bump foils are modeled with shell elements; hence, membrane, bending and shear effects are involved. The penalty method and node-to-segment contact scheme are used for contact modeling. To avoid numerical difficulties caused by discontinuity, the Coulomb law is modified by a regularization technique. Due to the nonlinear behavior of friction, the incremental iterative method is adopted for solution. The effect of misalignment on the static performance of the first and second generation foil bearings is studied. The results show that misalignment decreases the bearing load capacity and the second generation foil bearings exhibit better performance when misalignment occurs.

A new second-order modified hydrostatic reconstruction for the shallow water flows with a discontinuous topography

We consider a new second-order hydrostatic reconstruction (HR) based on the bottom topography, the depth-averaged velocity, and the water surface level for the shallow water flows with discontinuous bottom topography. The Audusse's second-order hydrostatic reconstruction (HR) (Audusse et al. (2004) [1]) fails to correctly reflect the wave pattern for the large discontinuous bottom topography. We propose a new modified HR method to overcome this numerical difficulty. The well-balancing property is obtained by using the information of the neighboring cells which is used to modify the reconstruction at the wet-dry front. The positivity preserving property is ensured by a “draining” time-step technique. A number of classical problems are successfully solved by the new scheme. Especially, the new scheme yields superior results for the shallow downhill flow over steps.

Differentiable scattering matrix for optimization of photonic structures

The scattering matrix, which quantifies the optical reflection and transmission of a photonic structure, is pivotal for understanding the performance of the structure. In many photonic design tasks, it is also desired to know how the structure's optical performance changes with respect to design parameters, that is, the scattering matrix's derivatives (or gradient). Here we address this need. We present a new algorithm for computing scattering matrix derivatives accurately and robustly. In particular, we focus on the computation in semi-analytical methods (such as rigorous coupled-wave analysis). To compute the scattering matrix of a structure, these methods must solve an eigen-decomposition problem. However, when it comes to computing scattering matrix derivatives, differentiating the eigen-decomposition poses significant numerical difficulties. We show that the differentiation of the eigen-decomposition problem can be completely sidestepped, and thereby propose a robust algorithm. To demonstrate its efficacy, we use our algorithm to optimize metasurface structures and reach various optical design goals.