Coupling of the Crank–Nicolson scheme and localized meshless technique for viscoelastic wave model in fluid flow

Abstract

This paper proposes an efficient localized meshless technique for approximating the viscoelastic wave model. This model is a significant methodology to explain wave propagation in solids modeled with a wide collection of viscoelastic laws. In the first method, a difference scheme with the second-order accuracy is implemented to obtain a semi-discrete scheme. Then, a localized radial basis function partition of unity scheme is adopted to get a full-discrete scheme. This localization technique consists of decomposing the initial domain into several sub-domains and constructing a local radial basis function approximation over every sub-domain. A well-conditioned resulting linear system and a low computational burden are the main merits of this technique compared to global collocation methods. Further, the stability and convergence analysis of the temporal discretization scheme are deduced using discrete energy method. Numerical results are shown to validate the accuracy and effectiveness of the proposed method.

Introduction

Wave propagation along with dispersion and attenuation in real-world media can be modeled using the viscoelastic wave equation. Numerous factors, such as fluid saturation and rock porosity, may lead to viscoelastic effects [1], [2]. Therefore, the detection of these viscoelastic impacts can directly apply to reservoir prediction and lithologic identification. Seismic migration and inversion techniques based on the wave equation and accounting for complex wave propagation processes are becoming more popular. The basis of these techniques is seismic wave modeling. As a result of intrinsic attenuation in real-world media, viscoelastic and viscoacoustic wave models can model the wave propagation more accurately than elastic and acoustic wave models. In this paper, we study a numerical approach for the approximated solution of the time viscoelastic wave equation (TVWE) as:

$$\frac{{\partial }^{2}u(\mathbf{x},t)}{\partial t^2}-\gamma \frac{\partial \Delta u(\mathbf{x},t)}{\partial t}-\sigma \Delta u(\mathbf{x},t)=f(\mathbf{x},t),(\mathbf{x},t)\in \Omega \times J,\mathbf{x}=(x,y).$$

The initial conditions are considered as

$$u(\mathbf{x},0)=\rho \left(\mathbf{x}\right),\mathbf{x}\in \Omega ,$$$$\frac{\partial u(\mathbf{x},0)}{\partial t}=\varpi \left(\mathbf{x}\right),\mathbf{x}\in \Omega ,$$

together with the boundary conditions

$$u(\mathbf{x},t)=h(\mathbf{x},t),(\mathbf{x},t)\in \partial \Omega \times \bar{J},$$

in which $\sigma$ and $\gamma$ are two given positive constants,$\Delta$ is the Laplacian operator, $\Omega \subset {\mathbb{R}}^d$ represents a bounded convex polyhedron region with boundary $\partial \Omega$, and $J=(0,T]$ denotes the temporal interval with $0<T<\infty$ ($\bar{J}$ is the closure of $J$). The forcing term $f(\mathbf{x},t)$, the two initial value functions $\rho \left(\mathbf{x}\right)$ and $\varpi \left(\mathbf{x}\right)$, and the boundary value function $h(\mathbf{x},t)$ denote prescribed functions with sufficient smoothness.

To our best knowledge, only a few numerical approaches have been presented for simulating the TVWE. Li et al. [3] established a space–time continuous finite element (STCFE) method, while Yuan [4] adopted the finite difference (FD) approach. Luo et al. [5] applied an optimized splitting positive definite mixed finite element (FE) extrapolation scheme based on proper orthogonal decomposition (POD) approach, whereas Luo and Xia [6], [7], [8] used the optimized FE extrapolating method and POD based on the FD method. Jin and Luo [9] adopted the Crank–Nicolson collocation spectral (CNCS) technique and Zhao et al. [10] proposed a space–time continuous Galerkin method with mesh modification. Oruç [11] formulated the barycentric rational interpolation and local radial basis function (BRI and LRBF, respectively) techniques based on the explicit method in time direction.

As seen from the articles cited above, researchers investigate the TVWE mostly using mesh-based techniques such as the FE and FD methods. To use the FE method in solving the considered problems, a mesh must be generated. However, it might be difficult to generate this mesh, especially over irregular solution domains. Moreover, it is required to calculate integrals over single elements in the FE, which might consume a lot of time. Nevertheless, in the FD, the computational burden might be lower, and implementing the technique may be less difficult compared to the FE; however, the FD cannot readily handle irregular domains. To avoid generating sophisticated meshes or calculating integrals over single elements and to easily handle irregular domains, researchers have introduced meshless techniques. These techniques are frequently used for numerically solving partial differential equations (PDEs).

Recently, considerable effort has been made to investigate meshfree techniques, also known as meshless methods in the literature. The aim of meshless techniques is to at least eliminate the mesh structure and to approximate the entire solution utilizing the nodes as a quasirandom or scattered set of nodes rather than using nodes in a grid/element-based discretization [12]. Indeed, these methods do not require to generate the underlying mesh for approximation [13], [14]. There is a class of meshless techniques that are based on RBF collocation approaches, that use linear combinations of the translates of one single basic univariate function for collocation [15], [16]. In recent years, meshfree RBF-based techniques have become popular in the scientific community for several reasons. Indeed, the hybridization of the RBF and interpolation leads to very good features in complex geometrical and high-dimension problems. The main features of these techniques are meshfree nature, spectral convergence, flexibility in handling sophisticated geometries, and simple extension to problems of several dimensions. An RBF depends only on the inter-center distances; hence, one typically uses RBF interpolants for approximating given functions. In addition, the RBF is used as kernel in support vector classification [17], [18]. Nevertheless, using globally-supported basis functions results in collocation matrices that are fully populated. These matrices become more ill-conditioned and computationally costly with an increase in the size of the dataset. Such constraints have motivated scientists and researchers to study various methods for restricting the influence domain of basis functions, hence reducing the issues of computational burden and numerical conditioning, whereas keeping the flexibility and efficiency of the previous formulation. A number of excellent books in this area have been written by various authors, who investigated the different characteristics and potential applications of the RBF [19], [20], [21], [22]. Wu [23] and Wendland [24] created compactly supported technique and proposed a group of positive definite polynomial functions for different continuity orders and spatial dimensions. Implementing these techniques is straightforward; however, they are not significantly flexible, and one can show that large support scales are usually needed for obtaining reasonable convergence rates and accuracy [25]. An alternate strategy is performing a direct domain decomposition of the RBF method, in which the initial domain is partitioned into a number of interrelating sub-domains that can or cannot overlap based on the formulation. Although domain decomposition tends to be attractive in principle, defining appropriate sub-domains can be very complex and highly problem-dependent. Here, we introduce the two main approaches to reduce the computational cost of RBF-based methods through localization. The first strategy, called the RBF-FD technique, is the natural generalization of the classical grid-based FD and is constructed by combining the RBF and FD. As far as we know, the RBF-FD was first published by Tolstykh [26] and investigated by [27], [28]. In fact, the local RBF-FD is more flexible since it works on various complex-shaped domains with any type of distribution node [29], [30], [31], [32], [33]. The second strategy, which is considered in the present work, is the RBF partition of unity method (RBF-PUM). This method involves the decomposition of the domain into an appropriate number of overlapping patches or sub-domains that form a cover of the initial domain. After creating a local RBF approximant on every sub-domain, one weights these approximations together using the compactly supported RBFs (CSRBFs) by the family of the Wendland functions in order to create the global fit. Based on this method, a large problem in different dimensions can be decomposed into numerous small problems. Hence, we will be able to work with a large number of scattered data-nodes. Moreover, one can use the convergence characteristics of the local approximations, while enforcing local couplings between approximations on various sub-domains by using the PU framework. Babuška and Melenk [34], [35] originally proposed the PUM for approximating PDEs and suggested the combination of the RBF approximation with the PUM. Subsequently, Wendland [36] proposed and examined a combination of the PUM and CSRBFs for the purpose of interpolation.

Recently, Cavoretto et al. [37], [38], [39], [40] developed a number of efficient algorithms based on the RBF-PUM interpolation. Cavoretto [41] studied the performance of product-type interpolants when they are applied as local approximants in the RBF-PUM. Larsson et al. [42], [43] published some works concerning the RBF-PUM for PDEs. Cavoretto et al. [44] proposed a new stable and accurate approximation technique based on the RBF-PUM. Mollapourasl et al. [45] adopted the RBF-PUM to price American options under Heston’s stochastic volatility model. Cavoretto [46] employed a fast algorithm for two-dimensional interpolation of large scattered data sets based on the RBF-PUM. Mollapourasl et al. [47] applied the RBF-PUM for option pricing problems in jump–diffusion model, while Fereshtian et al. [48] proposed for the price American and European options under Lévy model. Garmanjani et al. [49] used the RBF-PUM collocation scheme based on FD to approximate the convection–diffusion and pseudo-parabolic problems. Gholampour et al. [50], [51] presented a stable RBF-PUM for solving elliptic interface problems. Cavoretto and De Rossi [52] presented the RBF-PUM for solving time-independent elliptic PDEs such as Poisson problems. De Marchi et al. [53] performed a local computation using the PUM of rational RBF interpolant. They studied the well-posedness of the problem and provided the error bounds. Esmaeilbeig et al. [54] solved the stochastic PDEs based on the RBF-PUM, while Darani [55] applied the localized RBF-PUM for solving the Klein–Gordon equation. Cavorettoet al. [56] applied the RBF-PUM as a highly adaptive auxiliary tool for graph signal processing. Mirzaei [57] proposed a novel localized RBF-PUM for approximating boundary and initial–boundary value problems. Cavoretto [58] advanced a new adaptive algorithm for bivariate interpolation of large scattered data points through the RBF-PUM. Nikan and Avazzadeh [59] employed the RBF-PUM collocation scheme based on FD for the Sobolev equation.

The main objective of this work is to propose the localized RBF-PUM for finding the approximate solution of the TVWE. The proposed method is performed by using a Crank–Nicolson FD technique for the time discretization and then a spatial discretization is deduced through the localized RBF-PUM. This localized technique enables reducing the computational burden while keeping high accuracy. The major benefits of the RBF-PUM technique are its ability to maintain geometrical flexibility in problems of high dimensions and to make adaptive approximation easy. After that, we use an existing algorithm in the literature to select an optimal shape parameter that has a palpable effect on coefficient matrix. The unconditional stability and convergence of the time-discretized scheme are proven theoretically and verified numerically by doing its error analysis. Moreover, the obtained results are compared with those of other techniques in existing literature showing the accuracy and efficiency of the proposed approach.

The framework of this work has been organized as follows. Section 2 adopts a time-stepping scheme including second-order accuracy to derive the semi-discrete algorithm. Furthermore, the time-discrete approach in terms of the unconditional stability and convergence issues is analyzed. Section 3 illustrates the description of the RBF-PUM for space discretization and its numerical implementation. Section 4 implements the proposed method for three test problems that exhibit its efficiency and clarify the theoretical results. Finally, Section 5 allocates to the main conclusions of this work.