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.     

A starting algorithm for the numerical simulation of two-phase flow problems

The diferential equations governing two-phase flow through a porous medium are derived and the space discretization is carried out by the Galerkin finite element method. An examination of the eigenvalues of the resulting system follows in order to decide on a suitable stability criterion to be satisfied by the time-stepping algorithm. A class of time-stepping schemes is introduced and their behaviour when applied to a specific problem is analyzed. The aim is to achieve an adequate modelling of the water-oil interface and to preserve reasonable values of the capillary pressure. A technique to clustering the mesh around the interface is used to improved the modelling, and a boundary condition based on the relative mobilities of the two fiuids is imposed to enable the program to calculate the water-oil ratio after breakthrough. The effect of lumping the mass matrix is then examined with a view to reducting the surge in the values obtained for the capillary pressure. Finally an optimum method is suggested for use as a starting algorithm for more sophisticated time-stepping schemes.     

A lattice parameter method for the investigation of solid state precipitation

Over strictly limited composition ranges, the relationship between lattice parameter and solid solution composition can often be taken as linear even in multicomponent solid solution systems. The constants in the assumed linear equations relating the lattice parameter to the atomic percentage of each component may be calculated from lattice parameters measured for solution heat-treated alloys of known compositions. The subsequent changes in lattice parameters which occur during ageing of these alloys are shown to yield useful information about composition changes and the precipitate phases which occur during ageing. In particular, if the composition of the precipitate is known, then the linear equations may be used to calculate the composition of the remaining solid solution alloys. Even if both the composition of the precipitate phase and the remaining solid solution are unknown, changes in lattice parameter, combined with known lattice parameter versus composition functions, may be used to determine what average precipitate compositions are not allowed. These methods are illustrated in the case of Al-Cu, Cu-Ni-Al and Cu-Zn-Ni-Al alloys.     

An Efficient Harmonic Balance Method for Nonlinear Eddy-Current Problems

A method is presented which determines the steady-state solution of nonlinear eddy current problems. The unknown potentials are represented by Fourier-series and the nonlinear behavior of the material is split into a linear and a nonlinear term using a fixed-point technique. This approach leads to decoupled linear equations for each harmonic component. To take the nonlinearity into account, several fixed-point iterations have to be made. The method avoids calculating transient processes which normally have to be stepped through if using time-stepping methods. The present method is illustrated by two 2-D examples     

The flow of wet steam in a one‐dimensional nozzle

The formation of water droplets in a low‐pressure steam turbine, seriously degrades the efficiency of the generator. A model has been developed which includes the nucleation and subsequent growth of the droplets as the extra equations to the usual Euler equations for dry steam. A feature of this work is that all the equations are cast in Eulerian form compared to much of the previous work which considered the droplets in Lagrangian form. The ensuing equations are solved using a second‐order upwind TVD scheme which can cope with the steep gradients which occur in the solution. The results for a 1‐D nozzle are presented and compared with experimental results. Copyright © 1999 John Wiley & Sons, Ltd.     

Particular solutions for axisymmetric Helmholtz-type operators

In this paper, we consider the solution of the axisymmetric heat equation with axisymmetric data in an axisymmetric domain in R³. To solve this problem, we remove the time-dependence by various transform or time-stepping methods. This converts the problem to one of a sequence of modified inhomogeneous Helmholtz equations. Generalizing previous work, we consider solving these equations by boundary-type methods. In order to do this, one needs to subtract off a particular solution, so that one obtains a sequence of modified homogeneous Helmholtz equations. We do this by modifying the usual Dual Reciprocity Method (DRM) and approximating the right-hand sides by Fourier-polynomials or bivariate polynomials. This inevitably leads to analytical solving a sequence of ordinary differential equations (ODEs.) The analytic formulas and their precision are checked using mathematica. In fact, by using an infinite precision technique, the particular solutions can be obtained with infinite precision themselves. This work will form the basis for numerical algorithms for solving axisymmetric heat equation.     

ITERATIVE METHODS IN SOLVING NAVIER–STOKES EQUATIONS BY THE BOUNDARY ELEMENT METHOD

The solution of Navier–Stokes equations of time-dependent incompressible viscous fluid flow in planar geometry by the Boundary Domain Integral Method (BDIM) is discussed. The introduction of a subdomain technique to fluid flow problems is considered and improved in order to maintain the stability of BDIM. To avoid problems with flow kinematics computation in the sudomain mesh, a segmentation technique is proposed which combines the original BDIM with its subdomain variant and preserves its numerical stability. In order to reduce the computational cost of BDIM, which greatly depends on the solution of systems of linear equations, iterative methods are used. Conjugate gradient methods, conjugate gradients squared and an improved version of the biconjugate gradient method BiCGSTAB, together with the generalized minimal residual method, are used as iterative solvers. Different types of preconditioning, from simple Jacobi to incomplete LU factorization, are carried out and the performance of chosen iterative methods and preconditioners are reported. Test examples include backward facing step flow and flow through tubular heat exchangers. Test computation results show that BDIM is an accurate approximation technique which, together with the subdomain technique and powerful iterative solvers, can exhibit some significant savings in storage and CPU time requirements.     

A comparison of time marching schemes for the transient heat conduction equation

This paper investigates the phenomenon of ‘noise’ which is common in most time-dependent problems. The emphasis is on the achievable accuracy that is obtained with various time-stepping algorithms and how this can be improved if noise is artificially damped to an acceptable level. A series of experiments are made where the space domain is discretized using the finite element method and the variation with time is approximated by several finite difference methods. The conclusion is reached that the Crank–Nicolson scheme with a simple averaging process is superior to the other methods investigated.     

Elimination of Nonphysical Solutions and Implementation of Adaptive Step Size Algorithm in Time-Stepping Finite-Element Method for Magnetic Field–Circuit–Motion Coupled Problems

The time-stepping finite-element method (FEM) has become a powerful tool in solving transient electromagnetic fields. The formulation can include complex issues such as time harmonics and space harmonics, nonlinear magnetic property of iron materials, external circuit, and mechanical motion in the system equations. However, as the derivatives of physical quantities are usually unknown at the initial step of the time-stepping method, erroneous solutions might appear at the beginning of the transient process. To reduce the number of time steps, an adaptive step size algorithm can be used. In this paper, a method to eliminate the nonphysical or nonrealistic solutions at the start of the time-stepping finite-element analysis (FEA), when simulating the transient process of electric devices, is presented. A practical implementation of adaptive time step size algorithm for coupled problems is proposed. A matrix operation method, which can be understood clearly and implemented easily, that deals with matching boundary conditions in the study of mechanical motion, is also described.     

Time-domain transient elastodynamic analysis of 3-D solids by BEM

The numerical implementation of the Direct Boundary Element formulation for time-domain transient analysis of three-dimensional solids is presented in a most general and complete manner. The present formulation employs the space and time dependent fundamental solution (Stokes' solution) and Graffi's dynamic reciprocal theorem to derive the boundary integral equations in the time domain. A time-stepping scheme is then used to solve the boundary initial value problem by marching forward in time. Higher order shape functions are used to approximate the field quantities in space as well as in time, and a combination of analytical (time-integration) and numerical (spatial-integration) integration is carried out to form a system of linear equations. At the end of each time step, these equations are solved to obtain the unknown field quantities at that time. Finally, the accuracy and reliability of this algorithm is demonstrated by solving a number of example problems and comparing the results against the available analytical and numerical solution.     

A note on the grid movement induced by MFE

Moving-grid methods in one space dimension have become popular for solving several kinds of parabolic and hyperbolic partial differential equations (PDEs) involving fine scale structures such as steep moving fronts and emerging steep layers, pulses, shocks, etc. In two space dimensions, however, application of moving-grid methods is less trivial than in 1D. For some methods, e.g. those based on equidistributing principles, it is not even clear how to extend the underlying grid selection procedure to 2D. The moving-finite-element (MFE) method does not suffer from this drawback; its mathematical extension to 2D is trivial. However, because of the intrinsic coupling between the discretization of the PDE and the grid selection, the application of MFE, as for any other method, is not without difficulties. In this paper we describe the node movement induced by MFE for various PDEs and we indicate some problems concerning the grid structure that can result from this movement.     

A new moving finite element method based on quadratic approximation functions

Moving finite element methods are adaptive gridding procedures especially designed for systems of partial differential equations whose solutions contain steep gradients. A new moving finite element method based on quadratic approximation functions is presented. Both the theoretical and computational aspects are outlined. Performance of the method is illustrated with solutions to Burgers' equation. The solution is accurate and remarkably smooth in the entire domain.     

A new subregion boundary element technique based on the domain decomposition method

A new subregion boundary element technique based on the domain decomposition method is presented in this paper. This technique is applicable to the stress analysis of multi-region elastic media, such as layered-materials. The technique is more efficient than traditional methods because it significantly reduces the size of the final matrix. This is advantageous when a large number of elements need to be used, such as in crack analysis. Also, as the system of equations for each subregion is solved independently, parallel computing can be utilized. Further, if the boundary conditions are changed the only equations required to be recalculated are the ones related to the regions where the changes occur. This is very useful for cases where crack extension is modelled with new boundary elements or where crack faces come to contact. Numerical examples are presented to demonstrate the accuracy and efficiency of the method.     

A Legendre spectral element method for the Stefan problem

In this paper we present a Legendre spectral element method for solution of multi-dimensional unsteady change-of-phase Stefan problems. The spectral element method is a high-order (p-type) finite element technique, in which the computational domain is broken up into general (curved) quadrilateral macroelements, and the solution, data and geometry are expanded within each element in terms of tensor-product Lagrangian interpolants. The discrete equations are generated by a Galerkin formulation followed by Gauss–Lobatto Legendre quadrature, for which it is shown that exponential convergence to smooth solutions is obtained as the polynomial order of fixed elements is increased. The spectral element equations are inverted by conjugate gradient iteration, in which the matrix-vector products are calculated efficiently using tensor-product sum-factorization. To solve the Stefan problem numerically, the heat equations in the liquid and solid phases are transformed to fixed domains applying an interface-local time-dependent immobilization transformation technique. The modified heat equations are discretized using finite differences in time, resulting at each time step in a Helmholtz equation in space that is solved using Legendre spectral element elliptic discretizations. The new interface position is then computed using a variationally consistent flux treatment along the phase boundary, and the solution is projected upon the corresponding updated mesh. The rapid convergence rate and stability of the method are discussed, and numerical results are presented for a one-dimensional Stefan problem using both a semi-implicit and a fully implicit time-stepping scheme. Finally, a two-dimensional Stefan problem with a complex phase boundary is solved using the semi-implicit scheme.     

Fast Computation of Stationary Inviscid Flow Around an Airfoil

The parallel performance of an implicit solver for the Euler equations on a structured grid is discussed. The flow studied is a two-dimensional transonic flow around an airfoil. The spatial discretization involves the MUSCL scheme, a higher-order Total Variation Diminishing scheme. The solver described in this paper is an implicit solver that is based on quasi Newton iteration and approximate factorization to solve the linear system of equations resulting from the Euler Backward scheme. It is shown that the implicit time-stepping method can be used as a smoother to obtain an efficient and stable multigrid process. Also, the solver has good properties for parallelization comparable with explicit time-stepping schemes. To preserve data locality domain decomposition is applied to obtain a parallelizable code. Although the domain decomposition slightly affects the efficiency of the approximate factorization method with respect to the number of time steps required to attain the stationary solution, the results show that this hardly affects the performance for practical purposes. The accuracy with which the linear system of equations is solved is found to be an important parameter. Because the method is equally applicable for the Navier-Stokes equations and in three-dimensions, the presented combination of efficient parallel execution and implicit time-integration provides an interesting perspective for time-dependent problems in computational fluid dynamics.     

Application of polynomial preconditioners to conservation laws

Polynomial preconditioners which are suitable in implicit time-stepping methods for conservation laws are reviewed and analyzed. The preconditioners considered are either based on a truncation of a Neumann series or on Chebyshev polynomials for the inverse of the system-matrix. The latter class of preconditioner is optimal in a space of polynomials of certain degree if the matrix has only real eigenvalues and a non-singular system of eigenvectors. The preconditioning can be applied to any convergent splitting of the system matrix, i.e. to any classical implicit time-stepping method for conservation laws that is based on a quasi-Newton iteration. An efficient implementation based on SSOR is presented and the approach is applied to simulations of the viscous unsteady Burgers equation and to inviscid steady flow around an airfoil in two spatial dimensions to illustrate the method in large-scale computations. For viscous flows the efficiency increase due to preconditioning is considerable.     

Application of the Krein's method for determination of natural frequencies of periodically supported beam based on simplified Bresse-Timoshenko equations

Summary A new method is proposed for the stress analysis of an elastic space weakened by several arbitrarily located coplanar circular cracks under the action of an arbitrary normal pressure. The method is based on a new type of integral equations which has definite advantages over the existing methods: equations are non-singular, the iteration procedure is rapidly convergent even for very close interactions; there is no need to solve the integral equations if one is interested only in obtaining the upper and the lower bounds for the quantities of interest. In the case of the cracks far apart, these bounds are so close that provide, in fact, a sufficiently accurate solution to the problem. The method allows us also to obtain a practically exact numerical solution to the problem of very close interactions. Several illustrative examples are considered.     

Use of the over-relaxation technique in the simulation of large groundwater basins by the finite element method

The over-relaxation technique has been used in conjunction with the finite element method in a regional time-dependent simulation of the subsidence of Venice. One year of computer experiments have shown that in basins with a ratio h/R ? 10^?2 between the vertical and the horizontal dimension, the over-relaxation technique can lead to unsatisfactory results. In single precision and for relatively large time steps, the solution of the final linear system can be inaccurate if the optimum over-relaxation factor ω is not correctly assessed. In transient analyses, instability can also occur. The latter may be avoided by properly reducing the time step or by switching to double precision. Steady state simulations can also require double precision to provide accurate results even when the best ω is used. Instability and inaccuracy disappear in basins for which h/R ? 10^?1. In addition the processor time diminishes significantly as in this case the number of iterations necessary to obtain a good solution is considerably smaller than the order of the matrix.     

STABILITY ANALYSIS AND DESIGN OF TIME-STEPPING SCHEMES FOR GENERAL ELASTODYNAMIC BOUNDARY ELEMENT MODELS

In the literature there is growing evidence of instabilities in standard time-stepping schemes to solve boundary integral elastodynamic models. However, there has been no theory to support scientists and engineers in assessing the stability of their boundary element algorithms or to help them with the design of new, more stable algorithms. In this paper we present a general framework for the analysis of the stability of any time-domain boundary element model. We illustrate how the stability theory can be used to assess the stability of existing boundary element models and how the insight gained from this analysis can be used to design more stable time-stepping schemes. In particular, we describe a new time-stepping procedure that we have developed, which has substantially enhanced stability characteristics and greater accuracy for the same computational effort. The new scheme, which we have called ‘the half-step scheme’, is shown to have substantially improved performance for the displacement discontinuity boundary element method commonly used to model dynamic fracture interaction and propagation. © 1997 by John Wiley & Sons, Ltd.