Time Step Size Research Articles

Dynamic Phase‐Field Fracture in Viscoelastic Materials using a First‐Order Formulation

AbstractIn this contribution we present analytical results on a model for dynamic fracture in viscoelastic materials at small strains that have been obtained in full depth in [1]. In the model, the sharp crack interface is regularized with a phase‐field approximation, and for the phase‐field variable a viscous evolution with a quadratic dissipation potential is employed. A non‐smooth penalization prevents material healing. The viscoelastic momentum balance is formulated as a first order system and coupled in a nonlinear way to the non‐smooth evolution equation of the phase field. We give a full discretization in time and space using a discontinuous Galerkin method for the first‐order system. We discuss the existence of discrete solutions and, with the step size in space and time tending to zero, their convergence to a suitable notion of weak solution of the system. Eventually, we provide a numerical benchmark and compare it with simulation results found in [2].

A Stability Analysis of Modified Patankar–Runge–Kutta methods for a nonlinear Production–Destruction System

AbstractModified Patankar–Runge–Kutta (MPRK) methods preserve the positivity as well as conservativity of a production–destruction system (PDS) of ordinary differential equations for all time step sizes. As a result, higher order MPRK schemes do not belong to the class of general linear methods, i. e. the iterates are generated by a nonlinear map g even when the PDS is linear. Moreover, due to the conservativity of the method, the map g possesses non‐hyperbolic fixed points.Recently, a new theorem for the investigation of stability properties of non‐hyperbolic fixed points of a nonlinear iteration map was developed. We apply this theorem to understand the stability properties of a family of second order MPRK methods when applied to a nonlinear PDS of ordinary differential equations. It is shown that the fixed points are stable for all time step sizes and members of the MPRK family. Finally, experiments are presented to numerically support the theoretical claims.

First steps towards modeling the interaction of cardiovascular agents and smooth muscle activation in arterial walls

AbstractNumerical simulation of the response of healthy and pathological arteries to cardiovascular agents can provide valuable information to the physician in the treatment of diseases such as hypertension, atherosclerosis, and the Marfan syndrome. Here, we provide a first step towards a computational framework to model the effects of antihypertensive agents on the mechanical response of arterial walls. A material model is developed by extending an existing formulation for wall tissue to incorporate the effects of calcium‐ion channel blockers. The resulting coupled deformation‐diffusion problem is then solved using the finite element method. Simulation results with drug activity show that, indeed, an increased lumen diameter due to reduced contraction is obtained. Additionally, a decrease in the rate of arterial contraction is observed, which is also consistent with expected behavior. Finally, we compare results for an implicit or explicit treatment of the the deformation‐diffusion coupling, and we observe that both coupling schemes yield comparable results for a wide range of time step sizes.

Deep Neural Network Modeling for CFD Simulations: Benchmarking the Fourier Neural Operator on the Lid-Driven Cavity Case

In this work we present the development, testing and comparison of three different physics-informed deep learning paradigms, namely the ConvLSTM, CNN-LSTM and a novel Fourier Neural Operator (FNO), for solving the partial differential equations of the RANS turbulence model. The 2D lid-driven cavity flow was chosen as our system of interest, and a dataset was generated using OpenFOAM. For this task, the models underwent hyperparameter optimization, prior to testing the effects of embedding physical information on performance. We used the mass conservation of the model solution, embedded as a term in our loss penalty, as our physical information. This approach has been shown to give physical coherence to the model results. Based on the performance, the ConvLSTM and FNO models were assessed in forecasting the flow for various combinations of input and output timestep sizes. The FNO model trained to forecast one timestep from one input timestep performed the best, with an RMSE for the overall x and y velocity components of 0.0060743 m·s−1.

Spatial and temporal constraints of the cohesive modeling: A unified criterion for fluid‐driven fracture

AbstractWe present a unified criterion for cohesive modeling of fluid‐driven fracture based on the dimensional analysis to simultaneously provide the constraint for cohesive element and time step sizes. Complicated by the nonlinear interaction between solid deformation and fluid flow, the underlying correlation between discretization and physical parameters of fluid‐driven fracture is still unclear. This work studies this correlation through the dimensionless process of the governing equations that associate the cohesive element and time step sizes in a discrete regime. Three characteristic parameters (i.e., related to crack opening, fluid pressure, and fracture length) are introduced in the derivation, and two dimensionless parameters are proposed to construct the unified criterion. The criterion is validated by numerical tests of toughness‐dominated fracture with various conditions including the modulus of solid, injection rate of fluid, fracture energy, and in‐situ stress. The proposed criterion determines the spatial and temporal constraints of the cohesive zone model for modeling fluid‐driven fracture, which is often treated empirically in previous practices.

An energy stable, conservative and bounds‐preserving numerical method for thermodynamically consistent modeling of incompressible two‐phase flow in porous media with rock compressibility

SummaryIn this paper, we consider modeling and numerical simulation of incompressible and immiscible two‐phase flow in porous media with rock compressibility. Using the second law of thermodynamics, we rigorously derive a thermodynamically consistent mathematical model, which characterizes the two‐phase capillarity and rock compressibility through free energies, thereby following an energy dissipation law. We also derive a general and thermodynamically consistent formulation for the effective pore fluid pressure acting on rocks, which is a fundamental problem for two‐phase flow with rock compressibility. To solve the model effectively, we propose an energy stable numerical method, which can preserve multiple physical properties, including the energy dissipation law, full conservation law for both fluids and pore volumes, and positivity of porosity and saturations. Benefiting from the newly‐developed energy factorization approach and careful treatments for the effective pressure and porosity, the proposed scheme can inherit the energy dissipation law at the discrete level. The fully discrete scheme is constructed using a locally conservative cell‐centered finite difference method. The implicit strategy is applied to treat the upwind phase mobilities and porosity in the phase mass conservation equations and the porosity equation so as to conserve the mass of each phase as well as pore volumes. The positivity of porosity and saturations is proved without any restrictions on time step and mesh sizes. An efficient sequential iterative method is also developed to solve the nonlinear system resulting from the scheme. Finally, numerical results are given to verify the features of the proposed method.

Modeling earthquake cycles and wear on rough faults with the mortar finite element method

SUMMARYThis paper presents a mortar-based finite element formulation for modelling earthquake cycles and wear on rough faults governed by rate and state friction. The method allows large sliding on the fault and accounts for all stages in the earthquake cycle, using a variable time step size with a transition between quasi-static and fully dynamic time discretizations. Wear laws with linear and power law forms are discretized and implemented into the mortar method, as well as a minimum level of normal traction constraint to treat fault opening. We examine the effect of wear laws on the slip behaviour and near-fault stresses during simulations of earthquake cycles on rough faults. The simulations demonstrate that the implementation of wear allows more realistic modelling of the earthquake cycle without the development of unrealistically large stresses and with less reduction of earthquake magnitude with total fault slip. Moreover, the method enables to study the effects of roughness and fault slip on the gouge zone thickness.

Jacobian-free Newton–Krylov method for the simulation of non-thermal plasma discharges with high-order time integration and physics-based preconditioning

A preconditioning framework for the numerical simulation of non-thermal streamer discharges is developed using the Jacobian-free Newton-Krylov (JFNK) method. A reduced plasma fluid model is considered, consisting of electrons, one positive ion, one negative ion, and the electrostatic potential. The plasma kinetics model includes ionization, electron-ion recombination, electron attachment, electron detachment, and ion-ion recombination. The governing equations are made dimensionless, discretized in space with finite differences, and integrated in time with a fully implicit method based on high-order backward differentiation formulas. The preconditioning framework is based on a linearized form of the governing equations and physics-based operator splitting. The efficiency of the preconditioning strategy is assessed through two test cases: streamer propagation between parallel plates and an axisymmetric pin-to-pin discharge. The fully implicit approach overcomes traditional restrictions in the time step size due to processes such as electron drift, electron diffusion, and dielectric relaxation. Excellent performance is observed through relevant statistics of the JFNK solver, although the number of linear iterations increases for the pin-to-pin discharge when nonlinear numerical boundary conditions are imposed at the electrodes. Performance studies show scalability with O(100-1000) processors for O(10M) unknowns with ample room for optimization.

Nonlinear convolution finite element method for solution of large deformation elastodynamics

A new algorithm based on the convolution finite element method (CFEM) is proposed for the nonlinear wave propagation in elastic media. The formulation is developed in the context of the total Lagrangian framework, encompassing contributions due to both geometrical and material nonlinearities. As a basis, a counterpart of equations of motion – namely, the alternative field equations – is first established. The satisfaction of the alternative field equations is then realized in a weak sense. Next, the Newton–Raphson procedure and the consistent tangential matrix are applied to the weak formulation, where the CFEM is used as the linear solver in each iteration. Finally, several examples are carried out to examine the theoretical aspects and the feasibility of the proposed algorithm. In particular, problems of free vibration of Neo-Hookean and Saint Venant–Kirchhoff plates are explored. Also, a cantilever beam of the Neo-Hookean material is simulated for the case of forced vibrations. Conspicuously, in contrast to the existing time-step methods with finite order of accuracy, the new solution procedure obtains the accurate solution when the time-step size is increased.

Efficient dissipation-based nonlinear solver for multiphase flow in discrete fractured media

Solving nonlinear equation systems resulting from the Fully Implicit Method (FIM) remains a challenge for numerical simulation of multiphase flow in subsurface fractured media. The Courant numbers can vary orders of magnitude across discrete fracture-matrix (DFM) models because of the high contrasts in the permeability and length-scale between matrix and fracture. The standard Newton solver is usually unable to converge for big timestep sizes.Limited research has been conducted on nonlinear solver techniques for multiphase flow and transport in fractured media. We make an extension of a new dissipation-based continuation (DBC) method to DFM models. Our goal is to prevent timestep cuts and maintain efficient time-stepping for FIM. The DBC algorithm builds a homotopy of the discretized conservation equations through the addition of numerical dissipation terms. We introduce a continuation parameter for controlling the dissipation and ensuring that accuracy of the computed solution will not be reduced. Under the nonlinear framework of DBC, general dissipation operators and adaptive strategies are developed to provide the optimal dissipation matrix for the multiphase hyperbolic system.We assess the new nonlinear solver using multiple numerical examples. Results reveal that the damped-Newton solver suffers from serious restrictions on timestep sizes and wasted iterations. In contrast, the DBC solver provides excellent computational performance. The dissipation operators are able to successfully resolve main convergence difficulties. We also investigate the impact of star-delta transformation which removes the small cells at fracture intersections. Moreover, we demonstrate that an aggressive time-stepping does not affect the solution accuracy.

An efficient computational technique for solving a fractional-order model describing dynamics of neutron flux in a nuclear reactor

This paper focuses on the development of a computational technique for numerical solution of the time-fractional neutron diffusion (FND) equation, which describes dynamics of neutron flux in a nuclear reactor. The time-fractional derivatives are discretized by means of L1 scheme, while the space derivative is discretized by employing a collocation method based on quartic B-spline (QUBS) basis functions. Numerical experiments are presented to illustrate the performance of proposed technique. Results reveal that the proposed technique yields a reliable approximation to the solution of the underlying model problem and has O(k+h4)−order convergence, with k and h representing the step sizes in time and space directions, respectively. It is important to point out that the numerical method for the underlying problem has not yet been discussed in the literature. The effect of fractional order derivative on the behavior of neutron flux is investigated.

High-order numerical algorithm and error analysis for the two-dimensional nonlinear spatial fractional complex Ginzburg–Landau equation

In this paper, we first construct an appropriate new generating function, and then based on this function, we establish a fourth-order numerical differential formula approximating the Riesz derivative with order γ∈(1,2]. Subsequently, we apply the formula to numerically study the two-dimensional nonlinear spatial fractional complex Ginzburg–Landau equation and obtain a difference scheme with convergence order Oτ2+hx4+hy4, where τ denotes the time step size, hx and hy denote the space step sizes, respectively. Furthermore, with the help of some newly derived discrete fractional Sobolev embedding inequalities, the unique solvability, the unconditional stability, and the convergence of the constructed numerical algorithm under different norms are proved by using the discrete energy method. Finally, some numerical results are presented to confirm the correctness of the theoretical results and verify the effectiveness of the proposed scheme.

An Accurate, Controllably Dissipative, Unconditionally Stable Three-Sub-Step Method for Nonlinear Dynamic Analysis of Structures

An implicit truly self-starting time integration method for nonlinear structural dynamical systems is developed in this paper. The proposed method possesses unconditional stability, second-order accuracy, and controllable dissipation, and it has no overshoots. The well-known BN-stability theory is employed in the design of algorithmic parameters, ensuring that the proposed method can stably solve nonlinear structural dynamical systems without restricting the time step size. The spectral analysis shows that compared to existing second-order accurate time integration methods, the proposed method enjoys a considerable advantage in low-frequency accuracy. For nonlinear problems where the currently popular Generalized-[Formula: see text] method and [Formula: see text]-Bathe method fail, the proposed method shows strong stability and accuracy. Further, for nonlinear problems in which all methods’ results are convergent, the proposed method has greater accuracy, efficiency, and energy-conservation capability.

Adaptive energy-preserving algorithms for guiding center system

We develop two types of adaptive energy preserving algorithms based on the averaged vector field for the guiding center dynamics, which plays a key role in magnetized plasmas. The adaptive scheme is applied to the Gauss Legendre’s quadrature rules and time stepsize respectively to overcome the energy drift problem in traditional energy-preserving algorithms. These new adaptive algorithms are second order, and their algebraic order is carefully studied. Numerical results show that the global energy errors are bounded to the machine precision over long time using these adaptive algorithms without massive extra computation cost.

Bound-Preserving Discontinuous Galerkin Methods with Modified Patankar Time Integrations for Chemical Reacting Flows

In this paper, we develop bound-preserving discontinuous Galerkin (DG) methods for chemical reactive flows. There are several difficulties in constructing suitable numerical schemes. First of all, the density and internal energy are positive, and the mass fraction of each species is between 0 and 1. Second, due to the rapid reaction rate, the system may contain stiff sources, and the strong-stability-preserving explicit Runge-Kutta method may result in limited time-step sizes. To obtain physically relevant numerical approximations, we apply the bound-preserving technique to the DG methods. Though traditional positivity-preserving techniques can successfully yield positive density, internal energy, and mass fractions, they may not enforce the upper bound 1 of the mass fractions. To solve this problem, we need to (i) make sure the numerical fluxes in the equations of the mass fractions are consistent with that in the equation of the density; (ii) choose conservative time integrations, such that the summation of the mass fractions is preserved. With the above two conditions, the positive mass fractions have summation 1, and then, they are all between 0 and 1. For time discretization, we apply the modified Runge-Kutta/multi-step Patankar methods, which are explicit for the flux while implicit for the source. Such methods can handle stiff sources with relatively large time steps, preserve the positivity of the target variables, and keep the summation of the mass fractions to be 1. Finally, it is not straightforward to combine the bound-preserving DG methods and the Patankar time integrations. The positivity-preserving technique for DG methods requires positive numerical approximations at the cell interfaces, while Patankar methods can keep the positivity of the pre-selected point values of the target variables. To match the degree of freedom, we use $$Q^k$$ polynomials on rectangular meshes for problems in two space dimensions. To evolve in time, we first read the polynomials at the Gaussian points. Then, suitable slope limiters can be applied to enforce the positivity of the solutions at those points, which can be preserved by the Patankar methods, leading to positive updated numerical cell averages. In addition, we use another slope limiter to get positive solutions used for the bound-preserving technique for the flux. Numerical examples are given to demonstrate the good performance of the proposed schemes.

New stable, explicit, second order hopscotch methods for diffusion-type problems

The aim of this paper is to systematically construct and test novel odd–even hopscotch-type numerical algorithms solving the diffusion or heat equation. Among the studied explicit two-stage methods some of them are unconditionally stable and have second order convergence rate in time step size, which is proved analytically as well. We apply the best methods to the nonlinear Fisher’s equation to demonstrate that they work also for nonlinear equations. Then, in order to examine the competitiveness of the new algorithms, we test them for the heat equation against widely used numerical solvers in cases where the media are strongly inhomogeneous and thus the coefficients strongly depend on space. The results suggest that the new methods are significantly more effective than the widely used explicit or implicit methods, especially for extremely large stiff systems.

Parallel block preconditioners for virtual element discretizations of the time-dependent Maxwell equations

The focus of this study is the construction and numerical validation of parallel block preconditioners for low order virtual element discretizations of the three-dimensional Maxwell equations. The virtual element method (VEM) is a recent technology for the numerical approximation of partial differential equations (PDEs), that generalizes finite elements to polytopal computational grids. So far, VEM has been extended to several problems described by PDEs, and recently also to the time-dependent Maxwell equations. When the time discretization of PDEs is performed implicitly, at each time-step a large-scale and ill-conditioned linear system must be solved, that, in case of Maxwell equations, is particularly challenging, because of the presence of both H(div) and H(curl) discretization spaces. We propose here a parallel preconditioner, that exploits the Schur complement block factorization of the linear system matrix and consists of a Jacobi preconditioner for the H(div) block and an auxiliary space preconditioner for the H(curl) block. Several parallel numerical tests have been performed to study the robustness of the solver with respect to mesh refinement, shape of polyhedral elements, time step size and the VEM stabilization parameter.

Investigating the stability of an explicit ADE-FDTD scheme for modeling graphene: Avoiding erroneous conclusions

In recent years, an explicit auxiliary differential equation finite difference time domain (ADE-FDTD) scheme has been frequently used in graphene simulations. In this respect, a differencing scheme in which the magnetic field and the current density are collocated in time was developed and it has been recently stated that this explicit ADE-FDTD implementation is “conditionally stable with the maximum time-step size bounded by the common Courant–Friedrichs–Lewy (CFL) limit”. In this communication, the stability of this explicit ADE-FDTD scheme is revisited and it is shown that the conventional CFL constraint is not retained, and a time-step less than the CFL limit must be used in all domain grids to insure the stability in the whole system. Although this time-step stringent can be of little effect for some graphene applications, and in order to avoid drawing erroneous conclusions, it is shown that when the given ADE-FDTD scheme is used for simulating graphene structures that include noble and other plasmonic materials the derived time-step stringent can be of significant effect. To retain the CFL stability constraint, alternative formulation based on the Runge–Kutta (RK) time-differencing scheme is also presented in this communication. Contrary to the other existing CFL-stable ADE-FDTD approaches, the presented RK-FDTD method is equally applicable to both magnetized and unmagnetized graphene structures without incurring additional computational cost. Furthermore, the RK-FDTD scheme is extended for modeling the complex-frequency-shifted perfectly matched layer (CFS-PML) mesh truncating technique. Finally, the stability and accuracy of the explicit ADE/RK-FDTD schemes are demonstrated through numerical examples.

An improved multistage preconditioner on GPUs for compositional reservoir simulation

The compositional model is often used to describe multicomponent multiphase porous media flows in the petroleum industry. The fully implicit method with strong stability and weak constraints on time-step sizes is commonly used in mainstream commercial reservoir simulators. In this paper, we develop an efficient multistage preconditioner for the fully implicit compositional flow simulation. The method employs an adaptive setup phase to improve the parallel efficiency on GPUs. Furthermore, a multicolor Gauss–Seidel algorithm based on the adjacency matrix is applied in the algebraic multigrid methods for the pressure part. Numerical results demonstrate that the proposed algorithm achieves good parallel speedup while yielding the same convergence behavior as the corresponding sequential version.

Computational uncertainty and optimal grid size and time step of the Lax–Friedrichs scheme for the 1D advection equation

This paper examines truncation and round-off errors in the numerical solution of the 1D advection equation with the Lax–Friedrichs scheme, and accumulation of the errors as they are propagated to high temporal layers. The authors obtain a new theoretical approximation formula for the upper bound of the total error of the numerical solution, as well as theoretical formulae for the optimal grid size and time step. The reliability of the obtained formulae is demonstrated with numerical experimental examples. Next, the ratio of the optimal time steps under two different machine precisions is found to satisfy a universal relation that depends only on the machine precision involved. Finally, theoretical verification suggests that this problem satisfies the computational uncertainty principle when the grid ratio is fixed, demonstrating the inevitable existence of an optimal time step size under a finite machine precision.摘要本文对于应用Lax- Friedrichs 格式数值求解一维平流方程, 研究数值求解过程中产生的截断误差与舍入误差, 以及两种误差逐层向高时间层传播的累积, 得到新的数值解总误差上界的理论近似公式, 以及最优格距和最优时间步长的理论公式. 通过数值算例验证了所得公式的可靠性. 然后, 发现了两种不同机器精度下最优时间步长之比满足的一个仅与机器精度有关的普适关系. 最后, 理论验证了在网格比固定的情况下, 此问题满足数值计算的不确定性原理, 以及在机器有限精度下最优时间步长的必然存在.