Perfectly matched layers for transient elastodynamics of unbounded domains

One approach to the numerical solution of a wave equation on an unbounded domain uses a bounded domain surrounded by an absorbing boundary or layer that absorbs waves propagating outward from the bounded domain. A perfectly matched layer (PML) is an unphysical absorbing layer model for linear wave equations that absorbs, almost perfectly, outgoing waves of all non‐tangential angles‐of‐incidence and of all non‐zero frequencies. In a recent work [Computer Methods in Applied Mechanics and Engineering 2003; 192: 1337–1375], the authors presented, inter alia, time‐harmonic governing equations of PMLs for anti‐plane and for plane‐strain motion of (visco‐) elastic media. This paper presents (a) corresponding time‐domain, displacement‐based governing equations of these PMLs and (b) displacement‐based finite element implementations of these equations, suitable for direct transient analysis. The finite element implementation of the anti‐plane PML is found to be symmetric, whereas that of the plane‐strain PML is not. Numerical results are presented for the anti‐plane motion of a semi‐infinite layer on a rigid base, and for the classical soil–structure interaction problems of a rigid strip‐footing on (i) a half‐plane, (ii) a layer on a half‐plane, and (iii) a layer on a rigid base. These results demonstrate the high accuracy achievable by PML models even with small bounded domains. Copyright © 2004 John Wiley & Sons, Ltd.     

Boundary element‐free method for fracture analysis of 2‐D anisotropic piezoelectric solids

This paper considers a 2‐D fracture analysis of anisotropic piezoelectric solids by a boundary element‐free method. A traction boundary integral equation (BIE) that only involves the singular terms of order 1/r is first derived using integration by parts. New variables, namely, the tangential derivative of the extended displacement (the extended displacement density) for the general boundary and the tangential derivative of the extended crack opening displacement (the extended displacement dislocation density), are introduced to the equation so that solution to curved crack problems is possible. This resulted equation can be directly applied to general boundary and crack surface, and no separate treatments are necessary for the upper and lower surfaces of the crack. The extended displacement dislocation densities on the crack surface are expressed as the product of the characteristic terms and unknown weight functions, and the unknown weight functions are modelled using the moving least‐squares (MLS) approximation. The numerical scheme of the boundary element‐free method is established, and an effective numerical procedure is adopted to evaluate the singular integrals. The extended ‘stress intensity factors’ (SIFs) are computed for some selected example problems that contain straight or curved cracks, and good numerical results are obtained. Copyright © 2006 John Wiley & Sons, Ltd.     

Solving the advection–diffusion equation on unstructured meshes with discontinuous/mixed finite elements and a local time stepping procedure

Explicit schemes are known to provide less numerical diffusion in solving the advection–diffusion equation, especially for advection‐dominated problems. Traditional explicit schemes use fixed time steps restricted by the global CFL condition in order to guarantee stability. This is known to slow down the computation especially for heterogeneous domains and/or unstructured meshes. To avoid this problem, local time stepping procedures where the time step is allowed to vary spatially in order to satisfy a local CFL condition have been developed. In this paper, a local time stepping approach is used with a numerical model based on discontinuous Galerkin/mixed finite element methods to solve the advection–diffusion equation. The developments are detailed for general unstructured triangular meshes. Numerical experiments are performed to show the efficiency of the numerical model for the simulation of (i) the transport of a solute on highly unstructured meshes and (ii) density‐driven flow, where the velocity field changes at each time step. The model gives stable results with significant reduction of the computational cost especially for the non‐linear problem. Moreover, numerical diffusion is also reduced for highly advective problems. Copyright © 2009 John Wiley & Sons, Ltd.     

High‐order finite volume schemes for the advection–diffusion equation

A high‐order finite volume method based on piecewise interpolant polynomials is proposed to discretize spatially the one‐dimensional and two‐dimensional advection–diffusion equation. Evolution equations for the mean values of each control volume are integrated in time by a classical fourth‐order Runge–Kutta. Since our work focuses on the behaviour of the spatial discretization, the time step is chosen small enough to neglect the time integration error. Two‐dimensional interpolants are built by means of one‐dimensional interpolants. It is shown that when the degree of the one‐dimensional interpolant q is odd, the proper selection of a fixed stencil gives rise to centred schemes of order q+1. In order not to lose precision due to the change of stencil near boundaries, the degree of the interpolants close to boundaries is raised to q+1. Four test cases with small values of diffusion are integrated with high‐order methods. It is shown that the spatial discretization of the advection–diffusion equation with periodic boundary conditions leads to normal discretization matrices, and asymptotic stability must be assured to bound the spatial discretization error. Once the asymptotic stability is assured by means of the spectra of the discretization matrix, the spatial error is of the order of the truncation error. However, it is shown that the discretization of the advection–diffusion equation with arbitrary boundary conditions gives rise to non‐normal matrices. If asymptotic stability is assured, the spatial order of steady solutions is of the order of the truncation error. But, for transient processes, the order of the spatial error is determined by both the truncation error and the norm of the exponential matrix of the spatial discretization. The use of the pseudospectra of the discretization matrix is proposed as a valuable tool to analyse the transient error of the numerical solution. Copyright © 2001 John Wiley & Sons, Ltd.     

An improved continued‐fraction‐based high‐order transmitting boundary for time‐domain analyses in unbounded domains

A high‐order local transmitting boundary to model the propagation of acoustic or elastic, scalar or vector‐valued waves in unbounded domains of arbitrary geometry is proposed. It is based on an improved continued‐fraction solution of the dynamic stiffness matrix of an unbounded medium. The coefficient matrices of the continued‐fraction expansion are determined recursively from the scaled boundary finite element equation in dynamic stiffness. They are normalised using a matrix‐valued scaling factor, which is chosen such that the robustness of the numerical procedure is improved. The resulting continued‐fraction solution is suitable for systems with many DOFs. It converges over the whole frequency range with increasing order of expansion and leads to numerically more robust formulations in the frequency domain and time domain for arbitrarily high orders of approximation and large‐scale systems. Introducing auxiliary variables, the continued‐fraction solution is expressed as a system of linear equations in iω in the frequency domain. In the time domain, this corresponds to an equation of motion with symmetric, banded and frequency‐independent coefficient matrices. It can be coupled seamlessly with finite elements. Standard procedures in structural dynamics are directly applicable in the frequency and time domains. Analytical and numerical examples demonstrate the superiority of the proposed method to an existing approach and its suitability for time‐domain simulations of large‐scale systems. Copyright © 2011 John Wiley & Sons, Ltd.     

Implementation of meshless LBIE method to the 2D non‐linear SG problem

In this paper the meshless local boundary integral equation (LBIE) method for numerically solving the non‐linear two‐dimensional sine‐Gordon (SG) equation is developed. The method is based on the LBIE with moving least‐squares (MLS) approximation. For the MLS, nodal points spread over the analyzed domain are utilized to approximate the interior and boundary variables. The approximation functions are constructed entirely using a set of scattered nodes, and no element or connectivity of the nodes is needed for either the interpolation or the integration purposes. A time‐stepping method is employed to deal with the time derivative and a simple predictor–corrector scheme is performed to eliminate the non‐linearity. A brief discussion is outlined for numerical integrations in the proposed algorithm. Some examples involving line and ring solitons are demonstrated and the conservation of energy in undamped SG equation is investigated. The final numerical results confirm the ability of method to deal with the unsteady non‐linear problems in large domains. Copyright © 2009 John Wiley & Sons, Ltd.     

Efficient solution of time‐domain boundary integral equations arising in sound‐hard scattering

We consider the efficient numerical solution of the three‐dimensional wave equation with Neumann boundary conditions via time‐domain boundary integral equations. A space‐time Galerkin method with C^∞‐smooth, compactly supported basis functions in time and piecewise polynomial basis functions in space is employed. We discuss the structure of the system matrix and its efficient parallel assembly. Different preconditioning strategies for the solution of the arising systems with block Hessenberg matrices are proposed and investigated numerically. Furthermore, a C++ implementation parallelized by OpenMP and MPI in shared and distributed memory, respectively, is presented. The code is part of the boundary element library BEM4I. Results of numerical experiments including convergence and scalability tests up to a thousand cores on a cluster are provided. The presented implementation shows good parallel scalability of the system matrix assembly. Moreover, the proposed algebraic preconditioner in combination with the FGMRES solver leads to a significant reduction of the computational time. Copyright © 2015 John Wiley & Sons, Ltd.     

Indirect T‐Trefftz and F‐Trefftz methods for solving boundary value problem of Poisson equation

Abstract

The Poisson equation can be solved by first finding a particular solution and then solving the resulting Laplace equation. In this paper, a computational procedure based on the Trefftz method is developed to solve the Poisson equation for two‐dimensional domains. The radial basis function approach is used to find an approximate particular solution for the Poisson equation. Then, two kinds of Trefftz methods, the T‐Trefftz method and F‐Trefftz method, are adopted to solve the resulting Laplace equation. In order to deal with the possible ill‐posed behaviors existing in the Trefftz methods, the truncated singular value decomposition method and L‐curve concept are both employed. The Poisson equation of the type, ?² u = f(x, u), in which x is the position and u is the dependent variable, is solved by the iterative procedure. Numerical examples are provided to show the validity of the proposed numerical methods and some interesting phenomena are carefully discussed while solving the Helmholtz equation as a Poisson equation. It is concluded that the F‐Trefftz method can deal with a multiply connected domain with genus p(p > 1) while the T‐Trefftz method can only deal with a multiply connected domain with genus 1 if the domain partition technique is not adopted.     

Frequency‐domain analysis of time‐integration methods for semidiscrete finite element equations—part I: Parabolic problems

Time‐integration methods for semidiscrete equations emanating from parabolic differential equations are analysed in the frequency domain. The discrete‐time transfer functions of three popular methods are derived, and subsequently the forced response characteristics of single modes are studied in the frequency domain. To enable consistent comparison of the frequency responses of different algorithms, three characteristic numbers are identified. Frequency responses and L₂‐norms of the phase and magnitude errors are compared for the three time‐integration algorithms. The examples demonstrate that frequency‐domain analysis provides substantial insight into the time‐domain properties of time‐integration algorithms. Copyright © 2001 John Wiley & Sons, Ltd.     

A combined FIC‐TDG finite element approach for the numerical solution of coupled advection–diffusion–reaction equations with application to a bioregulatory model for bone fracture healing

Numerical schemes for the approximative solution of advection–diffusion–reaction equations are often flawed because of spurious oscillations, caused by steep gradients or dominant advection or reaction. In addition, for strong coupled nonlinear processes, which may be described by a set of hyperbolic PDEs, established time stepping schemes lack either accuracy or stability to provide a reliable solution. In this contribution, an advanced numerical scheme for this class of problems is suggested by combining sophisticated stabilization techniques, namely the finite calculus (FIC‐FEM) scheme introduced by Oñate et al. with time‐discontinuous Galerkin (TDG) methods. Whereas the former one provides a stabilization technique for the numerical treatment of steep gradients for advection‐dominated problems, the latter ensures reliable solutions with regard to the temporal evolution. A brief theoretical outline on the superior behavior of both approaches will be presented and underlined with related computational tests. The performance of the suggested FIC‐TDG finite element approach will be discussed exemplarily on a bioregulatory model for bone fracture healing proposed by Geris et al., which consists of at least 12 coupled hyperbolic evolution equations. Copyright © 2012 John Wiley & Sons, Ltd.     

GPU accelerated lattice Boltzmann model for shallow water flow and mass transport

A lattice Boltzmann method (LBM) for solving the shallow water equations (SWEs) and the advection–dispersion equation is developed and implemented on graphics processing unit (GPU)‐based architectures. A generalized lattice Boltzmann equation (GLBE) with a multiple‐relaxation‐time (MRT) collision method is used to simulate shallow water flow. A two‐relaxation‐time (TRT) method with two speed‐of‐sound techniques is used to solve the advection–dispersion equation. The proposed LBM is implemented to an NVIDIA ^® Computing Processor in a single GPU workstation. GPU computing is performed using the Jacket GPU engine for MATLAB ^® and CUDA. In the numerical examples, the MRT‐LBM model and the TRT‐LBM model are verified and show excellent agreement to exact solutions. The MRT outperforms the single‐relaxation‐time (SRT) collision operator in terms of stability and accuracy when the SRT parameter is close to the stability limit of 0.5. Mass transport with velocity‐dependent dispersion in shallow water flow is simulated by combining the MRT‐LBM model and the TRT‐LBM model. GPU performance with CUDA code shows an order of magnitude higher than MATLAB‐Jacket code. Moreover, the GPU parallel performance increases as the grid size increases. The results indicate the promise of the GPU‐accelerated LBM for modeling mass transport phenomena in shallow water flows. Copyright © 2010 John Wiley & Sons, Ltd.     

An adaptive fixed mesh method for modelling and simulating of unsteady two‐phase flows with sharp physical interfaces

We present in this paper a new computational method for simulation of two‐phase flow problems with moving boundaries and sharp physical interfaces. An adaptive interface‐capturing technique (ICT) of the Eulerian type is developed for capturing the motion of the interfaces (free surfaces) in an unsteady flow state. The adaptive method is mainly based on the relative boundary conditions of the zero pressure head, at which the interface is corresponding to a free surface boundary. The definition of the free surface boundary condition is used as a marker for identifying the position of the interface (free surface) in the two‐phase flow problems. An initial‐value‐problem (IVP) partial differential equation (PDE) is derived from the dynamic conditions of the interface, and it is designed to govern the motion of the interface in time. In this adaptive technique, the Navier–Stokes equations written for two incompressible fluids together with the IVP are solved numerically over the flow domain. An adaptive mass conservation algorithm is constructed to govern the continuum of the fluid. The finite element method (FEM) is used for the spatial discretization and a fully coupled implicit time integration method is applied for the advancement in time. FE‐stabilization techniques are added to the standard formulation of the discretization, which possess good stability and accuracy properties for the numerical solution. The adaptive technique is tested in simulation of some numerical examples. With the test problems presented here, we demonstrated that the adaptive technique is a simple tool for modelling and computation of complex motion of sharp physical interfaces in convection–advection‐dominated flow problems. We also demonstrated that the IVP and the evolution of the interface function are coupled explicitly and implicitly to the system of the computed unknowns in the flow domain. Copyright © 2005 John Wiley & Sons, Ltd.     

On the efficient computation of the stress components near a closed boundary in plane elasticity by using classical complex boundary integral equations

Complex boundary integral equations (Fredholm‐type regular or Cauchy‐type singular or even Hadamard–Mangler‐type hypersingular) have been used for the numerical solution of general plane isotropic elasticity problems. The related Muskhelishvili and, particularly, Lauricella–Sherman equations are famous in the literature, but several more extensions of the Lauricella–Sherman equations have also been proposed. In this paper it is just mentioned that the stress and displacement components can be very accurately computed near either external or internal simple closed boundaries (for anyone of the above equations: regular or singular or hypersingular, but with a prerequisite their actual numerical solution) through the appropriate use of the even more classical elementary Cauchy theorem in complex analysis. This approach has been already used for the accurate numerical computation of analytic functions and their derivatives by Ioakimidis, Papadakis and Perdios (BIT 1991; 31 : 276–285), without applications to elasticity problems, but here the much more complicated case of the elastic complex potentials is studied even when just an appropriate non‐analytic complex density function (such as an edge dislocation/loading distribution density) is numerically available on the boundary. The present results are also directly applicable to inclusion problems, anisotropic elasticity, antiplane elasticity and classical two‐dimensional fluid dynamics, but, unfortunately, not to crack problems in fracture mechanics. Brief numerical results (for the complex potentials), showing the dramatic increase of the computational accuracy, are also displayed and few generalizations proposed. Copyright © 2000 John Wiley & Sons, Ltd.     

Time‐domain integration methods of exponentially damped linear systems

In this paper, three new kinds of time‐domain numerical methods of exponentially damped systems are presented, namely, the simplified Newmark integration method, the precise integration method, and the simplified complex mode superposition method. Based on the traditional Newmark integration method and transforming the equation of motion with exponentially damping kernel functions into an equivalent second‐order equation of motion by using the internal variables technique, the simplified Newmark integration method is developed by using a decoupling technique to reduce the computer run time and storage. By transforming the equation of motion with exponentially damping kernel functions into a first‐order state‐space equation, the precise integration technique is used to numerically solve the state‐space equation. Based on a symmetric state‐space equation and the complex mode superposition method, a delicate and simplified general solution of exponentially damped linear systems, completely in real‐value form, is developed. The accuracy and efficiency of the developed numerical methods are compared and discussed by two benchmark examples.     

A semi‐analytic collocation technique for steady‐state strongly nonlinear advection‐diffusion‐reaction equations with variable coefficients

The aim of this paper is to present a new semi‐analytic numerical method for strongly nonlinear steady‐state advection‐diffusion‐reaction equation (ADRE) in arbitrary 2‐D domains. The key idea of the method is the use of the basis functions which satisfy the homogeneous boundary conditions of the problem. Each basis function used in the algorithm is a sum of an analytic basis function and a special correcting function which is chosen to satisfy the homogeneous boundary conditions of the problem. The polynomials, trigonometric functions, conical radial basis functions, and the multiquadric radial basis functions are used in approximation of the ADRE. This allows us to seek an approximate solution in the analytic form which satisfies the boundary conditions of the initial problem with any choice of free parameters. As a result, we separate the approximation of the boundary conditions and the approximation of the ADRE inside the solution domain. The numerical examples confirm the high accuracy and efficiency of the proposed method in solving strongly nonlinear equations in an arbitrary domain.     

Numerical solution of Cauchy type singular integral equations with logarithmic weight,based on arbitrary collocation points

Singular integral equations with a Cauchy type kernel and a logarithmic weight function can be solved numerically by integrating them by a Gauss-type quadrature rule and, further, by reducing the resulting equation to a linear system by applying this equation at an appropriate number of collocation points x _k. Until now these x _k have been chosen as roots of special functions. In this paper, an appropriate modification of the original method permits the arbitrary choice of x _k without any loss in the accuracy. The performance of the method is examined by applying it to a numerical example and a plane crack problem.     

Dynamic modelling of the flat 2‐D crack by a semi‐analytic BIEM scheme

We present an efficient numerical method for solving indirect boundary integral equations that describe the dynamics of a flat two‐dimensional (2‐D) crack in all modes of fracture. The method is based on a piecewise‐constant interpolation, both in space and time, of the slip‐rate function, by which the original equation is reduced to a discrete convolution, in space and time, of the slip‐rate and a set of analytically obtained coefficients. If the time‐step interval is set sufficiently small with respect to the spatial grid size, the discrete equations decouple and can be solved explicitly. This semi‐analytic scheme can be extended to the calculation of the wave field off the crack plane. A necessary condition for the numerical stability of this scheme is investigated by way of an exhaustive set of trial runs for a kinematic problem. For the case investigated, our scheme is very stable for a fairly wide range of control parameters in modes III and I, whereas, in mode II, it is unstable except for some limited ranges of the parameters. The use of Peirce and Siebrits' ε‐scheme in time collocation is found helpful in stabilizing the numerical calculation. Our scheme also allows for variable time steps. Copyright © 2001 John Wiley & Sons, Ltd.     

Time‐dependent fractional advection–diffusion equations by an implicit MLS meshless method

Recently, many new applications in engineering and science are governed by a series of fractional partial differential equations (FPDEs). Unlike the normal partial differential equations (PDEs), the differential order in an FPDE is with a fractional order, which will lead to new challenges for numerical simulation, because most existing numerical simulation techniques are developed for the PDE with an integer differential order. The current dominant numerical method for FPDEs is finite difference method (FDM), which is usually difficult to handle a complex problem domain, and also difficult to use irregular nodal distribution. This paper aims to develop an implicit meshless approach based on the moving least squares (MLS) approximation for numerical simulation of fractional advection–diffusion equations (FADE), which is a typical FPDE The discrete system of equations is obtained by using the MLS meshless shape functions and the meshless strong‐forms. The stability and convergence related to the time discretization of this approach are then discussed and theoretically proven. Several numerical examples with different problem domains and different nodal distributions are used to validate and investigate the accuracy and efficiency of the newly developed meshless formulation. It is concluded that the present meshless formulation is very effective for the modeling and simulation of the FADE. Copyright © 2011 John Wiley & Sons, Ltd.     

Boundary knot method for some inverse problems associated with the Helmholtz equation

The boundary knot method is an inherently meshless, integration‐free, boundary‐type, radial basis function collocation technique for the solution of partial differential equations. In this paper, the method is applied to the solution of some inverse problems for the Helmholtz equation, including the highly ill‐posed Cauchy problem. Since the resulting matrix equation is badly ill‐conditioned, a regularized solution is obtained by employing truncated singular value decomposition, while the regularization parameter for the regularization method is provided by the L‐curve method. Numerical results are presented for both smooth and piecewise smooth geometry. The stability of the method with respect to the noise in the data is investigated by using simulated noisy data. The results show that the method is highly accurate, computationally efficient and stable, and can be a competitive alternative to existing methods for the numerical solution of the problems. Copyright © 2005 John Wiley & Sons, Ltd.