An efficient method for solving implicit and explicit stiff differential equations

When the s‐stage fully implicit Runge–Kutta (RK) method is used to solve a system of n ordinary differential equations (ODE) the resulting algebraic system has a dimension ns. Its solution by Gauss elimination is expensive and requires 2s³n³/3 operations. In this paper we present an efficient algorithm, which differs from the traditional RK method. The formal procedure for uncoupling the algebraic system into a block‐diagonal matrix with s blocks of size n is derived for any s. Its solution is s²/2 times faster than the original, nondiagonalized system, for s even, and s³/(s−1) for s odd in terms of number of multiplications, as well as s² times in terms of number of additions/multiplications. In particular, for s=3 the method has the same precision and stability properties as the well‐known RK‐based RadauIIA quadrature of Ehle, implemented by Hairer and Wanner in RADAU5 algorithm. Unlike RADAU5, however, the method is applicable with any s and not only to the explicit ODEs My′=f(x, y), where M=const., but also to the general implicit ODEs of the form f(x, y, y′)=0. The block‐diagonal form of the algebraic system allows parallel processing. The algorithm formally differs from the implicit RK methods in that the solution for y is assumed to have a form of the algebraic polynomial whose coefficients are found by enforcing y to satisfy the differential equation at the collocation points. Locations of those points are found from the derived stability function such as to guarantee either A‐ or L‐stability properties as well as a superior precision of the algorithm. If constructed such as to be L‐stable the method is a good candidate for solving differential‐algebraic equations (DAEs). Although not limited to any specific field, the application of the method is illustrated by its implementation in the multibody dynamics described by both ODEs and DAEs. Copyright © 2000 John Wiley & Sons, Ltd.     

Weight‐adjusted discontinuous Galerkin methods: Matrix‐valued weights and elastic wave propagation in heterogeneous media

Weight‐adjusted inner products are easily invertible approximations to weighted L² inner products. These approximations can be paired with a discontinuous Galerkin (DG) discretization to produce a time‐domain method for wave propagation which is low storage, energy stable, and high‐order accurate for arbitrary heterogeneous media and curvilinear meshes. In this work, we extend weight‐adjusted DG methods to the case of matrix‐valued weights, with the linear elastic wave equation as an application. We present a DG formulation of the symmetric form of the linear elastic wave equation, with upwind‐like dissipation incorporated through simple penalty fluxes. A semidiscrete convergence analysis is given, and numerical results confirm the stability and high‐order accuracy of weight‐adjusted DG for several problems in elastic wave propagation.     

A two‐grid method for expanded mixed finite‐element solution of semilinear reaction–diffusion equations

We present a scheme for solving two‐dimensional semilinear reaction–diffusion equations using an expanded mixed finite element method. To linearize the mixed‐method equations, we use a two‐grid algorithm based on the Newton iteration method. The solution of a non‐linear system on the fine space is reduced to the solution of two small (one linear and one non‐linear) systems on the coarse space and a linear system on the fine space. It is shown that the coarse grid can be much coarser than the fine grid and achieve asymptotically optimal approximation as long as the mesh sizes satisfy H=O(h^1/3). As a result, solving such a large class of non‐linear equation will not be much more difficult than solving one single linearized equation. Copyright © 2003 John Wiley & Sons, Ltd.     

Convergence of the fixed point theorem with random parameters

Solving stochastic non‐linear dynamical problems represents a formidable task which, in many cases, can be achieved solely through numerical simulation techniques. This is true for high dimensional as well as low dimensional problems. One method to deal with the non‐linearity is to use the fixed point theorem which gives the convergence conditions of the iterative scheme towards the equilibrium point of the equation. In this paper we look at the particular case where the equilibrium equation depends on a random variable. This case arises for instance in the study of coupled non‐linear dynamical systems when structural uncertainties are introduced in the dynamical systems. We give almost sure and L ^p convergence conditions for the simulation iterative scheme. Copyright © 2003 John Wiley & Sons, Ltd.     

Solution of low péclet number heat transfer with viscous dissipation in the semi‐infinite axial region of a tube via functional analytic methods

Abstract

A solution of the extended Graetz problem with prescribed wall heat flux and viscous dissipation in a semi‐infinite axial region of a tube is obtained by functional analytic methods. The energy equation is split into a set of partial differential equations to obtain a self‐adjoint formulism. Then, an algebraic characteristic equation of the eigenvalue problem for an arbitrary velocity profile is obtained by an approximation method in L ₂[0, 1]. In addition, a backward recursive formula for calculating the expansion coefficients of the solution is developed.     

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.     

On the L2 and the H1 couplings for an overlapping domain decomposition method using Lagrange multipliers

In this paper, a comparison of the L² and the H¹ couplings is made for an overlapping domain decomposition method using Lagrange multipliers. The analysis of the local equations arising from the formulation of the coupling of two mechanical models shows that continuous weight functions are required for the L² coupling term whereas both discontinuous and continuous weight functions can be used for the H¹ coupling. The choice of the Lagrange multiplier space is discussed and numerically studied. The paper ends with some numerical examples of an end‐loaded cantilever beam and a cracked plate under tension and shear. It is shown that the continuity enforced with the H¹ coupling leads to a link with a flexibility that can be beneficial for coupling a very coarse mesh with a very fine one. To limit the effect of the volume coupling on the global response, a narrow coupling zone is recommended. In this case, volume coupling tends to a surface coupling, especially with a L² coupling. Copyright © 2006 John Wiley & Sons, Ltd.     

Refined global–local higher‐order theory and finite element for laminated plates

Based on completely three‐dimensional elasticity theory, a refined global–local higher‐order theory is presented as enhanced version of the classical global–local theory proposed by Li and Liu (Int. J. Numer. Meth. Engng. 1997; 40 :1197–1212), in which the effect of transverse normal deformation is enhanced. Compared with the previous higher‐order theory, the refined theory offers some valuable improvements these are able to predict accurately response of laminated plates subjected to thermal loading of uniform temperature. However, the previous higher‐order theory will encounter difficulty for this problem. A refined three‐noded triangular element satisfied the requirement of C¹ weak‐continuity conditions in the inter‐element is also presented. The results of numerical examples of moderately thick laminated plates and even thick plates with span/thickness ratios L/h = 2 are given to show that in‐plane stresses and transverse shear stresses can be reasonably predicted by the direct constitutive equation approach without smooth technique. In order to accurately obtain transverse normal stresses, the local equilibrium equation approach in one element is employed herein. Copyright © 2006 John Wiley & Sons, Ltd.     

Impact force on downstream bed of weir by free overfall flow

Abstract

According to previous experimental studies on free overfall, the main influencing factors of the local scour downstream of a drop structure were found to be the impact position of the nappe flow and the force striking on the channel bed. In this study, the pressure transducers and strain amplifiers, which did not disturb the flow field, were set up to measure the pressure distributions along the streamwise direction downstream of the weir. The data were recorded through an A‐D converter and compared with the data collected by Moore (1943) and Rand (1955). Expressions for the impact position and impact force of the flow nappe derived from the experimental results and the theoretical analysis. The relationship among weir height h, impact position of free overfall L_p , and drop number D were put into the following equation: L_p/h = 1.82D ^0.17.     

Different approaches for mixed LSFEMs in hyperelasticity: Application of logarithmic deformation measures

We present geometrically nonlinear formulations based on a mixed least‐squares finite element method. The L²‐norm minimization of the residuals of the given first‐order system of differential equations leads to a functional, which is a two‐field formulation dependent on displacements and stresses. Based thereon, we discuss and investigate two mixed formulations. Both approaches make use of the fact that the stress symmetry condition is not fulfilled a priori due to the row‐wise stress approximation with vector‐valued functions belonging to a Raviart‐Thomas space, which guarantees a conforming discretization of H(div). In general, the advantages of using the least‐squares finite element method lie, for example, in an a posteriori error estimator without additional costs or in the fact that the choice of the polynomial interpolation order is not restricted by the Ladyzhenskaya‐Babu?ka‐Brezzi condition (inf‐sup condition). We apply a hyperelastic material model with logarithmic deformation measures and investigate various benchmark problems, adaptive mesh refinement, computational costs, and accuracy.     

A fourth‐order compact scheme for the Helmholtz equation: Alpha‐interpolation of FEM and FDM stencils

We propose a fourth‐order compact scheme on structured meshes for the Helmholtz equation given by R(φ):=f( x )+Δφ+ξ²φ=0. The scheme consists of taking the alpha‐interpolation of the Galerkin finite element method and the classical central finite difference method. In 1D, this scheme is identical to the alpha‐interpolation method (J. Comput. Appl. Math. 1982; 8 (1):15–19) and in 2D making the choice α=0.5 we recover the generalized fourth‐order compact Padé approximation (J. Comput. Phys. 1995; 119 :252–270; Comput. Meth. Appl. Mech. Engrg 1998; 163 :343–358) (therein using the parameter γ=2). We follow (SIAM Rev. 2000; 42 (3):451–484; Comput. Meth. Appl. Mech. Engrg 1995; 128 :325–359) for the analysis of this scheme and its performance on square meshes is compared with that of the quasi‐stabilized FEM (Comput. Meth. Appl. Mech. Engrg 1995; 128 :325–359). In particular, we show that the relative phase error of the numerical solution and the local truncation error of this scheme for plane wave solutions diminish at the rate O((ξ?)⁴), where ξ, ? represent the wavenumber and the mesh size, respectively. An expression for the parameter α is given that minimizes the maximum relative phase error in a sense that will be explained in Section 4.5. Convergence studies of the error in the L² norm, the H¹ semi‐norm and the l^∞ Euclidean norm are done and the pollution effect is found to be small. Copyright © 2010 John Wiley & Sons, Ltd.     

Lattice Dislocations Enhancing Thermoelectric PbTe in Addition to Band Convergence

Phonon scattering by nanostructures and point defects has become the primary strategy for minimizing the lattice thermal conductivity (κ_L) in thermoelectric materials. However, these scatterers are only effective at the extremes of the phonon spectrum. Recently, it has been demonstrated that dislocations are effective at scattering the remaining mid‐frequency phonons as well. In this work, by varying the concentration of Na in Pb_0.97Eu_0.03Te, it has been determined that the dominant microstructural features are point defects, lattice dislocations, and nanostructure interfaces. This study reveals that dense lattice dislocations (≈4 × 10¹² cm^?2) are particularly effective at reducing κ_L. When the dislocation concentration is maximized, one of the lowest κ_L values reported for PbTe is achieved. Furthermore, due to the band convergence of the alloyed 3% mol. EuTe the electronic performance is enhanced, and a high thermoelectric figure of merit, zT , of ≈2.2 is achieved. This work not only demonstrates the effectiveness of dense lattice dislocations as a means of lowering κ_L, but also the importance of engineering both thermal and electronic transport simultaneously when designing high‐performance thermoelectrics.     

A new quadrilateral area coordinate method (QACM‐II) for developing quadrilateral finite element models

The quadrilateral area coordinate method proposed in 1999 (hereinafter referred to as QACM‐I) is a new and efficient tool for developing robust quadrilateral finite element models. However, such a coordinate system contains four components (L₁, L₂, L₃, L₄), which may make the element formulae and their construction procedure relatively complicated. In this paper, a new category of the quadrilateral area coordinate method (hereinafter referred to as QACM‐II), containing only two components Z₁ and Z₂, is systematically established. This new coordinate system (QACM‐II) not only has a simpler form but also retains the most important advantages of the previous system (QACM‐I). Hence, as an application, QACM‐II is used to formulate a new 4‐node membrane element with internal parameters. The whole process is similar to that of the famous Wilson's Q6 element. Numerical results show that the present element, denoted as QACII6, exhibits much better performance than that of Q6 in benchmark problems, especially for MacNeal's thin beam problem. This demonstrates that QACM‐II is a powerful tool for constructing high‐performance quadrilateral finite element models. Copyright © 2007 John Wiley & Sons, Ltd.     

On continuous,discontinuous, mixed,and primal hybrid finite element methods for second‐order elliptic problems

Finite element formulations for second‐order elliptic problems, including the classic H¹‐conforming Galerkin method, dual mixed methods, a discontinuous Galerkin method, and two primal hybrid methods, are implemented and numerically compared on accuracy and computational performance. Excepting the discontinuous Galerkin formulation, all the other formulations allow static condensation at the element level, aiming at reducing the size of the global system of equations. For a three‐dimensional test problem with smooth solution, the simulations are performed with h‐refinement, for hexahedral and tetrahedral meshes, and uniform polynomial degree distribution up to four. For a singular two‐dimensional problem, the results are for approximation spaces based on given sets of hp‐refined quadrilateral and triangular meshes adapted to an internal layer. The different formulations are compared in terms of L²‐convergence rates of the approximation errors for the solution and its gradient, number of degrees of freedom, both with and without static condensation. Some insights into the required computational effort for each simulation are also given.     

Fourier‐based numerical approximation of the Weertman equation for moving dislocations

This work addresses the numerical approximation of solutions to a dimensionless form of the Weertman equation, which models a steadily moving dislocation and is an important extension (with advection term) of the celebrated Peierls‐Nabarro equation for a static dislocation. It belongs to the class of nonlinear reaction‐advection‐diffusion integro‐differential equations with Cauchy‐type kernel, thus involving an integration over an unbounded domain. In the Weertman problem, the unknowns are the shape of the core of the dislocation and the dislocation velocity. The proposed numerical method rests on a time‐dependent formulation that admits the Weertman equation as its long‐time limit. Key features are (1) time iterations are conducted by means of a new, robust, and inexpensive Preconditioned Collocation Scheme in the Fourier domain, which allows for explicit time evolution but amounts to implicit time integration, thus allowing for large time steps; (2) as the integration over the unbounded domain induces a solution with slowly decaying tails of important influence on the overall dislocation shape, the action of the operators at play is evaluated with exact asymptotic estimates of the tails, combined with discrete Fourier transform operations on a finite computational box of size L; (3) a specific device is developed to compute the moving solution in a comoving frame, to minimize the effects of the finite‐box approximation. Applications illustrate the efficiency of the approach for different types of nonlinearities, with systematic assessment of numerical errors. Converged numerical results are found insensitive to the time step, and scaling laws for the combined dependence of the numerical error with respect to L and to the spatial step size are obtained. The method proves fast and accurate and could be applied to a wide variety of equations with moving fronts as solutions; notably, Weertman‐type equations with the Cauchy‐type kernel replaced by a fractional Laplacian.     

A finite element semi‐Lagrangian method with L2 interpolation

High‐order accurate methods for convection‐dominated problems have the potential to reduce the computational effort required for a given order of solution accuracy. The state of the art in this field is more advanced for Eulerian methods than for semi‐Lagrangian (SLAG) methods. In this paper, we introduce a new SLAG method that is based on combining the modified method of characteristics with a high‐order interpolating procedure. The method employs the finite element method on triangular meshes for the spatial discretization. An L² interpolation procedure is developed by tracking the feet of the characteristic lines from the integration nodes. Numerical results are illustrated for a linear advection–diffusion equation with known analytical solution and for the viscous Burgers’ equation. The computed results support our expectations for a robust and highly accurate finite element SLAG method. Copyright © 2012 John Wiley & Sons, Ltd.     

Application of multiple assessment items optimization method to determining the repair ranking of existing reinforced concrete bridge system

Abstract

The D. (Degree) E. (Extent) R. (Relevancy) evaluation method is widely used to assess the damage states of existing reinforced concrete (RC) bridges in Taiwan. The present study is first to distinguish between relevancy (R^a, a = 1.5, 2) and non‐relevancy (R^a, a= 1) for the repair ranking of the existing RC bridge system to be assessed. The multiple assessment items optimization method was mainly applied to seek the optimum repair ranking. The five existing RC bridges, i.e. Chi‐jou, Lan‐yang, Dah‐jea, Dah‐duh and Dah‐an in Taiwan are chosen as practical examples. The results show that when each bridge is judged to be non‐relevant by the system, the repair ranking predicted by the D.E.R. evaluation method is correct. Nevertheless, when the bridge system has one or more one important bridge, the repair ranking predicted by the D.E.R. evaluation method is not accurate. The proposed method may be remedied the defect of the D.E.R. evaluation method to predict the repair rankings of existing RC bridge system.     

Optimal modal reduction of dynamic subsystems: Extensions and improvements

The problem of optimal reduction of a linear dynamic subsystem is revisited. The subsystem may represent, for example, a complex but minor part of a large elastic structure. The goal is to drastically reduce the number of degrees of freedom of a subsystem attached through an interface to a main system in such a way as to affect the dynamic behavior of the main system in the least possible way. In a recent publication, the Optimal Modal Reduction (OMR) algorithm was developed to this end. This algorithm seeks a reduction of the subsystem that will have minimal effect, in the L₂ norm, on the Dirichlet‐to‐Neumann (DtN) map on the interface. Here this algorithm is extended and improved in a number of ways. First, a family of alternative formulations are derived for the DtN map, which lead to alternative OMR algorithms; one of them, called OMR _j=2, is shown to yield better results than the original formulation. Second, the OMR algorithm, which was originally developed for undamped subsystems, is extended to subsystems undergoing Rayleigh damping. In addition, an enhanced derivation of the original OMR algorithm is presented, and the good pointwise performance of OMR is explained by relating to minimization with respect to the H¹ norm. The extensions mentioned above are discussed theoretically as well as demonstrated via numerical examples. Experiments and discussion include comparison of the OMR algorithm to simple coarsening of the subsystem discretization. In all cases, central finite difference discretization in space and explicit time‐stepping are employed to solve the scalar wave equation. Copyright © 2010 John Wiley & Sons, Ltd.     

Evolution of an N-level system via automated vectorization of the Liouville equations and application to optically controlled polarization rotation

The Liouville equation governing the evolution of the density matrix for an atomic/molecular system is expressed in terms of a commutator between the density matrix and the Hamiltonian, along with terms that account for decay and redistribution. To find solutions of this equation, it is convenient first to reformulate the Liouville equation by defining a vector corresponding to the elements of the density operator, and determining the corresponding time-evolution matrix. For a system of N energy levels, the size of the evolution matrix is N^2?×^?N². When N is very large, evaluating the elements of these matrices becomes very cumbersome. We describe a novel algorithm that can produce the evolution matrix in an automated fashion for an arbitrary value of N. As a non-trivial example, we apply this algorithm to a 15-level atomic system used for producing optically controlled polarization rotation. We also point out how such a code can be extended for use in an atomic system with arbitrary number of energy levels.     

The effect of cooling load and thermal conductance on the local stability of an endoreversible refrigerator

The purpose of this paper is to present a local stability analysis of an endoreversible refrigerator operating at the minimum input power P for given cooling load R absorbed from the cold reservoir, for different thermal conductances α and β, in the isothermal couplings of the working fluid with the heat reservoirs T_H and T_L (T_H>T_L). An endoreversible refrigerator system that is modeled by the differential equation may depend on the numerical values of certain parameters that appear in the equation. From the local stability analysis we find that a critical point of an almost linear system is a stable node. After a small perturbation the system state exponentially decays to steady state with either of two relaxation times that are a function of α, β, T_L, R and the heat capacity C. We can exhibit qualitatively the behavior of solutions of the system by sketching its phase portrait. One eigenvector in a phase portrait is the nonzero constant vector, and the other is a function of α, β, R, T_H and T_L. Finally, we discuss the local stability and energetic properties of the endoreversible refrigerator.