Finite Difference Discretization Research Articles

A novel immersed boundary approach for irregular topography with acoustic wave equations

Irregular terrain has a pronounced effect on the propagation of seismic and acoustic wavefields but is not straightforwardly reconciled with structured finite-difference (FD) methods used to model such phenomena. Accurate wavefield simulation is paramount in subsurface imaging applications such as reverse time migration and full-waveform inversion, requiring suitable topography handling. Methods currently detailed in the literature generally are limited in scope application-wise or nontrivial to apply to real-world geometries. A general immersed boundary treatment capable of imposing a range of boundary conditions in a relatively equation-agnostic manner has been developed, alongside a framework implementing this approach to complement emerging code-generation paradigms. The approach is distinguished by the use of N-dimensional Taylor-series extrapolants constrained by boundary conditions imposed at some suitably distributed set of surface points. The extrapolation process is encapsulated in modified derivative stencils applied in the vicinity of the boundary, using hyperspherical support regions. This method ensures boundary representation is consistent with the FD discretization. Furthermore, high-dimensional and vector boundary conditions can be applied without approximation prior to discretization. A consistent methodology can thus be applied across free and rigid surfaces with first- and second-order acoustic wave equation formulations. Application to both equations is demonstrated, and numerical examples based on analytic and real-world topography implementing free and rigid surfaces in 2D and 3D are presented. Numerical examples and convergence tests demonstrate the accuracy of boundary treatments devised by the prescribed approach, their suitability to practical applications, and the feasibility of automatically generating treatments to suit each case.

A thermal energy analysis of binary (Go-Co/H2O) and ternary (Go-Co-Zro2/H2O) nanofluids based on characterization and thermal performance

Hybrid or binary nanofluids have superior mechanical and thermal characteristics but the tri-hybrid nanofluids comprise of more embellished thermal properties, better physical strength, and enhanced stability. The present work characterizes the thermal and physical aspects of the hybrid and tri-hybrid nanofluids. The nano-composition of graphene oxide ( Go) and cobalt ( Co) is used in the amalgamation of hybrid nanofluid Go-Co/H2O, whereas the addition of zirconium oxide ( ZrO2) in this mixture gives rise to the ternary Go-Co-ZrO2/H2O hybrid nanofluid. The activation energy and viscous dissipation terms are also amended in the governing equations. The mathematical framework consists of a complex natured dynamical system. However, a numerical algorithm based on finite-difference discretization is developed which can solve the system numerically via the MATLAB software. A comparison with the existing literature is provided in order to validate the numerical procedure. From the outcomes, it is noticed that the temperature of hybrid as well as tri-hybrid nanofuid increases rapidly with change in concentration of zirconium oxide and cobalt. Temperature increases up to 20% by taking 0.1 volume fraction of both zirconium oxide and cobalt. Porous medium and activation energy resist the flow and concentration respectively. A comparative judgment evidently reveals that tri-hybrid Go-Co-ZrO2/H2O nanofluid has a substantial effect on temperature as equated to hybrid or pure nanofluid.

Convergence analysis of a second order numerical scheme for the Flory–Huggins–Cahn–Hilliard–Navier–Stokes system

We present an optimal rate convergence analysis for a second order accurate in time, fully discrete finite difference scheme for the Cahn–Hilliard–Navier–Stokes (CHNS) system, combined with logarithmic Flory–Huggins energy potential. The numerical scheme has been recently proposed, and the positivity-preserving property of the logarithmic arguments, as well as the total energy stability, have been theoretically justified. In this paper, we rigorously prove second order convergence of the proposed numerical scheme, in both time and space. Since the CHNS is a coupled system, the standard ℓ∞(0,T;ℓ2)∩ℓ2(0,T;Hh2) error estimate could not be easily derived, due to the lack of regularity to control the numerical error associated with the coupled terms. Instead, the ℓ∞(0,T;Hh1)∩ℓ2(0,T;Hh3) error analysis for the phase variable and the ℓ∞(0,T;ℓ2) analysis for the velocity vector, which shares the same regularity as the energy estimate, is more suitable to pass through the nonlinear analysis for the error terms associated with the coupled physical process. Furthermore, the highly nonlinear and singular nature of the logarithmic error terms makes the convergence analysis even more challenging, since a uniform distance between the numerical solution and the singular limit values of is needed for the associated error estimate. Many highly non-standard estimates, such as a higher order asymptotic expansion of the numerical solution (up to the third order accuracy in time and fourth order in space), combined with a rough error estimate (to establish the maximum norm bound for the phase variable), as well as a refined error estimate, have to be carried out to conclude the desired convergence result. To our knowledge, it will be the first work to establish an optimal rate convergence estimate for the Cahn–Hilliard–Navier–Stokes system with a singular energy potential.

Analysis of a positivity-preserving splitting scheme for some semilinear stochastic heat equations

We construct a positivity-preserving Lie--Trotter splitting scheme with finite difference discretization in space for approximating the solutions to a class of semilinear stochastic heat equations with multiplicative space-time white noise. We prove that this explicit numerical scheme converges in the mean-square sense, with rate $1/4$ in time and rate $1/2$ in space, under appropriate CFL conditions. Numerical experiments illustrate the superiority of the proposed numerical scheme compared with standard numerical methods which do not preserve positivity. p, li { white-space: pre-wrap; }

Jump–diffusion productivity models in equilibrium problems with heterogeneous agents

In this paper we adopt a rational expectations framework to formulate general equilibrium models with heterogeneous agents. The productivity dynamics are characterized by a jump–diffusion model, thus allowing to account for sudden and impactful events. The modelling approach utilizes Hamilton–Jacobi–Bellman (HJB) formulations to represent the endogenous decision-making of firms to remain or exit the industry. When firms decide to exit, they are instantaneously replaced by new entrants. This dynamic leads to the development of a probability density function for firms, which satisfies a Kolmogorov–Fokker–Planck (KFP) equation with a source term. Both HJB and KFP formulations involve partial-integro differential operators due to the presence of jumps. Equilibrium models are completed with the household problem and feasibility conditions. Since (semi-)analytical solutions are not available, a numerical methodology is considered. This approach involves a Crank–Nicolson scheme for the time discretization, an augmented Lagrangian active set method and a finite difference discretization for the HJB formulation, and an appropriate finite difference method for the KFP problem. Moreover, Adams–Bashforth schemes are employed to handle integral terms explicitly. For the global equilibrium problem, we introduce a Steffensen algorithm. Numerical examples are provided to showcase the performance of our proposed numerical methodologies and to illustrate the expected behaviour of computed economic variables.

Fully conservative difference schemes for the rotation-two-component Camassa–Holm system with smooth/nonsmooth initial data

This paper derives a semi-discrete conservative difference scheme for the rotation-two-component Camassa–Holm system based on its Hamiltonian invariants. Mass, momentum and energy are preserved for the semi-discrete scheme. Furthermore, a fully discrete finite difference scheme is proposed without destroying any one of the conservative laws. Combining a nonlinear iteration with a threshold strategy, the accuracy of the scheme is guaranteed. Meanwhile, this scheme captures the formation and propagation of solitary wave solutions in long time behavior under smooth/nonsmooth initial data. Remarkably, a new type of asymmetric wave breaking phenomenon is revealed in the case of the nonzero rotational parameter.

A Novel Cell-Based Adaptive Cartesian Grid Approach for Complex Flow Simulations

As the need for handling complex geometries in computational fluid dynamics (CFD) grows, efficient and accurate mesh generation techniques become paramount. This study presents an adaptive mesh refinement (AMR) technology based on cell-based Cartesian grids, employing a distance-weighted least squares interpolation for finite difference discretization and utilizing immersed boundary methods for wall boundaries. This facilitates effective management of both transient and steady flow problems. Validation through supersonic flow over a forward-facing step, subsonic flow around a high Reynolds number NHLP airfoil, and supersonic flow past a sphere demonstrated AMR’s efficacy in capturing essential flow characteristics while wisely refining and coarsening meshes, thus optimizing resource utilization without compromising accuracy. Importantly, AMR simplified the capture of complex flows, obviating manual mesh densification and significantly improving the efficiency and reliability of CFD simulations.

Optimal finite-differences discretization for the diffusion equation from the perspective of large-deviation theory

When applying the finite-differences method to numerically solve the one-dimensional diffusion equation, one must choose discretization steps Δx, Δt in space and time, respectively. By applying large-deviation theory on the discretized dynamics, we analyze the numerical errors due to the discretization, and find that the (relative) errors are especially large in regions of space where the concentration of particles is very small. We find that the choice Δt=Δx2/(6D) , where D is the diffusion coefficient, gives optimal accuracy compared to any other choice (including, in particular, the limit Δt→0 ), thus reproducing the known result that may be obtained using truncation error analysis. In addition, we give quantitative estimates for the dynamical lengthscale that describes the size of the spatial region in which the numerical solution is accurate, and study its dependence on the discretization parameters. We then turn to study the advection–diffusion equation, and obtain explicit expressions for the optimal Δt and other parameters of the finite-differences scheme, in terms of Δx, D and the advection velocity. We apply these results to study large deviations of the area swept by a diffusing particle in one dimension, trapped by an external potential ∼|x| . We extend our analysis to higher dimensions by combining our results from the one dimensional case with the locally one-dimension method.

A comparative study on the acceleration techniques for solving finite difference discretization poisson’s equation in the PIC/MCC Method

A wide variety of plasma phenomena have been investigated during the past decades using the particle-in-cell/Monte Carlo collisions (PIC/MCC) method. As an important component of the PIC/MCC method, solving Poisson’s equation is crucial for the accuracy and efficiency of calculations. Different acceleration techniques for solving finite difference discretization Poisson’s equation are investigated and compared, including direct method, iterative method, multigrid (MG) method, parallel computing and inherited initial value. The charge density distribution with a known analytical solution is used to validate the algorithm and code. The optimal relaxation factor for the successive over-relaxation (SOR) method in 2D Poisson’s equation with unequal grid node numbers in different dimensions is derived, which is only related to the dimension with the largest grid number. Although there will be a ‘more optimal’ relaxation factor deviated from in some simulation cases, selecting the optimal relaxation factor derived always leads to a not slow solving speed. However, when SOR is used in MG for smoothing, the optimal relaxation factor will shift to 0.5–1.2 from the theoretical optimal value derived with the increase of MG levels. By comparing the convergence order under different relaxation factors and MG levels, the suitable MG level is proposed as log2[min(N x , N y )]−2. Combining the optimal SOR relaxation factor, MG, parallel computing and inherited initial values, the computational cost may decrease by 5 orders of magnitude than that by the simple Gaussian elimination (GE). Based on the optimal acceleration techniques mentioned above, a benchmark simulation case electron cyclotron drift instability (ECDI) in magnetized plasmas was run to further validate the developed PIC/MCC code. The distributions of electric field in the x-direction, electron density and electron temperature are all consistent with the literatures. This paper provides a reference for the acceleration strategy selection for solving Poisson’s equation quickly in plasma simulations.

Numerical solution of the boundary value problems for the biharmonic equations via quasiseparable representations

AbstractThe paper incorporates new methods of numerical linear algebra for the approximation of the biharmonic equation with potential, namely, numerical solution of the Dirichlet problem for $$ \left( \frac{d}{dx}\right) ^4u(x)+c(x)u(x)=\phi (x),\quad 0<x<1. $$ d dx 4 u ( x ) + c ( x ) u ( x ) = ϕ ( x ) , 0 < x < 1 . High-order discrete finite difference operators are presented, constructed on the basis of discrete Hermitian derivatives, and the associated Discrete Biharmonic Operator (DBO). It is shown that the matrices associated with the discrete operator belong to a class of quasiseparable matrices of low rank matrices. The application of quasiseparable representation of rank structured matrices yields fast and stable algorithm for variable potentials c(x). Numerical examples corroborate the claim of high order accuracy of the algorithm, with optimal complexity O(N).

PGD in thermal transient problems with a moving heat source: A sensitivity study on factors affecting accuracy and efficiency

AbstractThermal transient problems, essential for modeling applications like welding and additive metal manufacturing, are characterized by a dynamic evolution of temperature. Accurately simulating these phenomena is often computationally expensive, thus limiting their applications, for example for model parameter estimation or online process control. Model order reduction, a solution to preserve the accuracy while reducing the computation time, is explored. This article addresses challenges in developing reduced order models using the proper generalized decomposition (PGD) for transient thermal problems with a specific treatment of the moving heat source within the reduced model. Factors affecting accuracy, convergence, and computational cost, such as discretization methods (finite element and finite difference), a dimensionless formulation, the size of the heat source, and the inclusion of material parameters as additional PGD variables are examined across progressively complex examples. The results demonstrate the influence of these factors on the PGD model's performance and emphasize the importance of their consideration when implementing such models. For thermal example, it is demonstrated that a PGD model with a finite difference discretization in time, a dimensionless representation, a mapping for a moving heat source, and a spatial domain non‐separation yields the best approximation to the full order model.

3D numerical simulation of frequency-domain elastic wave equation with a compact second-order scheme

Finite-difference frequency-domain (FDFD) modeling has gained popularity in recent years. However, its main drawback lies in the need to solve large sparse linear systems, which is computationally expensive, especially for 3D elastic wave simulations. The existing elastic wave optimal schemes improve the simulation accuracy at the cost of using large-sized stencils or high-order finite-difference (FD) operators, and most of them only consider equal grid spacing and solid media. To address these issues, we have proposed a compact second-order scheme for simulating 3D elastic wave propagation. Starting from the second-order FD discretization of the first-order velocity-stress expression, the compact discrete scheme of the second-order elastic wave equation is obtained by using a parsimonious staggered-grid strategy to eliminate the stress variables. In this way, the scheme preserves the ability to handle high-contrast velocity and liquid-solid media, and has a more compact and sparse impedance matrix than existing schemes, which helps save computational costs. We also use an antilumped mass strategy to further improve the modeling accuracy. Dispersion analysis indicates that our compact second-order scheme has better accuracy than the classic second-order scheme, and even better accuracy than the average derivative 27-point scheme at large Poisson's ratio. Compared with the existing 3D elastic wave schemes, the compact second-order scheme has advantages in computational cost for the same grid. Finally, we employ these schemes to implement 3D numerical simulations, validating the effectiveness of our proposed scheme.

On iterative methods based on Sherman-Morrison-Woodbury splitting

We consider a new splitting based on the Sherman-Morrison-Woodbury formula, which is particularly effective with iterative methods for the numerical solution of large linear systems and especially for systems involving matrices that are perturbations of circulant or block circulant matrices. Such matrices typically arise in the discretization of differential equations using finite element or finite difference methods. We prove the convergence of the new iteration without making any assumptions regarding the symmetry or diagonal-dominance of the matrix, which are limiting factors for most classical iterative methods.To illustrate the efficacy of the new iteration we present various applications. These include extensions of the new iteration to block matrices that arise in certain saddle point problems as well as two-dimensional finite difference discretizations. The new method exhibits fast convergence in all of the test cases we used. It has minimal storage requirements, straightforward implementation and compatibility with nearly circulant matrices via the Fast Fourier Transform. Remarkably, the new method was tested against very large matrices demonstrating extremely fast convergence. For these reasons it can be a valuable tool for the solution of various finite elements and finite differences discretizations of differential equations.

Non-Overlapping Domain Decomposition for 1D Optimal Control Problems Governed by Time-Fractional Diffusion Equations on Coupled Domains: Optimality System and Virtual Controls

We consider a non-overlapping domain decomposition method for optimal control problems of the tracking type governed by time-fractional diffusion equations in one space dimension, where the fractional time derivative is considered in the Caputo sense. We concentrate on a transmission problem defined on two adjacent intervals, where at the interface we introduce an iterative non-overlapping domain decomposition in the spirit of P.L. Lions for the corresponding first-order optimality system, such that the optimality system corresponding to the optimal control problem on the entire domain is iteratively decomposed into two systems on the respective sub-domains; this approach can be framed as first optimize, then decompose. We show that the iteration involving the states and adjoint states converges in the appropriate spaces. Moreover, we show that the decomposed systems on the sub-domain can in turn be interpreted as optimality systems of so-called virtual control problems on the sub-domains. Using this property, we are able to solve the original optimal control problem by an iterative solution of optimal control problems on the sub-domains. This approach can be framed as first decompose, then optimize. We provide a mathematical analysis of the problems as well as a numerical finite difference discretization using the L1-method with respect to the Caputo derivative, along with two examples in order to verify the method.

Mancha3D Code: Multipurpose Advanced Nonideal MHD Code for High-Resolution Simulations in Astrophysics

The Mancha3D code is a versatile tool for numerical simulations of magnetohydrodynamic (MHD) processes in solar/stellar atmospheres. The code includes nonideal physics derived from plasma partial ionization, a realistic equation of state and radiative transfer, which allows performing high-quality realistic simulations of magnetoconvection, as well as idealized simulations of particular processes, such as wave propagation, instabilities or energetic events. The paper summarizes the equations and methods used in the Mancha3D (Multifluid (-purpose -physics -dimensional) Advanced Non-ideal MHD Code for High resolution simulations in Astrophysics 3D) code. It also describes its numerical stability and parallel performance and efficiency. The code is based on a finite difference discretization and a memory-saving Runge–Kutta (RK) scheme. It handles nonideal effects through super-time-stepping and Hall diffusion schemes, and takes into account thermal conduction by solving an additional hyperbolic equation for the heat flux. The code is easily configurable to perform different kinds of simulations. Several examples of the code usage are given. It is demonstrated that splitting variables into equilibrium and perturbation parts is essential for simulations of wave propagation in a static background. A perfectly matched layer (PML) boundary condition built into the code greatly facilitates a nonreflective open boundary implementation. Spatial filtering is an important numerical remedy to eliminate grid-size perturbations enhancing the code stability. Parallel performance analysis reveals that the code is strongly memory bound, which is a natural consequence of the numerical techniques used, such as split variables and PML boundary conditions. Both strong and weak scalings show adequate performance up to several thousands of processors (CPUs).

A positivity-preserving, linear, energy stable and convergent numerical scheme for the Poisson–Nernst–Planck (PNP) system

This article focuses on the convergence analysis for a fully discrete finite difference scheme for the time-dependent Poisson–Nernst–Planck system. The numerical scheme, a three-level linearized finite difference algorithm based on the reformulation of the Nernst–Planck equations, which was developed in [J. Sci. Comput. 81(2019),436–458]. The positivity-preserving property and energy stability were theoretically established. In this paper, we rigorously prove first-order convergence in time and second-order convergence in space for the numerical scheme in the ℓ∞(0,T;ℓ2)∩ℓ2(0,T;Hh1) norm under the condition that time step size is linearly proportional to spatial mesh size, in which the natural property of the exponential terms gives a lot of challenges. Moreover, the higher-order asymptotic expansion (up to third-order temporal accuracy and fourth-order spatial accuracy) has to be involved, due to the leading local truncation error will not be enough to recover an ℓ∞ bound for ion concentrations n and p. To our knowledge, this scheme will be the first linear decoupled algorithm to combine the following theoretical properties for the PNP system: ion concentration positivity preserving, unconditionally energy stability and optimal rate convergence. Numerical results are shown to be consistent with theoretical analysis.

Legendre wavelet collocation method for investigating thermo-mechanical responses on biological tissue during laser irradiation

Objective of the present work is to propose a numerical approach based on the Legendre wavelet to address a problem of biological tissue in the presence of superficial cancer and analyze the thermo-mechanical behavior during laser irradiation. The model of thermoelasticity with variable thermal conductivity is considered to investigate transient thermoelastic coupling response in the framework of Moore–Gibson–Thompson (MGT) bioheat transfer. The Kirchhoff mapping is first performed to linearize the nonlinear governing equations due to variable thermal properties, then finite difference discretization is applied for the time domain and subsequently space domain is approximated using Legendre wavelets. The collocation points are taken to transform the system of partial differential equations into a set of algebraic equations, from which the solution can be determined. A thorough parametric analysis is carried out to examine the effects of some material and geometrical parameters on the behavior of various fields such as the displacement, temperature, stress and the tumor-normal tissue interface. The outcomes of the present work are graphically shown and compared to the predictions of other existing models. The advantages of the current numerical procedure include its simplicity, effectiveness and low computational cost even with the few collocation points.

Solving inverse scattering problems via reduced-order model embedding procedures

We present a reduced-order model (ROM) methodology for inverse scattering problems in which the ROMs are data-driven, i.e. they are constructed directly from data gathered by sensors. Moreover, the entries of the ROM contain localised information about the coefficients of the wave equation. We solve the inverse problem by embedding the ROM in physical space. Such an approach is also followed in the theory of ‘optimal grids,’ where the ROMs are interpreted as two-point finite-difference discretisations of an underlying set of equations of a first-order continuous system on this special grid. Here, we extend this line of work to wave equations and introduce a new embedding technique, which we call Krein embedding, since it is inspired by Krein’s seminal work on vibrations of a string. In this embedding approach, an adaptive grid and a set of medium parameters can be directly extracted from a ROM and we show that several limitations of optimal grid embeddings can be avoided. Furthermore, we show how Krein embedding is connected to classical optimal grid embedding and that convergence results for optimal grids can be extended to this novel embedding approach. Finally, we also briefly discuss Krein embedding for open domains, that is, semi-infinite domains that extend to infinity in one direction.

A novel numerical note on the enhanced thermal features of water-ethylene glycol mixture due to hybrid nanoparticles (MnZnFe2O4 – Ag) over a magnetized stretching surface

In this article, we presented a numerical solution of steady-state magnetohydrodynamic (MHD) flow of an incompressible, electrically conducting Casson hybrid nanofluid through a stretching surface under the incitement of viscous dissipation, heat generation and induced magnetism. A viscous mixture of water and ethylene glycol ( W : E G = 60 % : 40 % ) is considered as a base fluid. Enhanced thermal aspects of the working fluid due to suspended manganese zinc ferrite ( MnZnFe 2 O 4 ) and silver ( Ag ) nanoparticles are analyzed by successive over relaxation iterative method. A self-similar procedure is adopted to obtain dimensionless system of governing equations. Order reduction and finite differences discretization of nonlinear coupled differential equations lead to the approximate solution of proposed problem. Impacts of various dominant parameters on the hybridized Casson fluid are examined through the graphs sketched for the functions of velocity, temperature, and induced magnetism. Skin friction coefficient and thermal transfer characteristics at the stretching surface have been simultaneously observed via tabular data prepared for pure nanofluid ( Ag-WEG ) as well as for hybrid nanofluid ( MnZnFe 2 O 4 -Ag/WEG ). It is perceived that induced magnetic field tends to decline the Nusselt number while exhibits a significant growth in the surface drag. Moreover, a significant enhancement in temperature of working fluid has been observed due to Eckert number and heat generation parameter. This decisive role of enhanced thermal features will help to improve industrial growth, waste-water treatment, and automobiles engine efficiency.

Numerical approximations and asymptotic limits of some nonlinear problems

In the present work, a numerical approach is dedicated to the approximation to the solutions of a time-independent nonlinear Schr¨odinger equation in a mixed case provided with numerical tests on the asymptotic limits of the solution according to some parameters. A finite difference discretization with calibrations is applied leading to a quasi-linear algebraic system. where the its solvability is investigated as well as its stability and convergence via Von Neumann method. Some numerical experiments are developed to validate the result, and to test the effect of some parameters on the asymptotic limit of the problem.